Insulin signaling regulates longevity through protein phosphorylation in Caenorhabditis elegans

Insulin/IGF-1 Signaling (IIS) is known to constrain longevity by inhibiting the transcription factor FOXO. How phosphorylation mediated by IIS kinases regulates lifespan beyond FOXO remains unclear. Here, we profile IIS-dependent phosphorylation changes in a large-scale quantitative phosphoproteomic analysis of wild-type and three IIS mutant Caenorhabditis elegans strains. We quantify more than 15,000 phosphosites and find that 476 of these are differentially phosphorylated in the long-lived daf-2/insulin receptor mutant. We develop a machine learning-based method to prioritize 25 potential lifespan-related phosphosites. We perform validations to show that AKT-1 pT492 inhibits DAF-16/FOXO and compensates the loss of daf-2 function, that EIF-2α pS49 potently inhibits protein synthesis and daf-2 longevity, and that reduced phosphorylation of multiple germline proteins apparently transmits reduced DAF-2 signaling to the soma. In addition, an analysis of kinases with enriched substrates detects that casein kinase 2 (CK2) subunits negatively regulate lifespan. Our study reveals detailed functional insights into longevity.

In this manuscript entitled "Lifespan regulation by insulin signaling through phosphorylation of proteins beyond FOXO" Dong and colleagues identified and characterized protein phosphorylation sites differentially regulated by insulin/IGF-1 signaling (IIS) in C. elegans. By performing a large scale quantitative phosphoproteome analysis, they first identified many differentially phosphorylated sites between long-lived daf-2/insulin/IGF-1 receptor mutant and wild-type and other short lived control animals. Following initial phosphoproteome analysis that validated their data, the authors employed iFPS, a prioritization algorithm based on multinomial logistic regression, which they developed for this study, to assess the functionality of the phosphosites. They then characterized three newly identified phosphosites in three proteins, AKT-1, EIF-2α and CDK-1. They showed that phosphorylation of AKT-1 at T492, which was upregulated in daf-2 mutants, inhibited the nuclear localization of DAF-16 and constituted a negative IIS feedback loop. They found that hyperphosphorylation of EIF-2α at S49 by GCN-2 reduced translation and contributed to longevity in daf-2 mutants. In addition, their data suggest that reduced phosphorylation of CDK-1 at T179 possibly by WEE-1.3 kinase in germline contributed to longevity in daf-2 mutants likely by impairing germline proliferation. Lastly, they proposed that reduced CK2 activity also underlies hypophosphorylation of many proteins in daf-2 mutants and contributes to longevity. Overall, their analysis of phosphoproteome data combined with subsequent functional study is important and of high quality. Although the paper has innate weaknesses of the lack of focus, because they characterized three proteins at the same time, the findings have many novel aspects and are exciting. Following are my concerns that the authors need to address.
Major comments: 1. Strictly speaking, it is difficult to state the author developed an entirely new machine-learningbased tool. They used a widely used multinomial logistic regression. The new thing is the combination of six biologically relevant variables in iFPS. Please manifest the most original feature of iFPS. In addition, they need to describe the methods and the results regarding the iFPS more in detail. 1) For example, are there weights on different criteria among the six variables? Either weight or no weight, please explain the logic why they did or did not employ the weights. Related to this issue, the authors also need to explain the biological relevance of the iFPS score better.
2) In addition, are phosphosites in figure 2D top 27? The author mentioned that the top 5% ranking as an initial cutoff for prioritizing the phosphosites regulated by daf-2. Based on their supplementary dataset, it does not seem to be the case.
3) What do the iFPS scores mean? The author did not include the iFPS scores of phosphosites in Fig.  2D. Some have negative scores (e.g. . I am confused with this part. 2. On page 8 (line 178 to 180), the authors mentioned that phosphorylation of human AKT1 T450 (equivalent to worm AKT-1 T492) is co-translationally phosphorylated by mTORC2. They need to test this in C. elegans by using rict-1 mutations. This is an important experiment because the results will be much more relevant for this phosphoproteome study than other experiments shown in figure 3. 3. Regarding Fig. 3E and Fig. S3E, with the model in Fig 3E, they need to test whether AKT-1::GFP levels are decreased in daf-16; daf-2 mutants compared to daf-2 mutants. In addition, they need to measure the levels of AKT-1 T492A::GFP in WT, daf-2 and daf-16; daf-2 double mutants. 4. They mentioned, "..among the 448 phosphoisoforms which were 266 differentially regulated in the daf-2 mutant (Fig. 2B), 124 apparently require daf-16, while 123 do not 267 (Supplementary Table  3)." Can they further analyze and discuss these two groups? For example, are canonical IIS kinase cascade proteins overrepresented in the DAF-16-independent hypophosphorylated ones? If not why is the case? Any potential crosstalk with other longevity pathways? These should be discussed in detail. 5. The authors performed lifespan assays with FUDR in main figures (Fig. 3, 4, and 6). To exclude possible confounding effects of FUDR on lifespan, they need to re-perform key lifespan assays without FUDR treatment. In addition, please specifically write the procedure of FUDR treatment including the concentrations in Methods. this study. It would be more correct to state that it promotes correct localization of AKT-1. Line 315: Please redraw figure 5D. This figure deviates from the other lifespan experiments, showing bars instead of entire lifespan curves. Why did the authors choose bars? Also the +/+ control is missing. The +/-signs are not explained in the figure legend and ambiguous (+ = not inhibited by auxin-induced degradation or + = auxin administered?). I assume the authors designated + to 'not inhibited by auxin'.

Discussion
Line 442-446: please rephrase this long sentence. It is not clear. Lines 465-473: In this section, the authors try to demonstrate the usefulness of their phosphoproteomics pipeline. However this section is not clear at all. Please rephrase this part.

Materials and methods
Lines 560-561: This is not an appropriate way to handle the data as it assumes that 1) phosphorylation will be directly related to the lifespan phenotype or 2) that single measurements are accurate enough for comparison. I understand the motivation of the authors to substitute the data, but this should be explicitly indicated in the tables and text where any of these values are used or interpreted. Line 673: Lifespans were run using the AID degron system. Were lifespan plates administered with auxin only once (with risk of degradation of auxin over time) or resupplemented at several occasions?
Typos and small errors Line39: Insulin/Insulin-like Growth Factor 1 Line 54-55: did the authors mean 'post-translational regulation'? Line 85: '…has not been surveyed rigorously. Seeking to increase…' Line 129: please explain PTM as post-translational modification at first use. Line 261: small mistake in figure 4F. At the bottom left, the text should be "Low translation rate and a long lifespan". Line 282: small mistake in Figure 5A. Worm CDK-1 has phosphorylation sites at T32 and Y33 (not 33 and 34). In the text, this is mentioned correctly. Also in the left drawing of inactive human CDK1, T161 should not be phosphorylated. Line 341: Shouldn't MNK-1 be bold and black in Fig 6A and be incorporated in the table? Line 345: …a motif X-analysis… Line 350: Explain TBB at first use in the manuscript. Line 377: …three phosphorylation sites we … Line 702: use micro symbol, not uM.
Reviewer #3 (Remarks to the Author): In this manuscript the authors have used phosphoproteomics to study the changes in protein phosphorylation in a C. elegans mutant of the insulin receptor (daf-2) as compared to the changes in phosphorylation in a FOXO mutant (daf-16), a double mutant (daf-2/daf-16) or wild type. The phosphorylation of each mutant or wild type was collected and compared by stable metabolic labelling. The authors first show that the daf-2 mutant has more coherent and different phosphorylation levels when compared to the other mutants or WT. From this the authors then focus on the changes occurring in the daf-2 mutant. As not all phosphorylation sites are likely to be equally important the authors then devised a computational method to rank phosphosites according to functional importance. From the phosphosites changing the daf-2 mutant they found 27 ranked highly by their method and from these they focused on 3 (AKT-1 T492, EIF-2α S49, and CDK-1 T179) that were characterized further in detail. The extensive experimental follow up work showed the relevance of these phosphorylation sites in C. elegans. There are many useful advances in this project including an expansion of the knowledge of protein phosphorylation in an important model system, the prediction of functionally important phosphosites, the knowledge of changes in phosphorylation in insulin signalling mutants and the detailed characterization of the phosphosites selected for analysis. This seems to be a useful for resource and study for a wide audience. I have some concerns that I will detail below, focused primarily on the large scale data and computational analysis part of the manuscript that is a better fit to my own background.

