Development of a 5-plex SILAC Method Tuned for the Quantitation of Tyrosine Phosphorylation Dynamics

The propagation of phosphorylation downstream of receptor tyrosine kinases is a key dynamic cellular event involved in signal transduction, which is often deregulated in disease states such as cancer. Probing phosphorylation dynamics is therefore crucial for understanding receptor tyrosine kinases' function and finding ways to inhibit their effects. MS methods combined with metabolic labeling such as stable isotope labeling with amino acids in cell culture (SILAC) have already proven successful in deciphering temporal phosphotyrosine perturbations. However, they are limited in terms of multiplexing, and they also are time consuming, because several experiments need to be performed separately. Here, we introduce an innovative approach based on 5-plex SILAC that allows monitoring of phosphotyrosine signaling perturbations induced by a drug treatment in one single experiment. Using this new labeling strategy specifically tailored for phosphotyrosines, it was possible to generate the time profiles for 318 unique phosphopeptides belonging to 215 proteins from an erlotinib-treated breast cancer cell line model. Hierarchical clustering of the time profiles followed by pathway enrichment analysis highlighted epidermal growth factor receptor (EGFR or ErbB1) and ErbB2 signaling as the major pathways affected by erlotinib, thereby validating the method. Moreover, based on the similarity of its time profile to those of other proteins in the ErbB pathways, the phosphorylation at Tyr453 of protein FAM59A, a recently described adaptor of EGFR, was confirmed as tightly involved in the signaling cascade. The present investigation also demonstrates the remote effect of EGFR inhibition on ErbB3 phosphorylation sites such as Tyr1289 and Tyr1328, as well as a potential feedback effect on Tyr877 of ErbB2. Overall, the 5-plex SILAC is a straightforward approach that extends sample multiplexing and builds up the arsenal of methods for tyrosine phosphorylation dynamics.

The propagation of phosphorylation downstream of receptor tyrosine kinases is a key dynamic cellular event involved in signal transduction, which is often deregulated in disease states such as cancer. Probing phosphorylation dynamics is therefore crucial for understanding receptor tyrosine kinases' function and finding ways to inhibit their effects. MS methods combined with metabolic labeling such as stable isotope labeling with amino acids in cell culture (SILAC) have already proven successful in deciphering temporal phosphotyrosine perturbations. However, they are limited in terms of multiplexing, and they also are time consuming, because several experiments need to be performed separately. Here, we introduce an innovative approach based on 5-plex SILAC that allows monitoring of phosphotyrosine signaling perturbations induced by a drug treatment in one single experiment. Using this new labeling strategy specifically tailored for phosphotyrosines, it was possible to generate the time profiles for 318 unique phosphopeptides belonging to 215 proteins from an erlotinib-treated breast cancer cell line model. Hierarchical clustering of the time profiles followed by pathway enrichment analysis highlighted epidermal growth factor receptor (EGFR or ErbB1) and ErbB2 signaling as the major pathways affected by erlotinib, thereby validating the method. Moreover, based on the similarity of its time profile to those of other proteins in the ErbB pathways, the phosphorylation at Tyr453 of protein FAM59A, a recently described adaptor of EGFR, was confirmed as tightly involved in the signaling cascade. The present investigation also demonstrates the remote effect of EGFR inhibition on ErbB3 phosphorylation sites such as Tyr1289 and Tyr1328, as well as a potential feedback effect on Tyr877 of ErbB2. Overall, the 5-plex SILAC is a straightforward approach that extends sample multiplexing and builds up the arsenal of methods for tyrosine phosphorylation dynamics. Molecular & Cellular Proteomics 12 The transmission of information in a cell is largely driven by protein phosphorylation (1,2). It is used to convey a signal within the cell and modulate different biological outcomes such as transcription regulation or cytoskeleton rearrangement. The receptor tyrosine kinases represent key players at the plasma membrane interface because they operate as initiators of phosphorylation cascades (3,4). Among the 58 receptor tyrosine kinases encoded by the human genome, the ErbB family members-composed of ErbB1 (also called epidermal growth factor receptor (EGFR)), 1 ErbB2, ErbB3, and ErbB4 -play a pivotal role because they control important processes such as cell differentiation, proliferation, and apoptosis (5). Deregulation of ErbB signaling networks is known to be connected to various diseases such as cancer (6,7). Thus, the ErbB family members represent an important target class for understanding and curing disease-indeed, both small-molecule tyrosine kinase inhibitors (TKIs) and monoclonal antibodies are under extensive development (8). In order to evaluate the proper effect of drugs on signaling pathways and better recognize their mode of action, one must assess tyrosine phosphorylation changes of ErbBs and their associated downstream proteins.
Recently, MS-based quantitative phosphoproteomics approaches have developed into attractive and powerful tools for studying signaling networks (9 -12). Moreover, it has become apparent that investigations based on "static" binary comparisons (control versus treated) need to be replaced by alternative ones in which the dynamics of the networks can be monitored (13)(14)(15)(16). Also, peptide-centric strategies are preferred because they permit deeper sampling of the phosphoproteome and discriminate between differentially regulated sites on a single protein (15). Typically, phosphoproteome dynamics screens are performed using one of two approaches: stable isotope labeling with amino acids in cell culture (SILAC) using three different versions of labeled lysine and arginine (3-plex SILAC) (17), or isobaric tags for relative and absolute quantitation (iTRAQ) using four (4-plex iTRAQ) (18) and eight (8-plex iTRAQ) (19) labeling reagents. MS signal quantitation is then performed using the precursor (MS1based) and reporter ions (MS2-based) for SILAC and iTRAQ, respectively. Although both strategies are highly valuable, they also have some shortcomings. The 3-plex SILAC approach requires the combination of two separate parallel experiments, including a common reference for normalization, to obtain a five-time-point kinetic. This implies that the two samples also need to be processed and measured individually, augmenting the analysis time per condition. The iTRAQ approach constricts the labeling and mixing steps at a late stage of the sample preparation procedure, introducing variability issues. Moreover, during data acquisition, co-fragmentation of interfering peptides present in the isolation window of the precursor ions perturbs the quantitative information encoded in the reporter ion region (20 -22). It has also been reported that iTRAQ derivatized phosphopeptides tend to undergo charge enhancement during ionization, leading to a significant drop in identification efficiency (23).
Here, we describe the development of a new 5-plex SILACbased strategy combining serial peptide immunoaffinity purification and LC-MS/MS detection for the examination of tyrosine phosphorylation dynamics. The approach relies on robust MS1-based quantitation and offers the advantage of increased sample multiplexing. The model breast cancer cell line KPL-4 was differentially labeled into five populations using a mixture of light and heavy isotope-encoded lysines, arginines, and tyrosines. After sample combination, each detected tyrosine-phosphorylated tryptic peptide was present as a quintuplet in MS, which allowed the monitoring of five time points in a single experiment. This method was applied to investigate the effects of the selective EGFR TKI erlotinib (Tarceva) on ErbB signaling over a 1-h treatment. It was possible to identify, quantitate, and generate the temporal drug response profiles of 318 unique tyrosine-phosphorylated peptides belonging to 215 proteins. After hierarchical clustering and pathway enrichment analysis, the data clearly show EGFR and ErbB2 as the main affected signaling networks.

