Research ArticleCell Cycle

Quantitative Phosphoproteomics Reveals Widespread Full Phosphorylation Site Occupancy During Mitosis

+ See all authors and affiliations

Sci. Signal.  12 Jan 2010:
Vol. 3, Issue 104, pp. ra3
DOI: 10.1126/scisignal.2000475

Abstract

Eukaryotic cells replicate by a complex series of evolutionarily conserved events that are tightly regulated at defined stages of the cell division cycle. Progression through this cycle involves a large number of dedicated protein complexes and signaling pathways, and deregulation of this process is implicated in tumorigenesis. We applied high-resolution mass spectrometry–based proteomics to investigate the proteome and phosphoproteome of the human cell cycle on a global scale and quantified 6027 proteins and 20,443 unique phosphorylation sites and their dynamics. Co-regulated proteins and phosphorylation sites were grouped according to their cell cycle kinetics and compared to publicly available messenger RNA microarray data. Most detected phosphorylation sites and more than 20% of all quantified proteins showed substantial regulation, mainly in mitotic cells. Kinase-motif analysis revealed global activation during S phase of the DNA damage response network, which was mediated by phosphorylation by ATM or ATR or DNA-dependent protein kinases. We determined site-specific stoichiometry of more than 5000 sites and found that most of the up-regulated sites phosphorylated by cyclin-dependent kinase 1 (CDK1) or CDK2 were almost fully phosphorylated in mitotic cells. In particular, nuclear proteins and proteins involved in regulating metabolic processes have high phosphorylation site occupancy in mitosis. This suggests that these proteins may be inactivated by phosphorylation in mitotic cells.

Introduction

The cell cycle is a highly regulated and evolutionarily conserved process that results in the duplication of the cell’s content and the equal distribution of the duplicated chromosomes into a pair of daughter nuclei. Progression through the cell cycle is governed by a complex network of regulatory proteins and signaling pathways, and deregulation of key players that coordinate this process can lead to carcinogenesis (1, 2). Global analyses of transcriptome dynamics during the cell cycle have been performed in microarray studies (3, 4). However, in addition to changes in messenger RNA (mRNA) abundance, protein phosphorylation and targeted protein degradation are also important regulators of cell cycle progression. Protein phosphorylation and dephosphorylation is a highly controlled biochemical process that responds to various intracellular and extracellular stimuli. Phosphorylation status modulates protein activity, influencing the tertiary and quaternary structure of a protein, controlling subcellular distribution, and regulating interactions with other proteins. Because phosphorylation is a dynamic process, quantitation of these phosphorylation events is crucial. Advances in phosphoproteomics, including phosphopeptide-enrichment methods, high-accuracy mass spectrometry (MS), and associated bioinformatic tools, make it feasible to obtain an unbiased view of phosphoproteomes at a “systems level” (58). Regulatory protein phosphorylation is a transient modification that is often of low occupancy or “stoichiometry,” meaning that only a fraction of a particular protein may be phosphorylated on a given site at any particular time, and that occurs on regulatory proteins of low abundance, such as protein kinases and transcription factors. Therefore, specific and efficient phosphopeptide-enrichment methods along with sensitive and accurate mass-spectrometric instrumentation are useful for analyzing complex protein mixtures, such as complete cell lysates.

We perform large-scale analysis of a phosphoproteome, which involves peptide fractionation by strong-cation exchange chromatography followed by phosphopeptide enrichment by titanium dioxide (TiO2) or immobilized metal affinity chromatography (IMAC) (9). The enriched phosphopeptide fractions are separated by online nanoscale reversed-phase high-pressure liquid chromatography (HPLC) directly coupled to a high-performance mass spectrometer. Phosphopeptides are ionized and analyzed by high-resolution MS using orbitrap mass analyzers in combination with fast sequencing by phosphorylation-specific multistage MS in a linear ion trap. Previously, combining all these methodologies with a quantitative method that relies on stable isotope labeling of amino acids in cell culture (SILAC) (10), we have obtained a global view of dynamic regulation of phosphorylation in mammalian cells as a function of a growth factor [epidermal growth factor (EGF)] stimulus in a time-resolved manner (11). Further developments in sample preparation strategies, high-resolution MS, computational proteomics, and associated bioinformatics have facilitated global proteome quantitation between different cellular states (1214). Here, we applied these phosphoproteomic and proteomic technologies to the cellular phosphorylation events and obtained a systems view of protein and phosphorylation-site dynamics during the cell cycle of HeLa S3 cells, a human cultured cell line. We furthermore introduce a method to determine the occupancy of thousands of phosphorylation sites, with which we found that a substantial proportion was nearly fully phosphorylated.

Results

Quantitative proteomics and phosphoproteomics to monitor changes over the cell cycle

Three HeLa S3 cell populations were labeled with light, medium, and heavy stable isotopic versions of arginine and lysine for SILAC analysis. Light and heavy populations were synchronized in six different stages of the cell cycle and mixed with medium-labeled, asynchronously growing cells that served as a common reference to allow comparison between all samples. HeLa S3 cells were synchronized in G1-S phase by a single thymidine block. One of the cell populations was harvested at this point, whereas the other cells were released into fresh medium and harvested at the indicated time points (Fig. 1A). A population of cells was released from the thymidine block into nocodazole-containing medium to enrich for mitotic cells. Subsequently, cells were released from the nocodazole block into fresh medium for 30 min to allow complete reformation of the bipolar spindle and, thus, to more closely resemble physiological mitotic conditions. The resulting mitotic phosphorylation events may not completely mimic M phase; therefore, we cannot rule out that at least some of the phosphorylation sites may be the result of spindle checkpoint activation that has not been completely reversed. Thymidine and nocodazole block and release methods are standard synchronization protocols that provide high proportions of synchronized cells and adequately reproduce cell cycle events and that have been used in other quantitative phosphoproteomics analyses (15, 16). We monitored synchronization of SILAC cell cultures by Western blotting and fluorescence-activated cell sorting (FACS) analyses (Fig. 1, B and C). After harvesting, cells from the three experiments were combined and processed as described in Materials and Methods for proteome and phosphoproteome analyses.

Fig. 1

Quantitative proteomic analysis of the human cell cycle. (A) HeLa S3 cells were SILAC-labeled with three different isotopic forms of arginine and lysine. Three individual populations of heavy and light SILAC cells were synchronized with a thymidine block and then collected at six different time points across the cell cycle after release from the thymidine arrest. Two samples were collected after a cell cycle arrest with nocodazole and release. Cells were lysed and mixed in equal amounts using an asynchronously growing cell population (medium SILAC) as the internal standard to allow normalization between experiments. Three independent experiments were performed to cover six cell cycle stages. (B) Immunoblot analysis of known cell cycle marker proteins in the different cell populations. (C) FACS profiles of the individual synchronized HeLa S3 populations. Cells were fixed and collected by centrifugation after which the DNA content of the cells was determined with propidium iodide. (D) Representative MS data showing how the abundance of the proteins was monitored in three experiments (Exp. 1, Exp. 2, Exp. 3) to obtain information from the six stages of the cell cycle. The data show the MS analysis of a tryptic SILAC peptide triplet derived from the cell cycle marker protein Geminin. Relative peptide abundance changes were normalized to the medium SILAC peptide derived from the asynchronously grown cells in all three experiments. The inset shows the combined six-time profile of Geminin over the cell cycle.