Major concerns
The authors start the work by generating a large collection of phosphorylation changes in 3 mutant backgrounds (daf-2, daf-16 and the daf-2;daf-16 double mutant). Right at the start I have difficulties in understanding exactly what was measured and reported in the Table S1. According to the schematic in Figure 1, the mutants were labelled with N14 and the WT with N15 suggesting that the quantifications are fold changes of mutant vs WT. However, then in the Table S1 the WT replicates have values and they appear to be also reported as ratios. Looking at the values in the table I suspect these are log2(ratios) but that is also not very clear. The authors then cluster the data from the WTs and mutants together as shown in Supp Figure 1D which again shows values for WT and mutant at the replicate level. The authors need to clarify: exactly how the quantifications were made; how the changes are reported in the table and what values exactly were clustered in the Supplementary Figure  1; how are the log2(ratios) calculated for each replicate.
The analysis of the data for the different types of mutants is superficial. I understand that the work was focused very much on the changes occurring in the daf-2 mutant but the clustering shown in Supp Figure 1D appears to suggest a more complex picture than simply saying that the daf-16 and daf-2;daf-16 double mutant is WT-like. Reading through the manuscript it seems like the authors could have easily just reported the WT and daf-2 data since they practically ignore the information from the other 2 mutants.
-Is there an obvious expectation in the field as to why the daf2;daf-16 double mutant should show such different patterns of phosphorylation than the daf-2 mutant alone ? If so, it would be useful to have this explained for those outside the field. -In the clustering of the data there appears to be 3 different clusters with one cluster containing all of the double mutants. The authors could try to investigate and devote more of the manuscript to the differences between these 3 clusters before delving specifically into the daf-2 mutant.
The training of the predictor for phosphosite functional importance is an important contribution in this work. There is however some concern that the number of true positives may be too small to generate a robust predictor. -The authors seem to report just one result for the area under the ROC curve that appear to be the result of training on a single negative dataset. It would be useful to have a sense of the variation in the AUC when doing multiple samples of the negative set.
-Please report the relative importance of each feature to the trained model and compare it with the importance of each feature when used as a stand-alone feature. It would be interesting to discuss briefly which features were most important for the model. -When coming up with the final score for each phosphosite, it could be useful to have each phosphosite receive an average score from several trained models via a cross-fold training procedure. This could avoid issues with lack of robustness based on the small set of positives and a single training on a random set of negatives.
-It is unfortunate that only 31 true positive phosphosites are in the dataset collected here. Are these 31 phosphosites within the full set of phosphosites collected for C. elegans in the past ? How confident are the authors that the TP list of phosphosites has been seen to be phosphorylated in C. elegans ? -I didn't fully understand how the feature "Protein-protein interactions (PPI) domain" was encoded. For a given phosphosite position found at position X of a domain type Y, did they count within 3DID what are all the different types of domains/motifs that are found to interact with this domain Y position X ? -Did I understand correctly that having a acetylation site within 15 amino-acid residues made it less likely that a phosphosite was considered to be functional? This is unexpected as one would think that having other PTMs nearby could be a positive sign that it would be more functional.
Minor comments -For the kinase substrate enrichment analysis it could have been more sensitive to run a test that compares the fold changes of the set that is linked to a kinase against the background of all phosphosite changes (z-score, KS tests for example). The outcomes are likely to be similar than the approach used but could be something to consider for future work.
-Some of the phosphosites chosen for follow up studies are fairly well characterized in other model systems. This does not detract from this study but it would be worth also mentioning this explicitly when defining the selection of sites.
-Prior studies comparing gene deletions with WT using phosphoproteomics have shown that it is difficult often interpret the changes in the context of signalling pathways. For example, a study of yeast kinase KOs (Bodenmiller et al. Sci Signalling 2010) has shown that the steady state differences, when compared with the WT, show many changes in phosphosites that are not direct substrates of these kinases. This is not unexpected as there are many indirect effects of deleting a gene and the cell adapts to the absence of that gene. This is a well known phenomena that limits the interpretation of such studies. It would be useful to add such points to the discussion and perhaps to discuss what could be done in the future. Maybe comparing stimulated vs unstimulated in the different mutants could be relevant.