5-plex SILAC Cell Culturing, Treatment, and Digest Preparation-
KPL-4 cells were grown as five different populations for at least six doublings in DMEM supplemented with 10% dialyzed fetal bovine serum and glutamine-containing antibiotics (Pen Strep Glutamine, Invitrogen) and with the following amino acid combinations: SILAC medium 1 (SM1): Lys0, Arg0, Tyr0; SM2: Lys4, Arg6, Tyr0; SM3: Lys8, Arg10, Tyr0; SM4: Lys8, Arg10, Tyr6; SM5: Lys8, Arg10, Tyr10 (see Fig. 1). Plated cells (2 ϫ 10 7 ) were serum starved for 24 h, and SM2-5 were treated with erlotinib (5 M end concentration) for 5, 10, 30, and 60 min, respectively. SM1 was treated with dimethyl sulfoxide to serve as the control. Cells were washed twice with cold PBS before the addition of lysis buffer for 10 min at room temperature. Lysis buffer consisted of 6 M urea, 2 M thiourea, 200 mM ammonium bicarbonate pH 8.5, and phosphatase inhibitors (PhosStop, Roche). Lysed cells were collected by scraping the plates, homogenized for 3 ϫ 10 s on ice using a tip sonifier (Branson), and centrifuged at 12,000 ϫ g for 10 min. A 10-l aliquot was taken from each sample for protein quantitation using a modified Bradford assay. After mixing in equal portions, 20 mg of total proteins were reduced with 5 mM DTT for 1 h and alkylated with 55 mM iodoacetamide for 1 h at room temperature. The sample was diluted 8-fold with 50 mM ammonium bicarbonate pH 8.5, 200 g of trypsin (Roche) was added, and proteins were digested overnight at 37°C. A supplementary 200 g of trypsin was added, and digestion continued for 4 h at 37°C. After acidification with formic acid at 5% v/v end concentration, samples were concentrated onto C18 Sep-Pak cartridges (100 mg, Waters, Etten-Leur, The Netherlands), and eluted peptides were dried down using a SpeedVac. Three independent biological replicates were prepared at one-week intervals following the same procedure.
Immunoprecipitation-Tyrosine phosphorylated peptides were sequentially immunoprecipitated with two antibodies following the procedure described by Boersema et al. (26), with some modifications. Briefly, samples were re-suspended in immunoprecipitation (IP) buffer (50 mM Tris-HCl pH 7.2, 150 mM NaCl, and 1 mM PMSF) and incubated with pre-equilibrated PY99-agarose beads (Santa Cruz Biotechnology, Santa Cruz, CA) overnight at 4°C. Following the first PY99 IP, the same samples were further incubated with pre-equilibrated 4G10-agarose beads (Millipore, Billerica, MA) for 4 h at 4°C. Both types of beads were washed three times with IP buffer and twice with water before elution with 0.15% TFA for 15 min at room temperature. Eluted peptide fractions from PY99 and 4G10 IPs were concentrated on stop-and-go extraction tips, dried down, and stored at Ϫ20°C. Separate PY99 and 4G10 samples were re-suspended in 5% formic acid/2% acetonitrile and analyzed in duplicate via LC-MS.
Online Nanoflow Liquid Chromatography Mass Spectrometry-Nanoflow LC-MS/MS was performed by coupling an Easy-nLC system (Proxeon, Odensa, Denmark) to an LTQ Orbitrap Velos (Thermo Fisher Scientific, Bremen, Germany) equipped with a nanoelectrospray source (Proxeon) and an active background-ion reduction device (ABIRD, ESI Source Solutions, Woburn, MA). Peptides were loaded onto an Aqua C18 (Phenomenex, Torrance, CA) trapping column (packed in-house, 10 mm length ϫ 100 m inner diameter, 5 m particle size) and separated on a Reprosil-Pur C18-AQ (Dr. Maisch GmbH, Ammerbuch, Germany) analytical column (200 mm length ϫ 75 m inner diameter, 3 m particle size, 120 Å). Trapping of the peptides was performed using the intelligent flow control at a maximal backpressure of 250 bar with 100% solvent A (0.6% acetic acid), and elution was carried out at a flow rate of 250 nl/min using the following gradient of solvent B (acetonitrile and 0.6% acetic acid): 0%-2% B in 2 min, 2%-35% B in 106 min, 35%-80% in 2 min, and 80% B for 10 min. The eluent was electrosprayed via a normal coated emitter tip butt-connected to the analytical column using a PicoClear union (New Objective, Cambridge, MA). The mass spectrometer was operated in the data-dependent mode to automatically switch be-tween MS and MS/MS. Survey full-scan MS spectra were acquired from m/z 350 to 1700 in the Orbitrap with a resolution of r ϭ 60,000 (at m/z 400) after accumulation to a target value of 1 ϫ 10 6 in the linear ion trap for a maximum time of 500 ms. The detected ions were recalibrated on the fly using the ambient air polysiloxane at m/z 445.120024 as lock-mass (27). The 10 most intense ions with a charge state Ͼ1ϩ and a threshold above 1000 were selected for collision-induced dissociation in the linear ion trap. Ions were accumulated to a target value of 10,000 using the predictive ion trap fill time mode, isolated with a width of 2 amu, and collision-induced dissociation was performed at a normalized collision energy of 35% with an activation Q of 0.25 for 10 ms and wideband activation mode enabled. Fragmented ions with an m/z width of Ϯ10 ppm were placed on an exclusion list for 30 s and a maximal size of 500.
LC-MS Data Processing-Raw files belonging to a biological replicate were combined for multi-dimensional protein identification technology and processed using Proteome Discoverer version 1.2 (Thermo Fisher Scientific), comprising a self-designed peptide/protein identification and quantitation workflow (all details can be found in supplemental Fig. S1). Data were searched using Mascot version 2.3.2 (Matrix Science, London, UK) against the human UniProt/Swis-sProt database (release 2011 04, 20,232 sequences) with trypsin as an enzyme, allowing a maximum of two miscleavages and 20 ppm and 0.5 Da as precursor and fragment mass tolerances, respectively. Carbamidomethylated cysteine was selected as a static modification, and the following were set as dynamic modifications: oxidized methionine; phosphorylated serine, threonine, and tyrosine, with the latter additionally in the Tyr6 labeled (ϩ85.986460 Da and ϩ85.9358 Da as monoisotopic and average mass shifts, respectively) and Tyr10 labeled (89.993559 and 89.9072 as monoisotopic and average mass, respectively) forms; Lys4 and Lys8 labeled lysine; and Arg6 and Arg10 labeled arginine. Protein grouping was enabled, and leucine and isoleucine were considered as equal. Individual top-ranked peptideto-spectrum matches (PSMs) were accepted for a Mascot score of Ն20 (26). PSMs were filtered by keeping only sequences containing a tyrosine residue and a phosphorylation, and the resulting initial dataset was used to address the identification efficiency of the method (see "Evaluation of the 5-plex SILAC Dataset"). The false discovery rate at this score was experimentally tested via a concatenated target-decoy database search that yielded Ͻ1% for phosphotyrosine containing peptides for the three biological replicates (supplemental Fig. S2). A common number of natural isotopes for all members of a quintuplet were used to create extracted ion chromatograms, compute the peak areas, and determine relative peptide ratios. To generate the extracted ion chromatograms, the precursor ion mass precision between successive scans was set to 2 ppm, the minimal signal-to-noise value was set at 1, and quintuplets were assembled allowing a maximum retention time deviation of 0.2 min. PSMs with no quantitation value (approximately 34% to 40% of PSMs) were removed. Quintuplet ratios from tyrosine containing unphosphorylated peptides nonspecifically co-immunoprecipitated in the experiment were used for data normalization with each separate biological replicate. Phosphopeptide quintuplet ratios originating from different biological replicates were log2-averaged, standard deviations and statistical significance using a paired two-tailed Student's t test were calculated in Excel, and only the peptide with the highest Mascot score when multiple isotopologues were identified was reported. Peptides bearing the same phosphorylation site and showing identical time profiles but with additional oxidized methionine and/or miscleavages were considered as redundant and were excluded. The false localization rate was tested using the Mascot delta score between two successive identical peptide sequences differing only in phosphorylation site assignment and accepted for values greater than or equal to 6 (see "Discussion") (28).
Data Clustering and Analysis-The phosphotyrosine time profiles were clustered hierarchically using Euclidean distance and Ward linkage. Applying a distance threshold of 0.32 resulted in nine time profile clusters determined to be optimal for subsequent pathway enrichment analysis. The correlation between a cluster's proteins and signaling pathways was assessed as follows: enrichment statistical significance of proteins belonging to each individual cluster (and combinations of adjacent clusters) in gene sets available from the Reactome pathway database (version 40, March 2012) was calculated via Fisher's exact test with the significance threshold set to 0.01. Gene sets were pre-filtered to consider only candidate proteins (N) complemented with the rest of the dataset (215-N "background" proteins), and at least three proteins per enriched pathway were required.