Using the MaxQuant software (13), applying unified statistical criteria, we analyzed liquid chromatography–tandem mass spectrometry (LC-MS/MS) data from 144 gel slices for the quantitative proteome and 333 fractions for the phosphoproteome. At a false discovery rate (FDR) of 1% on the protein level, this resulted in the identification of 6695 proteins (table S1), from which 6027 quantitative cell cycle profiles were obtained (Fig. 1D). The proteome identified here substantially overlaps with the proteins in two previous large-scale SILAC studies of the HeLa cell proteome (13, 17) and has more than two-thirds overlap with the proteome measured with another method (18) (fig. S1A). A total of 24,714 phosphorylation events were identified (FDR <1%, first applied at the peptide and then at the protein level), 20,443 of which were assigned to a specific residue with high confidence, with median localization confidence greater than 99.5% [class I sites (11)] (table S2). Comparison of these sites with three large-scale HeLa phosphoproteomes revealed that almost 10,000 sites are previously unknown (fig. S1B); 78% overlap with our previous time-resolved EGF phosphoproteome study (11) and 51% overlap with a study of phosphorylation sites on samples enriched for HeLa kinases (15).We identified 60% of the phosphorylation sites that were reported in a quantitative phosphoproteomics analysis of two stages of the HeLa cell cycle (16) and an additional 12,722 phosphorylation sites that were not identified in that study. The cell cycle phosphoproteome is by far the largest phosphoproteome that we or others have reported for any process (fig. S1B). Previously, 30% of total cellular proteins were estimated to be phosphorylated (19), but our data show that this proportion is at least 70% (5192 phosphoproteins in our proteome of 6695 proteins). This percentage could even increase as the phosphoproteome and proteome are probed in greater depth and transient phosphorylation sites on low-abundance proteins are captured more efficiently. The distribution of phosphoproteins across cellular compartments was essentially as observed before (11), with significant overrepresentation of nuclear proteins and underrepresentation of mitochondrial and secreted proteins (fig. S1C).

Regulatory proteins were well represented in our dynamic HeLa cell proteome. We found 220 protein kinases, 3% of the detected proteins. Because 2.5% of all genes encode kinases (table S3), this low-abundance protein class was not underrepresented in our proteome. Many macromolecular complexes were completely present and typically more than 70% of the members of low abundance cellular machineries were also quantified as estimated by the detection of 50 of 64 proteins in the “golden set” used in (20). The HeLa proteome contains at least 8827 distinct proteins (fig. S1A), which still does not encompass all known cell cycle factors (they cover 86% of the golden set). Therefore, the size of this cell line’s proteome is likely to be at least 10,000 proteins of which we quantified a majority. This is the largest quantified proteome to date, similar in size to transcriptome profiles of human cell lines (fig. S2).

The quality of the proteome-wide quantitation is demonstrated by the consistency of the dynamic MS-based profiles for the marker proteins (Fig. 1B) and key cell cycle–regulated proteins (fig. S3) with their abundance as shown by Western blotting. For example, we detected changes in the abundance of cyclins during the cell cycle (Fig. 1B and table S1), which ensures the oscillation in cyclin-dependent kinase (CDK) activity that is required for the global cell cycle control system. Furthermore, the abundance of a fifth of the proteome changed by at least fourfold throughout the cell cycle (difference between lowest and highest abundance). Because a fourfold change also best accounted for the dynamics of already described cell cycle components, we used this as a threshold for subsequent analysis (figs. S4 and S5). An RNA interference screen identified a cell cycle phenotype for more than 1000 proteins (21). Together, these findings implicate a large number of proteins in mediating the cell cycle, either by phenotype or by substantial changes in the abundance of proteins that make up the proteome or in their phosphorylation, many more than the known central cell cycle actors.

Analyzing the cell cycle–regulated proteins across the six time points revealed distinct up- and down-regulated clusters corresponding to each cell cycle stage (Fig. 2A). To understand the functional differences between the distinct cell cycle–regulated protein clusters, we assessed each protein cluster separately for overrepresented pathways, biological processes, and cellular components with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway using the proteomic phenotyping enrichment analysis previously described (22) (Materials and Methods). We retained each functional protein category that reached at least 95% statistical significance in one of the clusters and then performed one-way unsupervised clustering of the P values of the resulting categories. The most prominent cluster of proteins that were more abundant in mitosis than in other phases of the cell cycle was enriched in proteins involved in cell division and encompassed categories such as mitotic chromosome condensation, spindle organization and biogenesis, as well as mitotic cell cycle (fig. S6). The cellular compartments most overrepresented in proteins that were abundant during S phase were nucleoplasmic and nucleolar proteins (fig. S7). KEGG pathway mapping of proliferation- and growth-dependent signaling pathways, such as the insulin and mTOR (mammalian target of rapamycin) pathways, showed that they are activated in G1 and S phase cells (fig. S8).

Fig. 2

Dynamics of the proteome during the cell cycle. (A) Proteins whose abundance changed at least fourfold during the cell cycle were clustered in all cell cycle stages by calculating a time peak index by weighted mean of the ratio of maximal abundance. For each cell cycle stage, there are clear patterns of up- and down-regulation. (B) A circularized representation of the data shown in (A) was used to determine the angle in the cell cycle where the abundance of particular proteins peaks. Coordinately regulated protein complexes and organellar proteins at each cell cycle stage are indicated around the circle. (C) Comparison of mRNA and protein dynamics during the cell cycle. Measured protein dynamics were correlated to published mRNA data (4). Proteins were grouped on the y axis in four categories from top to bottom: unchanging mRNA and protein, changing mRNA and unchanging protein, unchanging mRNA and changing protein, and changing mRNA and changing protein. The x axis shows clustered gene ontology (GO) biological process terms enriched in at least one of the above four categories. High and low represent statistical over- or underrepresentation, respectively.

To determine when specific proteins and protein groups peaked in the cell cycle, we circularized the profiles of the grouped proteins on a cell cycle timeline (Fig. 2B). We observed statistically significant co-regulation of subunits within cell cycle–regulated complexes, on the basis of directional statistics using the Rayleigh test (Materials and Methods). Examples include the anaphase promoting complex/cyclosome (APC/C) (P < 0.02), mediator complex (P < 1 × 10−5), and DNA replication factor A complex (P < 0.007). Likewise, we found that proteins associated with particular subcellular organelles were significantly co-regulated (Fig. 2B), for example, mitochondrial proteins (P < 1.6 × 10−25), nucleolar proteins (P < 0.01), and components of the endoplasmic reticulum and Golgi (P < 0.001).

Comparison with microarray data set