REVIEWER COMMENTS
Reviewer #1: 1. In this manuscript entitled "Lifespan regulation by insulin signaling through phosphorylation of proteins beyond FOXO" Dong and colleagues identified and characterized protein phosphorylation sites differentially regulated by insulin/IGF-1 signaling (IIS) in C. elegans. By performing a large scale quantitative phosphoproteome analysis, they first identified many differentially phosphorylated sites between long-lived daf-2/insulin/IGF-1 receptor mutant and wild-type and other short lived control animals.
Following initial phosphoproteome analysis that validated their data, the authors employed iFPS, a prioritization algorithm based on multinomial logistic regression, which they developed for this study, to assess the functionality of the phosphosites. They then characterized three newly identified phosphosites in three proteins,  and CDK-1. They showed that phosphorylation of AKT-1 at T492, which was upregulated in daf-2 mutants, inhibited the nuclear localization of DAF-16 and constituted a negative IIS feedback loop. They found that hyperphosphorylation of EIF-2α at S49 by GCN-2 reduced translation and contributed to longevity in daf-2 mutants. In addition, their data suggest that reduced phosphorylation of CDK-1 at T179 possibly by WEE-1.3 kinase in germline contributed to longevity in daf-2 mutants likely by impairing germline proliferation. Lastly, they proposed that reduced CK2 activity also underlies hypophosphorylation of many proteins in daf-2 mutants and contributes to longevity. Overall, their analysis of phosphoproteome data combined with subsequent functional study is important and of high quality. Although the paper has innate weaknesses of the lack of focus, because they characterized three proteins at the same time, the findings have many novel aspects and are exciting. Following are my concerns that the authors need to address.
Ans: We are grateful that you appreciate the quality of this study and the discoveries it brings. Also, thanks for the critical but constructive comments. We have addressed all of the concerns raised by conducting experiments, re-analyzing data, and revising the text.
We believe that the revised manuscript is a much improved one and we want to thank you for your help.