RESULTS
Overview of the 5-plex SILAC Approach-Focusing on the phosphotyrosine fraction of the phosphoproteome enabled us to utilize Tyr as an additional amino acid for protein labeling (29), along with the more conventional Lys and Arg. By selecting a proper combination of light and heavy labeled Lys, Arg, and Tyr, we designed a new 5-plex SILAC strategy to distinguish tyrosine-containing peptides that originate from five different cell populations. The cell culture media were prepared with a mixture of Lys0/Arg0/Tyr0, Lys4/Arg6/Tyr0, Lys8/Arg10/Tyr0, Lys8/Arg10/Tyr6, or Lys8/Arg10/Tyr10 and numbered SM1 to SM5, respectively (Fig. 1). As a consequence, after digestion of the proteins with trypsin, every tyrosine-containing peptide holding a C-terminal Lys was detected in the MS with mass shifts of 4, 8, 14, and 18 Da (supplemental Fig. S3). When mixed in equal amounts, peptides of interest were present as SILAC quintuplets. The selection of the heavy labeled amino acids was done in such a way that a minimal mass shift of 4 Da between two peptide isotopologues was guaranteed. This constraint was of importance for minimizing the overlap between the isotope clusters within a quintuplet. A similar picture was observed when the peptide C terminus was an Arg, with mass shifts of 6, 10, 16, and 20 Da, respectively.
Four of the five KPL-4 cell populations were subjected to the selective EGFR TKI erlotinib for 5, 10, 30, and 60 min, and the last one was treated with dimethyl sulfoxide to serve as a control. In the present study, SM1-containing only light amino acids-served as the control, and SM2-5 were successively used for the longer treatment times, so that the "heaviest medium" represented the end point of the kinetic. After the lysates had been mixed in equal proportions and the proteins had been trypsinized, tyrosine phosphorylated peptides were enriched from the crude digest via the application of a serial anti-pTyr IP procedure. That is, peptides were initially incubated with anti-phosphotyrosine (PY99) affinity resin, and the unbound peptides were then additionally treated with the second anti-phosphotyrosine (4G10) affinity resin. Finally, both phosphotyrosine enriched peptide samples were analyzed separately via nano-LC-MS/MS on an Orbitrap Velos.
An illustration of the approach is shown at the peptide level with the quintuplet matched to the doubly charged peptide VTIADDpYSDPFDAK from the SH2 domain-containing adapter protein B (SHB) and phosphorylated on Tyr246 (Fig. 2). The corresponding ion was recorded at m/z 818.8421 for the isotopologue originating from SM1. Along with the unlabeled form, the co-eluting SM2-5 isotopologues at t R ϭ 53.6 min were found at m/z 820.8543 (ϩ4 Da), 822.8491 (ϩ8 Da), 825.8594 (ϩ14 Da), and 827.8630 (ϩ18 Da), respectively ( Fig. 2A). Because of the mass shifts observed within the quintuplet, the candidate phosphopeptide could be identified as containing a C-terminal lysine. Despite the complex nature of the sample, all five isotopologue precursor ions were selected for fragmentation, which demonstrates that the MS duty cycle was short enough to handle 5-plex SILAC labeled samples. Moreover, for this particular candidate, all five MS/MS events resulted in successful peptide identification when submitted to a database search (Fig. 2C). Relative quantitation was achieved by comparing quintuplets' precursor ion intensities after extracted ion chromatograms, followed by data normalization and averaging across biological replicates. For Tyr246 of SHB, erlotinib showed only a moderate down-regulative effect, with the minimum of the time profile observed after 30 min (approximately a 2-fold de-crease) followed by a slow return to the basal level after 60 min (Fig. 2D).
Evaluation of the 5-plex SILAC Dataset-Three independent biological replicates, each including two IPs that were both analyzed in duplicate via LC-MS, were performed (see Fig. 1). Searching the 12 raw files against a human protein database (see "Experimental Procedures" for the details) resulted in 6163 PSMs passing the Ն20 Mascot score criterion for phosphotyrosine-immunoprecipitated peptides (26). In order to get an idea of the number of identified isotopologues per quintuplet and how deep the peptide sequencing for a 5-plex SILAC experiment is, the PSMs were combined, and distinct phosphopeptides (i.e. same phosphorylation site with different versions of light or heavy Lys, Arg, and Tyr) were counted (supplemental Table S1). A total number of 846 phosphopeptides that mapped to 340 proteins was obtained. It is of importance to emphasize that this preliminary "uncurated" dataset was intended for method evaluation only, as it still contained redundancy. A consistent fraction (412 peptides representing 49% of the total) was obtained from the fragmentation of a single isotopologue (Fig. 3). Focusing on this subset, 260 (63%) were "unlabeled" SM1, whereas 37 (9%), 80 (19%), 12 (3%), and 23 (6%) were from labeled SM2-5 precursor ions, respectively. Also, examination of precursor ion masses indicated that the heavier the peptide labels were, the greater the bias toward the identification of low-m/z candidates was (supplemental Fig. S4). A possible explanation for this is the unavoidable increased overlap between isotopologues for large peptides leading to mixed (or "chimeric") MS/MS (30). Overall, about half of the identified phosphopeptides were related to a single isotopologue form, whereas as few as 7% had all five SM1-5 "versions" resulting in peptide identification.
In order to reduce the PSMs down to a list of unique phosphotyrosine peptides, the aforementioned candidate list was manually curated. For instance, entries containing no quantitation value, additional methionine oxidation, or miscleavage were removed. However, a phosphopeptide presenting an additional phosphorylation site was considered as non-redundant, and only the candidate with the highest Mascot score was kept when found in several biological replicates. Furthermore, phosphorylation site localization was tested using the Mascot delta score and considered confident for values greater than or equal to 6, based on the fact that the antibodies used are specific for phosphorylated tyrosine, rather than serine/threonine (28). This resulted in a collection of 318 unique identified and quantified phosphotyrosine peptides grouped into 215 proteins (Fig. 4, supplemental Table S2). All annotated MS/MS together with their time profiles are provided in supplemental Fig. S5. When the reproducibility of 5-plex SILAC was taken into account, 98 peptides (31%) were identified in all biological replicates, and 187 peptides (59%) were present in at least two replicates (Fig. 4A). It has to be mentioned that the 318 peptides in this final dataset were detected only as quintuplets in the MS-that is, each time    were removed during the manual curation process. The overlap between the three biological replicates was 31% (98 peptides in common), and around 42% to 50% between two of them, similar to what has been reported by others (12,13). Because two antibodies were used for serial IPs, data were examined to highlight whether a difference in IP efficiency was observed. Substantially more unique peptides (229 versus 89) originated from the 4G10 than from the PY99 agarose-conjugate fractions (see Fig. 4B). This further emphasizes that combining several pTyr peptide affinity enrichment steps allows deeper sample coverage (12,31).
Tyr-phosphorylated Sites Modulated by Erlotinib-Erlotinib is a reversible inhibitor of the EGFR kinase that competes with ATP for its binding site. Therefore, the phosphorylation of EGFR, its substrates, and downstream molecules is expected to be down-regulated (32)(33)(34). For the present experiment, an erlotinib concentration of 5 M was selected because it corresponds approximately to the measured plasma level in humans (35). Also, the last point of the investigated time window was selected as 60 min because limited EGFR degradation should be observed within this specific period of time (36). Hierarchical clustering of the 318 pTyr peptide time profiles revealed nine groups and readily allowed the recognition of the most affected phosphosites (Fig. 5, supplemental Table  S2). For instance, clusters 1 and 2, containing three and nine time profiles, respectively, exhibited an immediate, strong, and sustained decrease of phosphorylation. In the case of cluster 1 members, the reduction was on the order of 10-fold after only 10 min of treatment and concerned Tyr187 of MAPK1 and Tyr204 of MAPK3 (also known as ERK2 and ERK1, respectively), two markers of the Ras-Raf-Mek pathway connecting EGFR to MAPK/Erk signaling. A similar effect was reported for erlotinib-treated HN5 squamous carcinoma cells (32). Furthermore, a third member belonging to cluster 1, namely, Tyr659 of the GRB2-associated binding protein 1 (GAB1), was found among the most affected phosphorylation sites. The latter is a known substrate of EGFR (37) and is involved in EGF-mediated activation of ERK2 via binding to the protein-tyrosine phosphatase PTPN11 (alternatively named SHP2) (38). Also, a strong down-regulation of GAB1 Tyr659 phosphorylation was observed in H3255 non-small cell lung cancer cells treated with gefitinib, another selective EGFR TKI (39).
Time profiles found in cluster 2 displayed a similar rapid decrease after 10 min, but with a down-regulation magnitude of about 6-fold (Fig. 5B). Two EGFR autophosphorylation sites, Tyr1110 and Tyr1172, were present in this cluster and indicated a reduction of receptor kinase activity. These two