We next compared the dynamics of the proteome to a published cell cycle transcriptome from the same cell type and obtained with a related experimental procedure for cell synchronization (4, 20). We compared our proteomics data set with the microarray data set of Whitfield et al. (4) [experiment no. 4 (Thy-Noc)] because the experimental conditions paralleled our study. The complete microarray data set was categorized into changing (1096 probes; 893 Entrez gene IDs) and nonchanging (39,484 probes; 23,133 Entrez gene IDs) genes. Our complete identified proteome was divided into regulated (2857) and nonregulated (3169) distinct protein groups. The International Protein Index (IPI) identifiers of the proteomics data were mapped to the microarray probes with common Entrez gene identifiers. Some of the IPIs mapped to more than one Entrez gene identifiers resulting in 11,826 probe entries in common (4602 Entrez gene IDs).

We compared our proteome results to a steady-state HeLa transcriptome from another microarray study (23), which allowed us to correlate mRNA signals with the summed mass spectrometric peptide signals of each protein, which is a good proxy for protein abundance (12). We analyzed the transcriptome data with the same statistical rigor (P < 0.01) that we applied to our proteomic data. We found 8161 probes (5791 Entrez gene IDs) out of 22,283 probes profiled (15,769 Entrez gene IDs) on the chip were expressed as mRNA transcripts. There was a 63% overlap between our proteome and the transcriptome on the basis of common Entrez gene IDs and both data sets showed expression of a similar number of genes (fig. S2A).

Although steady-state transcript abundance (23) did not correlate well with protein abundance (fig. S2, B and C), as we have observed in several other experimental systems (24, 25); 59% of the genes that significantly changed as cells progressed through the cell cycle at the transcriptome level (4) also changed more than fourfold at the proteome level (highest to lowest abundance). A larger fraction of the proteome exhibited a fourfold change throughout the course of the cell cycle compared to the transcriptome (21% versus 5%), which suggests that the proteome may be more dynamic than the transcriptome. Although some of the differences in the proportion of the transcriptome and proteome regulation may reflect technical limitations, it is also likely that the greater change in proteome reflects additional posttranscriptional regulation of protein abundance compared to the regulation of transcript abundance. GO analysis confirmed a cell cycle function for the proteins encoded by genes regulated at both the mRNA and the protein levels, whereas genes regulated at neither of these levels encoded proteins preferentially involved in homeostasis and basic metabolic processes (Fig. 2C). A major functional class of genes that was specifically regulated at the protein level but not at the mRNA level encoded components of the transcriptional machinery. It is intriguing to speculate that this finding is due to another level of protein expression control, for example, by microRNAs that tightly regulate the protein abundance of key cell cycle players (26). Differential protein solubility in different cell cycle stages could, in principle, lead to artifactual differences between mRNA and protein regulation; for example, chromatin proteins in mitosis when the nuclear envelope is dissolved could appear to be more abundant compared to other stages in the cell cycle. However, we measured the abundance of chromatin proteins in both the soluble and the insoluble fractions and found that chromatin proteins changed to a similar extent in both fractions.

Analysis of protein degradation motifs

Together with CDK activation, proteolytic destruction of regulatory proteins is a key regulatory mechanism in cell cycle control that ensures progression through the cell cycle in a unidirectional and nonreversible manner. We analyzed the representation of genes encoding proteins containing degradation signals in our data, comparing the presence of the degradation signals in genes regulated at three levels: transcription, protein abundance, and protein phosphorylation. Degradation signals of the KEN box and PEST region types, but not D-box degradation signals, were significantly overrepresented in proteins that are cell cycle regulated. The enrichment was observed for those regulated at the level of transcription (KEN boxes, P < 10 × 10−18; PEST regions, P < 0.002), those regulated at the level of protein abundance (KEN boxes, P < 0.002; PEST regions, P < 0.02), and those regulated by periodic phosphorylation (KEN boxes, P < 10 × 10−4; PEST regions, P < 0.001) (Materials and Methods). As expected, proteins that are regulated at multiple levels were more enriched in degradation signals than the proteins regulated by a single mechanism; for example, proteins that were regulated by phosphorylation and transcription showed a 2.6-fold enrichment for KEN boxes (P < 10 × 10−4), whereas those regulated by phosphorylation alone showed only a 1.4-fold enrichment (P < 0.003) (table S4).

The regulated cell cycle phosphoproteome

To determine phosphorylation sites that show dynamic profiles due to changes in phosphorylation state rather than due to changes in protein abundance, we normalized the measured phosphopeptide ratios by the corresponding changes in protein abundance. Phosphorylation of more than half of the phosphorylation sites changed at least twofold over the cell cycle (highest to lowest amounts); of these, again about half were maximally phosphorylated in M phase (Fig. 3A). Compared to changes in phosphorylation in cells responding to a single stimulus, such as exposure to a growth factor (11), this is a much larger proportion of the phosphoproteome, highlighting the involvement of many more regulatory processes in the cell cycle than in a single signaling pathway.

Fig. 3

Dynamics of the phosphoproteome during the cell cycle. (A) Clustering of regulated phosphorylation sites in all cell cycle stages as described for the proteome in Fig. 2A. More than half of all identified regulated phosphorylation sites peak in mitosis. (B) Reproducibility of phosphopeptide ratios between biological replicates. Phosphopeptide ratios displayed are calculated as the relative difference between nocodazole-arrested cells and asynchronously grown cells. The graph is a density plot, where red represents high density and blue represents low density. (C) Dynamic profile of two CDK1 phosphopeptides during the cell cycle. The activating site T161 peaks in mitosis, whereas phosphorylation of the inhibitory sites T14 and Y15 is decreased in mitosis. (D) Heat map of cell cycle–regulated kinase substrates. The NetPhorest algorithm was used to predict kinase-substrate relationships of all serine and threonine phosphorylated proteins in our data set for all kinases available in this algorithm. The heat map shows over- (yellow) and underrepresentation (blue) of predicted kinase substrates during different stages of the cell cycle compared to a background of phosphorylation sites that did not change with the cell cycle. ATM/ATR kinases were found to be globally active during G1-S and late S phase.

To assess the biological reproducibility of our measured phosphopeptide ratios, we performed a second experiment in adherent HeLa cells where we compared a nocodazole trap overnight (mitotic arrest) with an asynchronously grown cell population. We observe a correlation (R = 0.7) between the phosphopeptide ratios of mitotic sites between HeLa S3 and standard adherent HeLa cells (Fig. 3B). Inspection of known cell cycle–regulated phosphorylation sites showed the expected kinetics, as exemplified for the activation loop phosphorylation site (pT161) and the inhibitory sites (pT14 and pY15) of CDK1, which show opposing kinetics (Fig. 3C). We also confirmed the kinetic profiles derived from MS for three of the key cell cycle–regulated phosphorylation sites—CDK1 T14, polo-like kinase 1 (PLK1) T210, and histone 3 (H3) S10—by Western blotting with phosphorylation site–specific antibodies (fig. S9).

We found 27 gene products that exhibit periodic regulation of phosphorylation states, protein abundance, and transcription (table S5). Among these are two genes LMNB1 (encoding lamin B1) and TMPO (encoding lamina-associated polypeptide 2), which encode lamina-related proteins that influence nuclear envelope stability. TMPO contains a lamin B–binding domain that mediates its interaction with LMNB1 (27). This interaction is disrupted by phosphorylation of TMPO during mitosis (28). We found an M phase–specific, proline-directed phosphorylation site (S306) within the lamin B–binding domain of TMPO that was phosphorylated with a 50.2% occupancy during mitosis. Kinase site specificity matching predicts that this would be a substrate for either a CDK or a mitogen-activated protein kinase (MAPK). Phosphorylation of TMPO at this site could contribute to the mitotic disruption in the interaction between TMPO and LMNB1.