Major comments:
Strictly speaking, it is difficult to state the author developed an entirely new machine-learning-based tool. They used a widely used multinomial logistic regression.
The new thing is the combination of six biologically relevant variables in iFPS. Please manifest the most original feature of iFPS. In addition, they need to describe the methods and the results regarding the iFPS more in detail.
Ans: We agree to your opinion that iFPS is not an entirely new machine learning-based tool. In iFPS, 6 sequence or structure features were encoded and integrated, and 5 of them have been used for similar purposes in previous studies. Only the number of interacting domains and/or motifs (IDMs) was a new feature introduced in this study. From an additional analysis, we found that each feature was weakly informative and not enough to achieve a superior accuracy (Fig. 2B). Combination of the 6 types of features markedly increase the performance. Based on your concern, we revised the manuscript as below: Page 9, paragraph 2, added, "…In iFPS, 5 frequently used sequence and structure features were integrated to evaluate the functionality of a candidate phosphosite, including the number of predicted upstream kinase families (UKFs) 18, 32 , the phosphorylation conservation (PhC) 18,30,32,33  We developed iFPS as a three-step method to computationally prioritize functionally important phosphosites in C. elegans, including individual feature encoding, feature integration and model training, and normalization of predicted scores. 1) Individual feature encoding. In this step, 6 sequence or structure features were encoded as below: i) UKF 18, 32 . Functional phosphosites are often regulated by multiple important kinases, and act as pivotal hubs to link various biological processes and signaling pathways. To assign potential UKFs for individual phosphosites, we used a previously developed tool named iGPS, which combined sequence profiles specifically recognized by difference kinases and protein-protein interactions (PPIs) between kinases and substrates to predict ssKSRs 50 . iGPS supported species-specific predictions for 5 model organisms, including H. sapiens, Mus musculus, Drosophila melanogaster, C. elegans and S. cerevisiae 50 . In iGPS, there were 44 and 15 specific predictors for S/T kinases and tyrosine kinases in C. elegans, respectively. To enable a higher coverage for phosphosite annotation, the low thresholds were adopted with false positive rates (FPRs) of 10% for S/T kinases and 15% for tyrosine kinases, respectively. From predicted ssKSRs, the number of UKFs was counted for each phosphosite.
ii) PhC 18, 30, 32, 33 . As previously described 69 , the residue conservation score (RCS) was calculated to measure the conservation of each phosphosite. First, the potential orthologs of worm phosphoproteins were pairwisely detected in other six eukaryotes, including H. sapiens, M. musculus, R. norvegicus, D. melanogaster, S. cerevisiae, and S. pombe. The classical method of reciprocal best hits (RBH) 81 was adopted, using the mainstream sequence alignment tool BLAST 81 . Then, protein sequences were multi-aligned by MUSCLE 82 for each worm phosphoprotein and its orthologs if available.
The RCS value was calculated for each worm phosphosite as below: For each phosphoprotein, its sequence in FASTA format was directly submitted to NetSurfP, and the RSA scores of known phosphosites were reserved.
vi) SS 18, 32 . The structural environment around phosphosites is also important for its function. Again, NetSurfP v1.1 88 was adopted to predict the secondary structures of phosphosites. The probability scores of α-Helix, β-strand and Coil were reserved for each worm phosphosite.
2) Feature integration and model training. For each worm phosphosite, the numerical values of the 6 types of sequence and structure features were separately obtained, and the initial weight value of each feature was equally assigned as 1. The MLR algorithm was implemented in Weka 3.8 80 for model training, in which the weight values were automatically determined based on the highest AUC value from the 10-fold cross-validation. Because the number of negative samples was much larger than the positive data set, we randomly picked out 121, 242, 605 or 1210 phosphosites from the negative samples, to form benchmark data sets with a positive vs. negative ratio of 1:1, 1:2, 1:5 or 1:10. By testing, the benchmark data set with a positive vs. negative ratio of 1:5 exhibited a higher AUC value ( Supplementary Fig. 8F). To avoid overfitting, we generated 10 different sets of benchmark data sets for model training ( its phosphorylation at S10 70 . However, acetyl group is hydrophobic and might interfere the phosphorylation of nearby residues. For example, acetylation of human Tau at K259, K321 and K353 markedly inhibits its phosphorylation at S262, S324 and S356, respectively 71 . In this work, the negative weight value of ASC indicated that functional phosphosites did not prefer to co-occur with nearby acetylation sites. For RSA, a previous analysis suggested that phosphosites located in disordered regions tended to be non-functional, because these phosphosites had higher accessibility to kinases and their phosphorylation might be more permissive but not rigorously regulated 72 . The negative weight value of RSA supported this hypothesis that phosphosites with lower accessibility had higher biological impacts. The results on SSs indicated that functional phosphosites were not preferentially located at β-strands." 3. In addition, are phosphosites in figure 2D top 27? The author mentioned that the top 5% ranking as an initial cutoff for prioritizing the phosphosites regulated by daf-2. Based on their supplementary dataset, it does not seem to be the case. Ans: Thank you for your sharp eyes. iFPS only considered static sequence and structure features of functional phosphosites, whereas the dynamic information of phosphoproteomic change in the daf-2 mutant was also included to filter potentially false positive predictions. To clarify this ambiguous point, we re-performed the statistical analysis, and the number of the daf-2 regulated phosphoisoforms were changed, as well as the number of predicted lifespan-related phosphosites (LiRPs). We revised the manuscript as below: Page 10, paragraph 2, revised, "Next, iFPS was applied to score all identified phosphosites, which covered 31 known Ans: Simply, a higher iFPS score denoted a higher probability of a predicted phosphosite to be truly functional. The raw scores were directly calculated by iFPS, ranged from -5.5649 to 24.9553 with any further normalization. In the revision, we normalized the iFPS scores into a range of 0 to 1, in a more readable manner. We added a step entitled "Normalization of predicted scores" in the section "The iFPS algorithm", and revised the manuscript as below: Page 40, paragraph 1, added, "3) Normalization of predicted scores. The raw scores directly predicted by iFPS ranged from -5.5649 to 24.9553 (Supplementary Data 2). Here, we normalized the scores into a range of 0 to 1. The IQ range was calculated, while upper fence and lower fence were defined as below: Where Q1 was the lower 25% quantile, and Q3 was the upper 25% quantile. To eliminate the influence of extremely higher or lower iFPS scores, the 3*IQ was adopted in this study. Phosphosites with iFPS scores higher or lower than upper or lower fence were scored to 1 and 0, respectively. The highest and lowest iFPS scores within the upper and lower fences were denoted as Smax and Smin, whereas other iFPS scores were normalized as below: Page 12, paragraph 2, added, "Utilizing the same targeted quantitation assay, we found that the C. elegans TOR complex 2 (CeTORC2) is involved in phosphorylating AKT-1 T492. RICT-1, the only homolog of human RICTOR defines CeTORC2 38 ( Supplementary Fig. 4A). We found that a loss-of-function mutation of rict-1 reduced AKT-1 T492 phosphorylation by 40% without affecting the AKT-1 protein level ( Supplementary Fig. 4B). In line with this result, unphosphorylated AKT-1 T492 peptide, which was undetectable in WT worms, became detectable in the rict-1(lf) mutant ( Supplementary Fig. 4C). We thus conclude that AKT-1 T492 is a substrate phosphorylation site of CeTORC2 and that AKT-1 may exist stably in the absence of this constitutive phosphorylation on T492." As the reviewer noted, worm AKT-1 T492 is equivalent to human AKT1 T450, for which we cited the literature finding that it is reportedly co-translationally phosphorylated by mTORC2. Our study, though, concerns IIS only. Whether or not worm AKT-1 T492 is phosphorylated by worm TORC2 changes not the conclusion of this manuscript, i.e., T492 phosphorylation is part of a feedback loop of IIS. Whether the loss of function of daf-2 alters mTOR activity, and if so in what direction, is up to debate and is a research topic on its own.
In addition, they need to measure the levels of AKT-1 T492A::GFP in WT, daf-2 and daf-16; daf-2 double mutants.  Table 3)." Can they further analyze and discuss these two groups?