FIG. 4. Identified and quantitated Tyr-phosphorylated peptides.
A, overlap between the three biological replicates. Peptides unique to a biological replicate (keeping the highest Mascot score candidate) are indicated between brackets. B, distribution among agarose-conjugate resins used for peptide pTyr IPs. Left-hand and right-hand bars represent all and unique peptides, respectively. EGFR tyrosine residues, along with Tyr1092 and Tyr1197 (found in cluster 4 in the present study), were inhibited by erlotinib in a concentration-dependent manner (33). Other highly down-regulated phosphosites were Tyr771 of phospholipase C-gamma-1 (PLCG1) and Tyr427 of SHC-transforming protein 1 (SHC1). The latter two phosphosites have been described as affected by TKI treatment (32,39) and are known to specifically bind to phosphorylated Tyr1172 and Tyr1197 of EGFR via their SH2 binding domain (40 -42). Another interesting phosphorylation site located in cluster 2 appeared to be Tyr453 of protein FAM59A. This protein (also called GAREM) is a recently discovered EGFR adaptor that interacts with GRB2 and SHP2 and participates in ERK activation (43).
At the other end of the dendogram, a large proportion of the time profiles (68% (215 out of 318)) populate clusters 7-9 (Fig.  5A). These particular phosphorylation sites showed little inhibitory effect, while other sites appeared to increase upon treatment with erlotinib. In particular, phosphorylation of Tyr877 of ErbB2 (member of cluster 9) became gradually up-regulated and showed a 2-fold increase in signal after 60 min. This particular site, located in the activation loop of ErbB2, is known to be a substrate of the Src family kinases (c-Src, Fyn, and Yes) and to modulate the receptor's catalytic activity (44 -46). Up-regulated Tyr-phosphorylation sites for c-Src (Tyr530) and Fyn (Tyr185, Tyr213, Tyr420, and Tyr440) were highlighted in the present experiment, and their time profiles were found in clusters 8 and 9 (see also Fig. 6).
Signaling Pathways Affected by Erlotinib-To better understand whether there was a connection between proteins found in a specific cluster and the erlotinib effect, we performed pathway enrichment analysis. For this, the Reactome pathway database was used and the over-representation of members of each cluster, as well as of two adjacent clusters, was calculated (see "Experimental Procedures"). As a result, when members of clusters 1 and 2 were combined, corresponding to a total of 13 candidate proteins, 6 proteinsnamely, EGFR, GAB1, PLCG1, MAPK1, MAPK3, and SHC1were mapped to 51 pathways (supplemental Table S3). In particular, EGFR and ErbB2 signaling were found to be the top enriched pathways. Proteins identified in the present experiment, together with their time profiles, are represented schematically for the EGFR and ErbB2 signaling pathways (Fig. 6). Other proteins containing phosphotyrosine sites that were strongly affected by erlotinib, such as Tyr453 FAM59A and Tyr1135 INPPL1, but which were not previously described as members of EGFR/ErbB2 signaling were included. Our data confirm that SHC1, PLCG1, and GAB1 are central molecules involved in transmission of the signal that regulates the MAPK/Erk signaling pathway, and they therefore validate the 5-plex SILAC approach as a valuable method for monitoring phosphotyrosine dynamics.
In addition to the aforementioned proteins, we examined the behavior of other phosphorylation sites in the ErbB family. Because KPL-4 cells highly express ErbB2 and ErbB3 (no ErbB4 expression) (25), we looked at the crosstalk between EGFR and ErbB2-ErbB3, especially knowing that ErbB2 needs no ligand for activation and ErbB3 cannot signal by itself (47). Grouping time profiles for each individual ErbB family member enabled the identification of ErbB3 phosphorylation sites that were systematically affected by erlotinib relative to ErbB2 (supplemental Fig. S6). That is, four ErbB3 phosphotyrosine sites, Tyr1159, Tyr1289, Tyr1307, and Tyr1328, all located in the C-terminal tail, displayed a sus-tained down-regulated time profile (clusters 4 -6). To better visualize the remote effect of erlotinib on ErbB2 and ErbB3, sequences of the C termini were aligned (supplemental Fig.  S6C). The same outcome had already been reported for ErbB3 Tyr1159 and Tyr1328 in the case of dual EGFR/ErbB2 TKI lapatinib-treated BT-474 breast cancer cells (46). Because ErbB3 phosphorylation down-regulation is an indirect result of EGFR inhibition, the observed levels were less than that observed for EGFR, which was in the range of 3-fold reduction. Interestingly, detected ErbB3 downstream interacting partners, such as GRB7 (Tyr284, Tyr388) or the phosphatidylinositol 3-kinase (PI3K) members PIK3R2 (Tyr464) and PIK3R3 (Tyr199, Tyr341), also exhibited decreased phosphorylation over time (supplemental Fig. S6B). The latter could FIG. 6. Representative Tyr phosphorylation sites affected by erlotinib in the EGFR and ErbB2 canonical pathways. The proteins are colored according to their association with the different pathways. Other known undetected members are indicated to complement the pathways. Even though ErbB3 is reported as an ErbB2 signaling pathway member, it is colored gray to avoid confusion. Only EGFR-ErbB2 and EGFR-ErbB3 heterodimers are represented for clarity. Arrows representing selected known interactions between proteins (black) and ErbB's phosphosites and adapter (blue) are indicated. point to an inhibitory effect on the PI3K/Akt signaling pathway (32), which was further tested by immunoblot detection of phospho-Akt (Thr308 and Ser473). However, the connections between EGFR inhibition and PI3K/Akt signaling via phospho-ErbB3 need to be investigated in more detail. DISCUSSION The examination of phosphorylation perturbations caused by a drug treatment has evolved from single to multiple time experiments. To allow this, recent developments in MS-based quantitation methods have focused on higher multiplexing (48). We have established an innovative 5-plex SILAC approach involving peptide IP to monitor Tyr-phosphorylation dynamics and compare five different cellular conditions within a single experiment. In addition to Lys and Arg, stable isotopically labeled Tyr ( 13 C 6 and [ 13 C 9 , 15 N 1 ]) was combined to introduce the necessary peptide mass shifts. A similar 5-plex SILAC method was proposed by Molina et al. to monitor protein expression changes during adipocyte differentiation; however, their approach was limited to the detection of Argcontaining peptides (49). Also, the choice of a minimum mass shift of 4 Da in the method described here results in MSseparated isotopologues, which was considered to deliver more reliable quantitation data than strategies involving overlapping isotopologues and necessitating deconvolution. Our protocol has the advantage of combining all five SILAC labeled samples at the beginning of the workflow and thus keeps variability between experimental conditions at low levels. Also, the reduction to a single sample for monitoring five conditions represents a clear advantage in terms of the shorter analysis time. Because MS signal quantitation is performed at the precursor ion level, the method does not suffer from the weaknesses inherent to isobaric labeling, such as reporter ion signal distortion (due to co-eluting peptide fragmentation) (20 -22) and peptide charge enhancement (23).
Using three biological replicates and 20 mg (5 ϫ 4 mg) total starting material, it was possible to identify 846 peptides, quantitate 318 unique phosphotyrosine sites found on 215 proteins, and reveal their dynamics over 60 min of erlotinib treatment. These numbers are in good agreement with what has previously been reported in the literature for temporal (15,50) and dose-response (51) studies involving peptide anti-pTyr IP. As a typical example, Wolf-Yadlin et al. were able to quantify the dynamics of 222 phosphotyrosines found on 146 proteins from EGF-stimulated 184A1 human mammary epithelial cells (52). In another study, Boersema et al. obtained the dynamics of 73 unique phosphotyrosines (52 proteins) from EGF-treated HeLa cells, albeit using as little as 2 mg of starting material (26). Also, the reproducibility of the 5-plex SILAC approach was similar to that typically achieved by other methods (12,52), as the overlap represented approximately 30% between two replicates (see Fig. 4A).
Despite the a priori rather complex nature of the 5-plex SILAC encoded samples, more than half the peptides (i.e. 434 out of 846; see Fig. 3A) were identified as Ն2 isotopologue forms. From a quantitation point of view, it is obvious that the identification of only one member of a quintuplet is sufficient for subsequent generation of its tyrosine phosphorylation time profile. As anticipated, the reduction of complexity owing to the IP fractionation resulted in a sample containing a number of peptide quintuplets that is manageable in MS. A certain bias toward the identification of SM1 peptides was observed, presumably due to the higher abundance of untreated control phosphotyrosines relative to their SM2-5 erlotinib down-regulated isotopologues. To better address and circumvent this effect, SM and treatment times could have been randomized among the replicates, but this was outside the scope of the present study. It has to be mentioned that our final phosphopeptide dataset was composed of quintuplets only; that is, no peptide was observed with less than five isotopologues due to complete signal absence. For instance, Tyr204 of MAPK3, one of the most affected phosphotyrosines, was still recorded with minute signals for the SM3-5 isotopologues used for quantitation (supplemental Fig. S7). It is accepted that anti-pTyr pan-specific antibodies possess an intrinsic specificity (12,31). We chose to have two serial IP steps in our analytical workflow to maximize phosphotyrosine coverage; that is, first we precipitated the lysates with PY99, and then we followed with an incubation of the unbound peptide fraction with 4G10agarose conjugates. In the present experiment, 89 of the unique identified peptides were from PY99, and 229 (72%) were related to 4G10 (see Fig. 4B). It should be emphasized that this reflects better performance of the latter antibody, because the PSMs with the highest Mascot scores were kept to populate the final list of 318 unique phosphopeptides. Furthermore, the observation should be taken with care, as the two IPs were performed in serial and not in parallel.
The decision to apply the newly introduced 5-plex SILAC approach to EGFR signaling was motivated by the fact that the latter is well documented (see Refs. 12,15,16, and 40 for a review) and could serve as good model for evaluation. Also, using KPL-4 as the model cell line allowed us to probe the effect of EGFR inhibition on the other two family members ErbB2 and ErbB3. By performing phosphorylation time profile clustering followed by pathway enrichment analysis, we highlighted the involvement of EGFR and ErbB2 signaling as the main pathways affected by erlotinib and thereby confirmed the validity of our approach. The analysis unfortunately did not reveal any additional significant pathways perturbed by erlotinib, as the set of 215 genes was quite small. Despite the fact that many phosphosites identified in the present study were already known interactors of these pathways, additional ones appeared to be also linked, as suggested by their phosphorylation dynamics. For instance, Tyr453 of FAM59A, possessing a time profile similar to those of Tyr1110 and Tyr1172 of EGFR (all three are found in cluster 2; see "Results"), is a confirmed candidate of the EGFR signaling pathway. The effect of erlotinib on other ErbB family members, in particular the inactive kinase ErbB3, was also quantified. Based on the time profiles and their location within the C-terminal tail of ErbB3, Tyr1289, Tyr1307, and Tyr1328 seem to act as alternative docking sites for adaptor proteins. Interestingly, based on the sequence alignments of the three receptor tyrosine kinases, phosphorylated Tyr1307 of ErbB3 might have a distinct role, such as affinity for a specific adaptor, because neither EGFR nor ErbB2 possesses a homologous tyrosine (see supplemental Fig. S6C).
The 5-plex SILAC approach developed here is a valid alternative to existing methods. It has the clear advantage of allowing the analysis of five different conditions within a single sample. It joins the robustness of precursor-ion-based quantitation with higher multiplexing than in conventional iTRAQor SILAC-based strategies. Use of the 5-plex SILAC to study EGFR inhibition resulted in the elucidation of phosphotyrosine time profiles specifically affected by a treatment. This information is highly valuable, as it can lead to better understanding of the drug's mode of action and the finding of potential new biomarkers. Future applications of the strategy will be focused on the study of more complex therapies, such as combined small molecule-antibody or antibody-antibody treatments. Finally, the reported 5-plex SILAC approach has the potential to be further improved to increase its overall performance in phosphotyrosine dynamics quantification. For instance, analysis of the samples with a newer generation of MS, such as Q-Exactive, would improve the analytical depth, as thousands of phosphotyrosine sites are now routinely amenable to quantitation. Furthermore, data processing using software tailored for 5-plex SILAC including feature-rather than identification-centric capability would potentially also increase phosphopeptide sampling.