Kinase-substrate predictions

We predicted possible kinase-substrate relationships in our cell cycle phosphoproteome with NetPhorest (29) and constructed a heat map of the dynamic phosphoproteome. For each stage of the cell cycle, this map indicates the degree of over- or underrepresentation of substrates of different kinase groups (Fig. 3D). As expected, predicted CDK2 and CDK3 substrates were most highly phosphorylated in M phase. Similarly, overrepresentation of NIMA-related kinase (NEK) family substrates in G2 phase was observed, in agreement with an essential role of some of members of this kinase family in mitotic entry (30). In addition, substrates of the Ste20-like MST family appeared overrepresented in G2, pinpointing a possible global role of this family kinase at this cell cycle stage similar to the one observed for other Ste20-like kinases, such as SLK (31). Substrates of the DNA damage response (DDR) kinases ATM, the related kinase ATR (referred to jointly as ATM/ATR), and DNA-dependent protein kinase (DNA-PK) are significantly overrepresented in S phase, which is likely due to coupling between DNA replication and repair. Although this has been reported previously for individual phosphoproteins, our data show that it is a general phenomenon because 124 of the 1113 sites that peak in S phase match the ATM/ATR and DNA-PK kinase motifs.

To assess if the DDR kinase activation that we observe in S phase is due to the use of thymidine to synchronize the cells, we repeated the phosphoproteome analysis of cells arrested in S phase with aphidicolin instead. We decided to use aphidicolin because it mediates S-phase synchronization through a mode of action different from the mechanism of thymidine block; aphidicolin specifically inhibits DNA polymerase α. From the ratio plot (Fig. 4A) between the DDR kinase substrates induced by aphidicolin and thymidine, respectively, we observed a correlation between the two experiments (R = 0.7, Pearson correlation), showing that the DDR kinase activation in S phase is reproducible and suggesting that it is a universal replication stress response. The DDR kinase substrates were mapped to protein-protein interaction networks with STRING (32) and clusters in the graph were determined by MCODE (33) (Fig. 4B).

Fig. 4

Protein-protein interaction network of DDR kinase substrates. (A) Reproducibility of S phase–induced ATM/ATR/DNAPK substrates. The correlation of phosphopeptide ratios using two different S phase synchronizing agents is displayed. The x-axis phosphopeptide ratios are calculated as the relative intensity difference between aphidicolin-arrested cells and asynchronously grown cells, whereas phosphopeptide ratios on the y axis are from the thymidine-arrested cells. (B) The full network extracted from STRING database is shown and the color-coded nodes belong to 10 significant protein clusters as obtained from MCODE analysis. (C) The four highest-ranking and regulated network clusters from MCODE analysis [shown with yellow background in (B)]. The genes are color coded by their encoding protein’s maximal DDR phosphosite log2 ratio in S phase and the identified phosphorylation sites are shown below the gene names.

One tightly connected cluster of phosphoproteins contains many proteins that are involved in DNA repair, such as RAD50, PARP1, TP53BP1, and CLSPN (Fig. 4C), or in DNA replication, such as the MCM proteins (fig. S10). Supporting a role for these phosphorylation sites in DNA damage repair, another SILAC-based study of this process also identified most of these sites (6). To investigate further connections of cell cycle and DNA damage repair in our data set, we used the NetworKIN algorithm, which increases the accuracy in classifying kinase-substrate relationships by combining prediction of kinase-substrate motifs with contextual information (34). The DDR kinases (35) coordinate cell cycle checkpoints and DNA repair mechanisms and almost exclusively phosphorylate substrates with the linear sequence motif pS/T-Q (36). We identified more than 700 pS/T-Q sites in our data set, which could be phosphorylated by DDR kinases or other phosphoinositide 3-kinases (PI3Ks) or phosphoinositide 4-kinases (PI4Ks) or FAT domain–containing kinases (table S6). Approximately 350 of the sites that we identified were previously not in the 1705 recorded in phospho.ELM (37) and PhosphoSite (38). We focused on the potential ATM or ATR sites in the 27 proteins that were regulated at multiple levels (table S5 and Materials and Methods). Of the 21 sites obtained in this way, this analysis tied 14 substrates to ATM/ATR and DNA-PK (fig. S11), of which five are known to be part of the DDR apparatus. At least part of the activation of these DDR kinases may be due to DNA damage as a result of the cell cycle synchronization procedures.

Determination of absolute phosphorylation site occupancy

Large-scale phosphoproteomics has been successful at relative quantitation of phosphorylation sites between different cellular states. However, an inherent challenge in global phosphoproteomic screens is separating the change in the amount of the observed phosphorylation at any particular site from possible changes in protein abundance. Another challenge is the estimation of the functionally important parameter of phosphorylation site occupancy (the fraction of a given protein phosphorylated at a given site), which has so far eluded measurement in large-scale studies. Until now, it has only been feasible to analyze both phosphopeptides and unmodified peptides together in targeted studies of a single protein or protein complex. For example, Steen et al. estimated the ionization factor of phosphopeptides compared to the same nonphosphorylated peptide and used this knowledge to determine the absolute amount of phosphorylation for purified proteins and complexes (39, 40). In our quantitative proteome data set, we noticed that many proteins had one or a few peptides with SILAC ratios that were significantly different from the median of the other peptides of the same protein. In their phosphorylated form, these peptides had the opposite ratios of their unmodified forms. Thus, we could calculate absolute stoichiometry of phosphosites from any two SILAC states without knowing the ionization efficiency of the phosphopeptide. This approach makes use of the information encoded in the SILAC ratios of a phosphopeptide, the corresponding nonphosphorylated peptide, as well as the SILAC protein ratio (Fig. 5A). Briefly, the proportion between phosphorylated and nonphosphorylated peptide (termed “a” for the light SILAC state and “b” for the heavy SILAC state) is calculated from the observed SILAC ratio for the phosphopeptide (x), the SILAC ratio for the nonphosphorylated peptide (y), and the SILAC ratio of the protein (z), with the constraint that the sum of phosphorylated peptides and nonphosphorylated peptides divided by the number of protein molecules must remain constant between SILAC states (fig. S12). The calculation of the occupancies of an example peptide is illustrated in Fig. 5A. This calculation depends only on SILAC ratios for peptides covering the phosphorylation site and applies equally in case of missed cleavages induced by phosphorylation. However, our model assumes that (i) all three ratios are measured accurately; (ii) there are two states for the peptide, phosphorylated on a single site and nonphosphorylated; and (iii) the ratios change between the different states. As developed here, the model applies to singly phosphorylated peptides; multiply phosphorylated peptides require more complicated calculations.

Fig. 5