Ans
For example, are canonical IIS kinase cascade proteins overrepresented in the DAF-16-independent hypophosphorylated ones? If not why is the case? Any potential crosstalk with other longevity pathways? These should be discussed in detail.
Ans: Thank you! This is a very good point. We conducted a KEGG-based enrichment analysis for the DAF-16-dependent and DAF-16-independent phosphoisoforms, and added a new Supplementary Fig. 8A and Supplementary Data 3 to present the results.
We did not see enrichment of the canonical IIS proteins in either the DAF-16-dependent or DAF-16-independent group. The phosphoisoforms that were quantified at least three times in both the daf-16 and daf-16; daf-2 worms were subject to differential analysis.
The phosphoisoforms of the IIS cascade proteins thus quantified were too few to yield meaningful results in the enrichment analysis. We revised the manuscript as below: Page 22 (Fig. 3, 4, and 6). In brief, we found that under FUDR-free conditions:
In lifespan assays performed with FUDR, we added 50 ng/μl FUDR to NGM agar before pouring plates and initiated the FUDR treatment from adult day one (see the revised Methods). 50 ng/μl FUDR completely eliminated the bagging phenotype, reduced bacterial or fugal contamination, and showed no significant effect on the WT lifespan.
We have revised the manuscript as below: Page 13, paragraph 1, revised, "…Consistently, the T492A mutation moderately but significantly extended the lifespan of WT worms by 8-17% hours. Synchronized worms were cultured on normal NGM plates until lifespan assays initiated and adult day one worms were transferred to FUdR-containing plates."

Minor comments:
In Supplementary Fig. 1E, the KEGG enrichment revealed that phosphorylation on proteins involved in FOXO signaling and longevity regulation were down-regulated in the daf-2 mutant. In Fig. 1D Ans: We have explained in details as below: Page 11, paragraph 2, revised, "To experimentally validate iFPS predictions and to flesh out the mechanism of lifespan extension by protein phosphorylation in response to reduced insulin signaling, we focused on phosphosites within the 3 prominent protein function groups for in-depth functional analysis. Among the FoxO signaling group, we were interested in AKT-1 pT492 because of its unexpected hyper-phosphorylation upon reduction of daf-2 activity.
In the rest 2 groups, we chose to validate phosphorylation changes on EIF-2α -a key component of translation initiation machinery, and CDK-1 -a master regulator of cell  Supplementary Fig. 1F). Additionally, we newly identified a phosphorylation hotspot in the HSF-1 C-terminus and two phosphosites on FKH-9 ( Supplementary Fig. 1F for regulating these phosphosites might be also important for longevity and act as potential LiRKs. Based on the site-specific kinase-substrate relations (ssKSRs) predicted with the iGPS algorithm 50 , we implemented a statistical method named pLiRK to assess whether the predicted ssKSRs of a given kinase were enriched in hyper-or hypo-phosphorylated sites of the daf-2 mutant against all identified phosphosites (One-sided hypergeometric test, p < 0.05). In total, pLiRK predicted 27 potential LiRKs with enriched ssKSRs on the hypo-phosphorylated sites, whereas no LiRKs were predicted on the hyper-phosphorylated sites (Fig. 6A). 10 of the 27 potential LiRKs have been reported to regulate lifespan. There were 8 kinases reported to extend lifespan with the treatment of RNAi or loss-of-function mutation, including the worm mTOR kinase LET-363, the MAPK activated kinase MAK-2, the cell cycle kinases CDK-1, CDK-2, CHK-1, PAR-1, PAK-1, and PDHK-2 (Fig. 6A).
Among the predicted LiRKs, the CK2 kinase was composed of the catalytic subunit KIN-3 and the regulatory subunit KIN-10. Using a motif discovery tool pLogo 51 , we found that two CK2-specific phosphorylation motifs were significantly overrepresented in the hypo-phosphorylated sites found in the daf-2 mutant (Fig. 6B), hinting a role of CK2 in daf-2 longevity. It was reported that CK2 accelerates both chronological and replicative ageing in Saccharomyces cerevisiae 52, 53 . However, a previous study in C.
Knockdown of kin-3 or kin-10 during adulthood moderately but significantly extended worm lifespan in independent trials ( Fig. 6C

1) Proteins should be written in Roman capital and genes should be written in lowercase
Italic. Please doublecheck everything regarding this issue. e.g. Fig. 6A needs to be properly changed to proteins.
Ans: Sorry about it. We went through the manuscript several more times and corrected all the errors found.
Ans: The text has been deleted and the Supplementary Fig. 1 has been updated.

In the volcano plot of
Ans: It has been added.  Fig. 1A Fig. 1C)." is better than using Supplementary Fig. 1A-C at the end.