Absolute phosphorylation-site stoichiometry calculations. (A) Three SILAC ratios, those for phosphopeptide, nonphosphopeptide (the unphosphorylated version of the phosphopeptide), and protein ratio (the total amount of the phosphoprotein in both phosphorylated and nonphosphorylated forms), determine absolute phosphorylation-site stoichiometry (fig. S12). (B) Distribution of phosphosite occupancy in the different cell cycle stages. Fifty percent of all mitotic phosphorylation sites have occupancy of 75% or more. (C) Absolute phosphorylation-site stoichiometry of regulated substrate sites of known mitotic kinases and phosphorylation sites that conform to the Polo box–binding motif.

In our data set, 5869 sites had the required data, from which we calculated fractional phosphorylation site occupancy. For about half of all of these sites we determined a fractional occupancy of 70% or higher during mitosis (Fig. 5B). In contrast, most S-phase phosphosites show less than 10% fractional occupancy (Fig. 5B). In contrast to high fractional occupancies, low fractional occupancies are more complex to interpret because they could be caused by a low fractional occupancy in a cellular compartment in which the abundance of the protein is high but in which the phosphorylated form is not functional, coupled with a high fractional occupancy in compartments in which the protein is not particularly abundant but in which the phosphorylated form of the protein is functionally important. In Table 1, stoichiometry of selected phosphorylation sites during mitosis of functionally important phosphosites is given. Of the 1917 phosphorylation sites that both are up-regulated during mitosis and show more than 50% occupancy (table S7), 43% are substrates of proline-directed kinases (such as CDKs). Furthermore, these sites are enriched in PLK1 substrates (P < 0.05). During mitosis, CDK1 inactivates key anabolic proteins, such as components of the transcription complex TFIIIB, TFIID, and RNA polymerase II, and ribosomal S6 kinase 1 (41). Most of the 243 dynamic sites identified in our data set that are predicted to be CDK1 substrates show high occupancy during mitosis with a mean of 88% (Fig. 5C). Substrates of the other known mitotic kinases, such as PLK1, Aurora kinases, CHEK, and NEKs, also exhibit high phosphorylation-site stoichiometry in mitosis.

Table 1

List of selected examples of well-known, cell cycle–regulated phosphorylation sites and their phosphorylation-site stoichiometry during mitosis. Abbreviations for the amino acids are as follows: A, Ala; C, Cys; D, Asp; E, Glu; F, Phe; G, Gly; H, His; I, Ile; K, Lys; L, Leu; M, Met; N, Asn; P, Pro; Q, Gln; R, Arg; S, Ser; T, Thr; V, Val; W, Trp; and Y, Tyr.

View this table:

To obtain functional insights into the biological processes and cellular organization affected by a similar magnitude of phosphorylation-site occupancy, we analyzed our data set with the proteomic phenotyping approach (22). Briefly, the phosphoproteins were divided into five quantiles on the basis of their maximum phosphorylation-site occupancy (Materials and Methods) and analyzed for GO category (biological process and cellular compartment) enrichment by hypergeometric testing. Subsequently, the significant categories (P < 0.05) were collated across the quantiles and clustered (Fig. 6). This analysis was performed separately on the regulated proteins (sites for which phosphorylation site stoichiometry increased during M or S phase) and the complete set. In the regulated mitotic set, we observe high phosphorylation-site occupancy on proteins involved in mitosis and related cell cycle functions (Fig. 6A, left panel). Surprisingly, phosphorylation sites on proteins regulating metabolic processes show low occupancy in S-phase cells, which are metabolically active (Fig. 6A, right panel); whereas the same protein categories are enriched in high-occupancy sites in mitotic cells (Fig. 6A, left panel), where protein synthesis and related processes are comparatively quiescent. Therefore, we speculate that protein phosphorylation functions as an inhibitory switch on proteins that regulate metabolic processes. In terms of the phosphorylation-site occupancy of proteins grouped according to cellular compartment, we find that regulated phosphorylation of cytoplasmic proteins shows low occupancy, whereas nuclear proteins have high occupancy during mitosis (Fig. 6B, left panel). Conversely, cytoskeletal proteins with dynamic S-phase phosphorylation sites have phosphorylation sites that are highly occupied, whereas nuclear phosphoproteins are least occupied (Fig. 6B, right panel).

Fig. 6

Proteomic phenotyping of phosphorylation-site stoichiometry. (A) The upper panels show the phenotypic phosphoproteome comparison organized by GO biological process for mitotic (left) and S phase (right) cells. Proteins involved in metabolic processes have high-occupancy phosphorylation sites during mitosis, but low-occupancy sites during S phase (color scale: yellow, high overrepresentation; dark blue, high underrepresentation). (B) Proteomic phenotype analysis of GO CC level [similar analysis as in (A)]. (Color scale: red, high overrepresentation; blue, high underrepresentation).

Discussion

Here, we present the combined global analysis of proteome and phosphoproteome dynamics during the cell cycle representing ~6000 proteins and quantifying more than 20,000 unique phosphorylation sites. The quantitative proteomics approach we describe is generic and the data should be valuable resources for large-scale studies of in vivo phosphorylation dynamics at a systems-biology level. The data set is freely accessible through the PHOSIDA and MAPU databases (42, 43), as well as through the DAS system in ENSMBL (44). As illustrated with several examples here, we expect it to be useful for the cell cycle and cancer communities because it directly connects gene expression changes with protein regulatory information at a proteome-wide level.

In this study, we have investigated both relative and absolute changes in phosphorylation-site stoichiometry of more than 5000 phosphorylation sites in different stages of the cell cycle. Our analyses revealed that phosphorylation-site stoichiometry varies significantly between phosphosites on the same protein and that it is also highly dependent on the specific cell cycle phase. Additionally, sites phosphorylated by particular kinases showed unexpectedly high phosphorylation-site stoichiometry. For example, most of the regulated substrates of CDK1 and other mitotic kinases are almost fully phosphorylated, which is compatible with a model in which these proteins are inactivated by phosphorylation upon nuclear envelope breakdown.

Reversal of phosphorylation dynamics—high occupancy during mitosis and low occupancy during S phase—of nuclear matrix proteins suggests an important regulatory function for these sites because the cellular roles of these proteins are significantly different during cell division and interphase. Inactivation of protein activity by phosphorylation may be a widespread phenomenon in mitotic cells, which requires a high fractional occupancy. This mechanism would be biologically plausible, because activation by phosphorylation might be effective even if a small fraction of the protein pool is modified, whereas inactivation needs to affect the entire population to prevent the interphase activity of the protein. The high phosphorylation-site occupancy on mitotic proteins also has implications for phosphoproteomic analysis in general. About 5% of asynchronously growing HeLa cells are mitotic, a number similar to that of most other mammalian cell lines. Therefore, the mitotic, high-occupancy sites present in the cellular population could contribute many of the sites identified in large-scale proteomics experiments. We conclude that both relative quantitation between cellular states and quantitation of occupancy are desirable to help distinguish functional from nonfunctional phosphorylation sites (45) and that obtaining such information is now becoming feasible.

Materials and Methods

Cell culture and sample preparation

HeLa S3 cells were labeled for SILAC analysis, as previously described, with three different isotopic versions of lysine and arginine (10). Cells were synchronized in G1-S overnight with a thymidine block at a concentration of 4 mM (Sigma). Cells were then released from thymidine and subsequently collected at four different time points: 0 hours (G1-S phase), 2.5 hours (S phase early release), 5.5 hours (S phase late release), and 7.5 hours (S-G2 phase) after removal of thymidine. Two sets of cells were arrested overnight with nocodazole, after the 7.5-hour release from thymidine. The next morning, these cells were released for either 0.5 hour (M phase) or 3 hours (G1 phase). Western blotting and FACS analyses were performed to monitor the efficiency of the cell cycle arrest (Fig. 1, B and C). Culture dishes (5 cm by 15 cm) of synchronized HeLa S3 cells per condition (10 to 15 mg of protein) were harvested by lysis in modified radioimmunoprecipitation assay (RIPA) buffer containing 1% NP-40, 0.1% sodium deoxycholate, 150 mM NaCl, 1 mM EDTA, 50 mM tris (pH 7.5), phosphatase inhibitors (1 mM sodium orthovanadate, 5 mM sodium fluorate, 5 mM β-glycerophosphate), and protease inhibitors (Complete tablets, Roche Diagnostics), and were left on ice for 15 min. The cells were scraped, collected, and then vortexed for 2 min after which protein extracts were clarified by centrifugation at 17,000g (12,000 rpm Sorval SS-34) to pellet chromatin and other insoluble material. This insoluble pellet was redissolved in 8 M urea–1% N-octylglucoside supplemented with phosphatase inhibitors and benzonase. The soluble proteins in the RIPA extract were precipitated overnight at −20°C by addition of four volumes of ice-cold acetone. Following centrifugation, precipitated proteins were redissolved in 8 M urea–1% N-octylglucoside supplemented with phosphate inhibitors. The protein concentrations of all the fractions were determined with the Bradford assay and evaluated by SDS–polyacrylamide gel electrophoresis (SDS-PAGE) (fig. S13). Protein extracts derived from the different cell cycle arrest stages were then mixed 1:1:1 accordingly with an asynchronous cell population as the internal standard. Twenty percent of the protein mixtures were separated by one-dimensional (1D) SDS-PAGE, sliced in 20 gel plugs, and digested with trypsin in-gel (46). Thirty percent of the extracted peptide mixtures were used for quantitative proteome analysis by LC-MS, whereas the other 70% of the extracted peptides were subjected to titanium dioxide enrichment (47) in the presence of 2,5-dihydroxybenzoic acid (DHB) (48) and analyzed by LC-MS. The remaining 80% of protein mixtures were not fractionated by 1D SDS-PAGE but were directly reduced with dithiothreitol (DTT), alkylated with iodoacetamide, and subsequently digested with endoproteinase Lys-C and trypsin as described (11). The resulting peptide mixtures were either directly subjected to titanium oxide enrichment (5% of sample) or first fractionated (95% of the sample) by strong-cation chromatography followed by titanium dioxide enrichment.

Fluorescence-activated cell sorting analysis

Cell suspensions were fixed with 80% ethanol, permeabilized by treatment for 5 min with 0.25% Triton X-100 in phosphate-buffered saline (PBS), and incubated with 0.1% ribonuclease (RNase) and propidium iodide (10 μg/ml). Cellular DNA content was determined by flow cytometry with FACSCalibur (BD Biosciences Clontech) system and CellQuest software (Becton-Dickinson).

Western blotting

Cells were washed once with ice-cold PBS containing 1 mM phenylmethylsulfonyl fluoride, scraped off the plate, and resuspended in ice-cold Hepes lysis buffer [50 mM Hepes (pH 7.4), 150 mM NaCl, and 0.5% Triton X-100] containing 1 mM DTT, RNase A (30 μg/ml), deoxyribonuclease (30 μg/ml), protease, and phosphatase inhibitors (1 mM sodium orthovanadate, 5 mM sodium fluorate, and 5 mM β-glycerophosphate). After 15 min on ice, lysed cells were centrifuged at 13,000 rpm in a microcentrifuge for 15 min at 4°C. Protein concentrations in the cleared lysate were determined with the Dc protein assay (Bio-Rad), and equal protein amounts were loaded on SDS-PAGE gels. Separated proteins were transferred to nitrocellulose membranes (Whatman Schleicher & Schuell). For Western blot analysis, rabbit antibody against cyclin D (Sta. Cruz Biotechnology), mouse monoclonal antibody (mAb) against cyclin E (a gift from J. Bartek), goat antibody against cyclin A (Sta. Cruz Biotechnology), mouse antibody against cyclin B (BD Transduction Laboratories), mouse mAb antibody against α-tubulin (Sigma-Aldrich), rabbit antibody against Geminin (a gift from R. A. Lasckey), mouse mAb against Bub1 [an antibody against Bub1 hybridoma (clone 61-22) was produced after mice were injected with Bub1 recombinant protein spanning residues 1 to 318, following standard protocols; the antibody is of the immunoglobulin G1 subclass], rabbit antibody against Eg5 (49), mouse antibody against Aurora B (BD Transduction Laboratories), mouse antibody against Plk1 (50), mouse antibody against securin (Abcam), mouse antibody against TPX2 (Abcam), rabbit antibody against Kif20A (a gift from T. U. Mayer), mouse antibody against pT210 Plk1 (BD Transduction Laboratories), antibody against pT14 Cdk1, and rabbit antibody against pS10 histone 3 (Upstate Biotechnology) were used and detected by ECL Supersignal (Pierce Chemical) with a digital Fujifilm LAS-1000 camera attached to an Intelligent dark box II (Raytest).