Please state figures properly in the
Ans: It has been assigned accordingly.

On page 4 (line 73), the authors mentioned that their identification of 15,443
phosphosites was a doubling of the current collection of C. elegans phosphorylation database. However, the number of phosphosites registered in dbPAF was 10,767 in Fig.   1B. Please specify which is doubled by their identification.
Ans: 9,949 phosphosites identified in this study were not covered by dbPAF. This number was nearly equal to 10,767 known phosphosites maintained in dbPAF. Fig. 2A-H to match the sentence.

On page 6 (line 139), change the orders of Supplementary
Ans: The orders of Supplementary Fig. 2A-H are matched with the sentence. We have modified the sentence to avoid misunderstanding.

On page 7 (line 144), predicted secondary structures included alpha helix, beta
strand, and coiled-coil, but they were not explicit.
Ans: We have changed this sentence.

On page 7 (line 157), WT plus daf-16 mutants are confusing. Please rewrite this.
Ans: We have rephrased it into "WT control". We revised the manuscript to provide details on preparation of the WT control data as below: Page 34, paragraph 3, changed, "To determine the daf-2 regulated phosphorylation, phosphoisoforms quantified at least 3 times in daf-2 samples were extracted to compare with WT control. 14 N/ 15 N ratios of WT was adopted as control-1 if the number of quantitation ratios in WT samples were more than 2. Alternatively, 14 N/ 15 N ratios of WT, daf-16 and daf-16; daf-2 samples was adopted as control-2 if the phosphoisoform was quantified only once or twice in WT samples. The control-1 and control-2 were merged into a single WT control data, and Log 2 (median of daf-2/ median of control) distribution was plotted to estimate the median and normalized interquartile (NIQ) ranges…"

On page 15 (line 342), explicitly mention that KIN-3 and KIN-10 are CK2
components in the first sentence.
Ans: We have changed this sentence.

Change the bar graphs in
Ans: It has been changed. Supplementary Fig. 1D, replicates 1 and 4 of daf-2 samples were different from replicate 2 and 3. The authors should state how to handle the difference between

In
replicates. In addition, please use another description rather than "WT or WT-like lifespan". daf-16(-) and daf-16(-); daf-2(-) mutants tend to live shorter than the wild-type animals. It is inappropriate to call these animals "WT-like". Use ticks in the color bar.
Ans: The difference reflects biological variations, the source of which isn't clear at the moment. The original Supplementary Fig. 1D has been revised and moved to a new Supplementary Fig. 2D. Supplementary Fig. 2I Fig. 2E). Similarly, the enrichment of ribosome, glycerolipid metabolism, and glycerophospholipid metabolism were not evident in the daf-16 or daf-16; daf-2 mutants, suggesting that phosphorylation in these pathways are regulated by daf-2 in a DAF-16-dependent manner ( Supplementary Fig. 2E).

In
By analyzing the hypo-or hyper-phosphoproteins separately, we found that proteins involved in RNA transport had low phosphorylation levels in both the daf-2 and daf-16; daf-2 mutant, indicating a group of DAF-16-independent phosphorylation that required daf-2 ( Supplementary Fig. 2F). Notably, hypo-phosphorylated sites may be either directly or indirectly targeted by IIS, whereas hyper-phosphorylated sites are surely indirectly related to IIS. Glycerolipid metabolism and glycerophospholipid metabolism were enriched for proteins with up-regulated phosphorylation in the daf-2 mutant ( Supplementary Fig. 2F). Up-regulation of lipid metabolisms was a major phenotype of daf-2 mutants 4 . Phosphorylation changes, along with changes of gene expression 24, 29 and protein abundance 8 , might attribute to the phenotypic alteration of lipid metabolisms due to daf-2 deficiency. Taken together, these results supported the high quality of phosphoproteomic quantification, and could be highly informative for further characterization of biological processes regulated by IIS, as well as the extension of the mechanistic understanding of this field to emphasize the importance of phosphorylation-mediated regulation." The original Line 123 describes DAF-16-dependent changes. There are also DAF-16-independent changes. Lines 262-270 describe both.

Line 156-159: please rephrase this long sentence. It is not clear.
Ans: Thank you for pointing it out. We have revised this sentence. Ans: You are right. The +/-signs could be interpreted differently by people. We had intended to use the +/-signs to indicate the presence/absence of target proteins, respectively. We have changed this and the entire lifespan curves have been shown in the revised Fig. 5D.