Western blotting images were quantified with the ImageJ software (http://rsb.info.nih.gov/ij/). For quantification of Western signal, a set area around the band of interest was selected, the average pixel intensity was measured, and the average background signal from an area of identical size at a nonspecific region of the membrane was subtracted. The results were normalized to the signal for the asynchronous population.

Mass spectrometry

All experiments were performed on an LTQ-Orbitrap instrument connected to an online nanoflow HPLC system (Agilent 1100) via a nanoelectrospray ion source (Proxeon Biosystems). The tryptic peptide mixtures were autosampled onto a 15-cm-long 75-μm inside-diameter column packed in-house with 3-μm C18-AQUA-Pur Reprosil reversed-phase beads (Dr. Maisch) and eluted with a linear gradient from 8 to 40% MeCN in 2 hours. The separated peptides were electrosprayed directly into the LTQ-Orbitrap mass spectrometer (Thermo Fisher Scientific), which was operated in the data-dependent acquisition mode to automatically switch between one orbitrap full-scan and five ion trap tandem mass spectra. The tandem mass spectra were acquired with the multistage activation enabled for neutral loss of phosphoric acid [32.66, 48.99, and 97.97 atomic mass units (amu)] (51). All full-scan spectra were recalibrated in real time with the lock-mass option (52).

Data processing and analysis

Mass spectrometric data were analyzed with the in-house–developed software MaxQuant version 1.0.12.0 (13), which performs peak list generation, SILAC- and extracted ion chromatogram (XIC)–based quantitation, estimation of FDRs for search engine results, peptide to protein group assembly, as well as data filtering and presentation. The MS/MS spectra were searched against the human International Protein Index sequence database (IPI version 3.37) supplemented by frequently observed contaminants, concatenated with reversed versions of all sequences. Mascot (version 2.2.04) was used for the database search. Enzyme specificity was set to trypsin, allowing for cleavage N-terminal to proline and between aspartic acid and proline. Carbamidomethyl cysteine was set as fixed and oxidized methionine, N-acetylation, loss of ammonia from N-terminal glutamine, as well as phosphorylation of serine, threonine, and tyrosine as variable modifications. Spectra determined to result from medium- or heavy-labeled peptides by presearch MaxQuant analysis were searched with the additional fixed modifications Arg6 and Lys4 or Arg10 and Lys8, respectively, whereas spectra with a SILAC state not determinable a priori were searched with Arg10 and Lys8 as additional variable modifications. A maximum of three missed cleavages and three labeled amino acids (arginine and lysine) were allowed. The required FDR was set to 0.01 at the peptide and at the protein level and the minimum required peptide length to six amino acids. If the identified peptide sequence set of one protein was equal to or contained another protein’s peptide set, these two proteins were grouped together by MaxQuant and not counted as independent protein hits. Protein SILAC ratios are reported as the median of the ratios derived from SILAC triplets assigned to the protein. For phosphopeptides, the phosphorylation site(s) were assigned with a modified version of the PTM score (11) in MaxQuant. All high-confidence phosphosites (FDR <0.01) together with their cell cycle–dependent ratios were uploaded to PHOSIDA (http://www.phosida.com), which is a freely available phosphorylation site database and analysis (42).

PHOSIDA, MAPU, and Ensembl data upload

We uploaded identified phosphorylated peptides and phosphorylated sites along with further information ranging from evolutionary conservation to structural constraints and Swissprot annotations to PHOSIDA (42). Furthermore, we reassigned all determined peptides to gene transcripts that are annotated in the genome database Ensembl (http://www.ensembl.org) (44). With phosphorylation sites assigned to genes, the genome annotation section of PHOSIDA links directly to the Ensembl database, which displays our proteomic data on the basis of the DAS source management. The setup of PHOSIDA DAS layers for a given gene of interest in Ensembl is described in the “Background” section of PHOSIDA and also in Ensembl. We used the DAS/Proserver (53) technology to integrate our proteomic data into Ensembl. Moreover, the proteomic data (identified peptides) were uploaded to MAPU, the Max-Planck Unified Proteome database (http://www.mapuproteome.com). The detected peptides, which have been assigned to gene entries, are also available through the MAPU DAS source in the Ensembl database. PHOSIDA can be accessed freely at http://www.phosida.com.

Peak time index calculation for (phospho)-proteomic temporal profiles

The workflow of the analysis is shown in fig. S14. The fold ratios (r1 through r6) for each protein over the six time points (t1 = 1 through t6 = 6) were scaled between range [0, 1]. Then, for each protein, we calculated a time peak index (tpeak) by weighted mean of the expression ratio of maximal expression (ri = 1) at time point ti with respect to its adjacent time points (ti − 1 and ti + 1). To maintain the cyclicity, we made two assumptions: (i) If the maximal expression was at t1 (r1 = 1), then t1 was preceded by t0 = 0 with expression r6; and (ii) if the maximal expression was at t6 (r6 = 1), then t6 was followed by t7 = 7 with expression r1. The equations for the peak time (tpeak) calculation are as follows:

tpeak={ti1×ri1+tiri+ti+1×ri+1ri1+ri+ri+1, if max (ri at i [2,5] tiri+ti+1×ri+1+0×r6ri+ri+1+r6, if max (ri at i =1ti1×ri1+tiri×ri+7riri1+ri+r1, if max (ri at i =6 

The protein expression profiles were subsequently ordered in increasing order to get a temporal map of the cell cycle. For the purpose of rendering according to their increasing tpeak, the original (unscaled) expression profile for each protein was z-transformed before rendering as in Fig. 2A.

Cyclic angular peak calculations based on peak time index of (phospho)-proteomic temporal profiles

The time peak measure for any protein j, tpeak(j) was further converted to an angular peak measure θpeak(j) in the range [0°, 360°] by the following equation:θpeak(j)=tpeak(j)mink[1,N](tpeak(k))maxk[1,N](tpeak(k))mink[1,N](tpeak(k))×360°where N is the number of proteins in the data set.

Thus, the time peak measures tpeak(j) were converted to a polar coordinate system with the radial coordinate r = 1 and the polar angle θpeak(j). By definition, the polar coordinate system has an anticlockwise orientation with 0° ray as the polar axis (fig. S15A). But to represent and analyze our proteomics data according to the standard cell cycle stages beginning with “Mitosis,” we further wished to choose a transformed polar coordinate system having clockwise orientation with 90° ray as the polar axis (fig. S15B). In such a polar coordinate space, θpeak(j) has to be transformed by the following equation:

θpeak(j)*=90°θpeak(j)

These values were used to render the z-transformed protein profiles in the transformed polar coordinate space as shown in Fig. 2B.

Enrichment analysis for GO cellular component based on circular statistics

Here, we explain the circular statistics based testing procedure for GO cellular component (CC) categories (54). For each of the GO category C, the circular peak angles of the complete protein set (N proteins) was used to derive a 1 × N vector θC such that its jth entry was θ*peak(j) if protein j was annotated with C else “NA.” This θC vector was tested for nonhomogeneous distribution across the unit circle in the transformed polar coordinate system (θ*) with the “Rayleigh test” (55). Only categories that had a P value < 0.05 were considered significant. The mean (θC) provided mean direction of enrichment for the category C and was used to render the category at particular angles in Fig. 2B.

Comparison with Whitfield microarray data set

The microarray data set of (4) [experiment no. 4 (Thy-Noc)] was chosen for comparison with our proteomics data set because the experimental conditions therein paralleled our study. The complete microarray data set was categorized into changing (1100 probes) and nonchanging (39,484 probes). The complete identified proteome was divided into regulated (2857 IPIs) and nonregulated (3169 IPIs) proteins. The IPI identifiers of the proteomics data were mapped to the microarray probes with common Entrez gene identifiers. Some of the IPIs could be mapped to more than one Entrez gene identifiers and hence were multiplicatively mapped for each identifier, thereby resulting in 11,826 probe entries in common. This final mapped data set was categorized into four classes as shown in Table 2.

Table 2

mRNA compared to proteome over the cell cycle. Genes were divided into four different clusters; the number of proteins in each cluster is indicated. Data are from (4) and this study.

View this table:

These sets of clusters were then analyzed with GO-based clustering method for enriched biological processes (BPs) by a method similar to the one described in the later section “GO and KEGG pathways enrichment-based clustering for protein groups based on peak time.”

Comparison with steady-state HeLa microarray data

We used a published microarray data set (23), which used the Affymetrix HGU133A GeneChip with 22,283 probe sets that map to 12,999 Entrez gene identifiers. We downloaded the data from the four control data sets representing the normal HeLa transcriptome. In accordance with our proteomics experiment, we defined a transcript as present if the MAS5 P values were at least 0.01 in three of the four experiments. (The MAS5 probability values are a standard measure of the presence of a transcript. They are calculated from the signals of the different elements in each probe set.) We then mapped the 8161 probe sets with a “present call” on the Affymetrix probe set to 5791 unique Entrez gene identifiers. Our quantified proteome of 6026 IPI entries was mapped to 5455 unique Entrez gene identifiers. Subsequently, the overlap between the two data sets was calculated by common Entrez gene identifiers. Figure S2A shows their overlap with the proteomic HeLa data set.

GO and KEGG pathways enrichment-based clustering for protein groups based on peak time

The enrichment analysis for GO (54) BP and CC were done separately for each of the peak clusters (M peak, G1 peak, G1-S peak, early S peak, late S peak, G2 peak) derived from peak time index clustering (Fig. 2A) with respect to the whole quantified proteome by the “conditional hypergeometric test” available in the GOstats (56) package in the R statistical environment (57). For further hierarchical clustering based on GO terms, we first collated all the categories obtained after enrichment along with their P values, and then filtered for those categories which were at least enriched in one of the clusters with P value < 0.05. This filtered P value matrix was transformed by the function x = −log10(P value). Finally these x values were z-transformed for each GO category. These z scores were then clustered by one-way hierarchical clustering (Euclidean distance, average linkage clustering) in Genesis (58). KEGG pathway (59) enrichment analysis was done in the same way, except that the hypergeometric test was used and the reference set was complete human KEGG annotations.

Analysis of kinase-substrate relationships during phases of the cell cycle

To predict kinase-substrate relationships, all identified class I serine and threonine phosphorylation sites (pS/T) were scored with the NetPhorest algorithm (29). The overall score distribution of the cycling pS/T was significantly different from the score distribution of the noncycling counterparts (P < 10-10, χ2 test). For each time point in the cell cycle, we investigated which kinases contribute most to the observed differences in score distributions compared to noncycling pS/T. This was plotted with the “heat map” package in R (the time points in the cell cycle were z-score scaled). Over- and underrepresented kinase predictions were colored yellow and blue, respectively (Fig. 3D).

New candidates in the DDR network

We found 587 pS/T sites in the data set that matched the pS/T-Q consensus motif for kinases of the DDR, PI3K, or FAT families. Out of these sites, we defined a dynamic subset that was regulated on multiple levels. This subset consisted of 13 sites regulated by phosphorylation (high confidence) and protein abundance, as well as 8 sites regulated by phosphorylation (intermediate confidence) and mRNA abundance. We next used the NetworKIN algorithm (34) to classify which upstream kinases target the particular phosphorylation sites and found that 15 pS/T-Q sites in 14 proteins were contextually linked to ATM or DNA-PK (fig. S11).

STRING network analysis for the DDR network

The HGNC identifiers corresponding to our protein IPI identifiers were searched against the STRING database version 8 (32) for protein-protein interactions. Only interactions between the proteins belonging to the searched data set were selected, thereby excluding external candidates. STRING defines a metric called “confidence score” to define interaction confidence; we fetched all interactions that had a confidence score ≥ 0.7 (high confidence). The resulting interactome had 242 nodes and 560 interactions that we call DNA damage response network (DDRN).

MCODE-based node cluster analysis

The DDRN was analyzed for densely connected regions with a graph theoretical clustering algorithm, molecular complex detection (MCODE) (33). MCODE is part of the plug-in toolkit of the network analysis and visualization software Cytoscape (60). The four highest-ranking modules that contained up-regulated phosphorylation sites were selected for further analysis and rendering.

Cell-cycle phosphoproteome occupancy stoichiometry-based proteomic phenotyping

For the phosphopeptides shared between multiple protein identifiers, the identifier with the maximum occupancy stoichiometry was used in the analysis. The identifiers were then divided into five stoichiometry-based quantile classes corresponding to percentage cutoffs of 0, 15, 30, 70, 85, and 100%. The subsequent proteomic phenotyping was done for these quantiles for GO BP and CC as described earlier (22).

The enrichment analysis for GO BP and CC were done separately for these quantiles with respect to the whole quantified proteome by hypergeometric test available in the GOstats package (56) in the R statistical environment (57). For hierarchical clustering, we first collated all the categories obtained after enrichment along with their P values and then filtered for those categories with at least two membership counts and that were at least enriched in one of the quantiles with P value < 0.05. Categories that did not have a P value after collation in any quantile were assigned a conservative P value of 1. This filtered P value matrix was transformed by the function x = −log10(P value). The x values were z-transformed for each GO category. These z scores were then clustered by one-way hierarchical clustering (Euclidean distance, average linkage clustering) in Genesis (58).

Acknowledgments

We thank members of the department of proteomics and signal transduction at the Max Planck Institute for Biochemistry for helpful discussions and comments on the manuscript. The Protein Research Center (CPR) is supported by a generous donation from the Novo Nordisk Foundation. This work is supported by a European Union research directorate sixth framework grant (Interaction Proteome grant LSHG-CT-2003-505520). M.V. is supported by a grant from the Dutch Cancer Society. A.S. is supported by a grant from the Spanish Education and Science Ministry. T.S.J. is supported by an FP6 Network of Excellence grant (ENFIN LSHG-CT-2005-518254). This work was in part funded by the Bundesministerin für Bildung und Forschung QuantPro grant (0313831D).

Supplementary Materials

www.sciencesignaling.org/cgi/content/full/3/104/ra3/DC1

Fig. S1. GO analysis of protein and phosphoproteins subcellular localization.

Fig. S2. Comparison of steady state HeLa transcriptome with quantified proteome.

Fig. S3. Peptide and protein ratio reproducibility.

Fig. S4. Benchmark curves for periodically regulated proteins.

Fig. S5. Benchmark curves for periodically phosphorylated proteins.

Fig. S6. Proteomic phenotyping of biological processes (BP).

Fig. S7. Proteomic phenotyping of cellular components (CC).

Fig. S8. Proteomic phenotyping of KEGG pathways.

Fig. S9. WB profiles and MS data of key cell cycle regulated phosphorylation sites.

Fig. S10. MS identification of S-phase induced MCM phosphopeptides.

Fig. S11. NetworKIN analysis to identify potential substrates of the DDR kinases.

Fig. S12. SILAC based absolute phosphorylation site stoichiometry equations.

Fig. S13. 1D-SDS PAGE loading control of lysates.

Fig. S14. Data analysis workflow to cluster and circularize proteome and phosphoproteome changes.

Fig. S15. Coordinate system for circularization.

Table S1. List of all identified proteins.

Table S2. List of all identified phosphopeptides.

Table S3. List of identified protein kinases.

Table S4. Enrichment of degradation signals for cell cycle-regulated proteins.

Table S5. Protein regulated at multiple levels during the cell cycle.

Table S6. List of DNA damage response kinase substrates identified according to the kinase motif pS/pT-Q.

Table S7. List of regulated mitotic phosphorylation sites with high occupancy.

Table S8. Quantification of Western blots using ImageJ software.

Table S9. Overlap of quantified HeLa proteome with known cycle protein complexes.

References

References and Notes

View Abstract

Related Content

Navigate This Article