Discussion
Line 442-446: please rephrase this long sentence. It is not clear.
Ans: It has been rewritten as follows "…As an example, we showed that hypo-phosphorylated CDK-1 T179, which inactivates CDK-1 and thereby delays the cell cycle, potentially contributes to daf-2 longevity. It indicates that hypo-phosphorylated Ans: We have modified the section as below: Page 28, paragraph 2, revised, "Notably, we also tested 2 phosphosites -EIF-5 pT376 and pS380 -that were among top 10-13% ( Supplementary Fig. 6G Ans: We apologize that the original description on preparation of WT control data was not clearly presented. The section entitled "Phosphorylation changes in IIS mutants" were carefully revised as below: Page 34, paragraph 3, changed, "To determine the daf-2 regulated phosphorylation, phosphoisoforms quantified at least 3 times in daf-2 samples were extracted to compare with WT control. 14 N/ 15 N ratios of WT was adopted as control-1 if the number of quantitation ratios in WT samples were more than 2. Alternatively, 14 N/ 15 N ratios of WT, daf-16 and daf-16; daf-2 samples was adopted as control-2 if the phosphoisoform was quantified only once or twice in WT samples. The control-1 and control-2 were merged into a single WT control data, and Log 2 (median of daf-2/ median of control) distribution was plotted to estimate the median and normalized interquartile (NIQ) ranges. Here, Q 1 was the lower 25% quantile, and Q 3 was the upper 25% quantile. Then interquartile (IQ) and NIQ ranges were calculated as below: The 14 N/ 15 N ratios of each phosphoisoforms were subjected to Wilcoxon rank-sum test.
Then a flexible filter was applied to define the daf-2 regulated phosphoisoforms, which met the criteria of either [log 2 (daf-2/control) out the range of median ± 1.5 * NIQ with one-tailed p < 0.05, Wilcoxon rank-sum test] or [log 2 (daf-2/control) out the range of median ± NIQ with two-tailed p < 0.05, Wilcoxon rank-sum test].
Similarly, phosphoisoforms quantified at least 3 times in both daf-2 and daf-16; daf-2, were subjected to statistical comparison. The DAF-16-dependent phosphorylation were defined using the same filter as the daf-2 regulated phosphoisoforms." 11. Line 673: Lifespans were run using the AID degron system. Were lifespan plates administered with auxin only once (with risk of degradation of auxin over time) or resupplemented at several occasions?
Ans: We seeded OP50 on auxin plates 12 hours before use and transferred living worms to freshly prepared plates every four days. This information has been included in the revised methods.
Page 44, paragraph 2, revised, "Auxin treatment was performed as previous description 47 . Briefly, auxin, which was dissolved in ethanol was added to NGM agar before pouring plates. The final concentration of auxin and ethanol per plate was 1 mM and 0.25%, respectively. 0.25% ethanol was used as control. Freshly prepared auxin plates were maintained at 4°C in the dark for up to two weeks. Auxin plates were seeded with OP50 12 hours before used. Living worms were transferred to fresh plates every four days. Auxin treatment was imitated from adult day one."

Typos and small errors
Line39: Insulin/Insulin-like Growth Factor 1 Ans: We have added "-like". 17. Line 282: small mistake in Figure 5A. Worm CDK-1 has phosphorylation sites at T32 and Y33 (not 33 and 34). In the text, this is mentioned correctly. Also in the left drawing of inactive human CDK1, T161 should not be phosphorylated.
Ans: We have corrected the text accordingly. Human CDK1 with phosphorylation on T14, Y15, and T161 is inactive (PMID: 21900495, 32508972). We have ordered the phosphorylation events during CDK1 activation in the revised Fig. 5A. Fig 6A and be incorporated in the table?

Line 341: Shouldn't MNK-1 be bold and black in
Ans: The lifespan phenotype of mnk-1 was controversial and not recorded in WormBase.

Line 345: …a motif X-analysis…
Ans: We have re-performed the motif analysis with pLogo because the Motif-X webpage is no longer accessible after 2019.

Line 350: Explain TBB at first use in the manuscript.
Ans: We have added the information.

Line 377: …three phosphorylation sites we …
Ans: We have rewritten this sentence.

Line 702: use micro symbol, not uM.
Ans: We have changed this symbol. Thank you for the positive comments and the helpful suggestions. We tried very hard in this study. We are glad to know that you think highly of this work.

Major concerns
The authors start the work by generating a large collection of phosphorylation changes in 3 mutant backgrounds (daf-2, daf-16 and the daf-2;daf-16 double mutant). Right at the start I have difficulties in understanding exactly what was measured and reported in the Table S1. According to the schematic in Figure 1, the mutants were labelled with N14  Figure 1; how are the log2(ratios) calculated for each replicate.
Ans: In this study, we sampled four 14 N-labeled strains: WT, daf-2, daf-16 and daf-16; daf-2 mutants. Each of the 14 N worm samples was added with an equal volume of 15 N-labeled reference worms. The 15 N-labeled reference worms were prepared by culturing WT worms for generations on 15 N-labeled food and harvesting in mixed stages.
They are intrinsically different from 14 N-labeled WT. Ans: We agree to your opinion, and re-performed all data analyses on the phosphoproteomic data. The section entitled "Phosphorylation changes resulted from genetic disruption of IIS" was largely re-written. We revised the manuscript at below: Page 7, paragraph 3, revised, "More than 15,000 phosphopeptides were quantified against their 15 N-labeled cognate peptides, which were introduced as an internal reference standard by feeding C. elegans However, we found more of the DAF-16-independent phosphorylation changes than the DAF-16-dependent ones in the daf-2 mutant.
We have discussed with this issue in the revised manuscript: Page 16, paragraph 2, "…Pursuing this, it was surprising when we found that among the 408 phosphoisoforms which were differentially regulated in the daf-2 mutant, 100 apparently required daf-16, 6. The training of the predictor for phosphosite functional importance is an important contribution in this work. There is however some concern that the number of true positives may be too small to generate a robust predictor.
Ans: We agree to your opinion. However, the phosphorylation in C. elegans was less studied in contrast to mammalians, and much fewer functional phosphosites were reported in C. elegans. We added a section entitled "Preparation of benchmark data sets for iFPS" to provide more details on this point as below: Ans: Yes. We clarified this point by revising the manuscript as below: Page 28, paragraph 1, added, "…For ASC, phosphorylation might be promoted by adjacent acetylation sites. For example, acetylation of yeast histone H3 at K9 by GCN5 facilitates its phosphorylation at S10 70 . However, acetyl group is hydrophobic and might interfere the phosphorylation of nearby residues. For example, acetylation of human Tau at K259, K321 and K353 markedly inhibits its phosphorylation at S262, S324 and S356, respectively 71 . In this work, the negative weight value of ASC indicated that functional phosphosites did not prefer to co-occur with nearby acetylation sites. …"

Minor comments
-For the kinase substrate enrichment analysis it could have been more sensitive to run a test that compares the fold changes of the set that is linked to a kinase against the background of all phosphosite changes (z-score, KS tests for example). The outcomes are likely to be similar than the approach used but could be something to consider for future work.
Ans: We agree to your opinion, and actually this method was named KSEA andpublished in 2017 (Bioinformatics, 2017, 33, 3489-3491). However, this study was started in 2015, and at that time we didn't have such a method to prioritize potential lifespan-related kinases (LiRKs). So we had to develop an alternative one named pLiRK.
Although this method is quite simple, 10 of 27 predicted LiRKs were previously reported to be involved in longevity, indicating the accuracy of pLiRK. Using KSEA, we also predicted 50 potential LiRKs, and only 9 were supported by experimental evidence (Supplementary Fig. 8C and Source data). We carefully describe the pLiRK method as below: Page 46, paragraph 3, revised, "The ssKSRs predicted by iGPS 50 were adopted for computational detection of potential LiRKs. The one-sided hypergeometric test was used to determine whether targets of any kinases were statistically enriched in the daf-2 hyper-or hypo-phosphorylated data set.
For each kinase k i , we defined the following: The enrichment ratio (E-ratio) of k i was computed, and the p value was calculated based on the hypergeometric distribution as below: Potential LiRKs were identified with significantly enriched hyper-or hypo-phosphorylated sites (p < 0.05 and E-ratio > 1)." We also added a brief discussion about the KSEA result: Page 24, paragraph 2, added, "…We also adopted the Kinase-Substrate Enrichment Analysis (KSEA) 57 to estimate kinase activities in the daf-2 mutant ( Supplementary Fig. 8C). 50 kinases were significantly (p < 0.05) enriched and all showed a predicted downregulation of kinase activity in daf-2 worms. However, only 9 of the 50 KSEA-enriched kinases are known lifespan regulators, a proportion lower than that from the pLiRK enrichment (10/27).
KIN-3 was also enriched by KSEA but not among the top 27. Our pLiRK method is thus efficient to infer kinases that are involved in daf-2 longevity.…" 14. -Some of the phosphosites chosen for follow up studies are fairly well characterized in other model systems. This does not detract from this study but it would be worth also mentioning this explicitly when defining the selection of sites.
Ans: We have added the information in the revised manuscript as below: Page 11, paragraph 2, revised: "…To experimentally validate iFPS predictions and to flesh out the mechanism of lifespan extension by protein phosphorylation in response to reduced insulin signaling, we focused on phosphosites within the 3 prominent protein function groups for in-depth functional analysis. Among the FoxO signaling group, we were interested in AKT-1 pT492 because of its unexpected hyper-phosphorylation upon reduction of daf-2 activity.
In the rest 2 groups, we chose to validate phosphorylation changes on EIF-2α -a key component of translation initiation machinery, and CDK-1 -a master regulator of cell Ans: This is an excellent discussion point! For this we have added the following paragraph in the revised manuscript.
Page 21, paragraph 2, added, "In a genetic mutant, the organism makes many adjustments to cope with the mutation until it reaches a steady state. As such, there are often many differences between the mutant and WT, some of which are direct and some, likely most, are not. A comprehensive phosphoproteomic study of the budding yeast nonessential kinases and phosphatases found that 32% and 53% of the 8,814 regulated phosphorylation events resulted from direct actions of the kinases or phosphatases examined, respectively 55 . In this study, the three functionally validated phosphosites AKT-1 pS492, EIF-2α S49, and CDK-1 pT179 are all indirectly regulated by the tyrosine kinase DAF-2. To differentiate direct targets from indirect ones, we envision that a kinase could be degraded rapidly in vivo using the AID 47 or proteolysis-targeting chimera (PROTAC) 56 method, followed by a time course phosphoproteomic analysis. For cultured cells or unicellular organisms, a kinase inhibitor could be used in place of AID or PROTAC. Since direct targets should change earlier than indirect targets, the early responding phosphosites are more likely to be direct targets of a kinase."