Temporal phosphoproteomics reveals circuitry of phased propagation in insulin signaling

Abstract

Insulin is a pleiotropic hormone that elicits its metabolic and mitogenic actions through numerous rapid and reversible protein phosphorylations. The temporal regulation of insulin’s intracellular signaling cascade is highly complex and insufficiently understood. We conduct a time-resolved analysis of the global insulin-regulated phosphoproteome of differentiated human primary myotubes derived from satellite cells of healthy donors using high-resolution mass spectrometry. Identification and tracking of ~13,000 phosphopeptides over time reveal a highly complex and coordinated network of transient phosphorylation and dephosphorylation events that can be allocated to time-phased regulation of distinct and non-overlapping subcellular pathways. Advanced network analysis combining protein-protein-interaction (PPI) resources and investigation of donor variability in relative phosphosite occupancy over time identifies novel putative candidates in non-canonical insulin signaling and key regulatory nodes that are likely essential for signal propagation. Lastly, we find that insulin-regulated phosphorylation of the pre-catalytic spliceosome complex is associated with acute alternative splicing events in the transcriptome of human skeletal muscle. Our findings highlight the temporal relevance of protein phosphorylations and suggest that synchronized contributions of multiple signaling pathways form part of the circuitry for propagating information to insulin effector sites.

Introduction

Insulin exerts its action through activation of the insulin receptor that triggers a complex intracellular signaling cascade involving protein phosphorylation of multiple target proteins¹. Insulin has pleiotropic effects such as modulating energy metabolism, cell differentiation, survival, and growth, resulting in complex, tissue-specific responses to the hormone stimulation². Compromised insulin action in skeletal muscle precedes the onset and development of type 2 diabetes³, but the cellular signaling pathways regulated by the hormone, and their relevance for the etiology of the disease are not well understood⁴.

Insulin receptor-mediated tyrosine phosphorylation of insulin receptor substrate (IRS) proteins and other cellular adapter molecules leads to the activation of several Ser/Thr kinases, including AKT/PKB, mTORC1, and S6K, which are part of a wide protein phosphorylation network regulating metabolism and various other vital cell functions⁵. Conversely, phosphorylation of specific Ser/Thr residues in IRS and other downstream targets is associated with the inhibition of insulin signaling, which in turn may contribute to counterregulatory feedback mechanisms and/or play roles in the pathophysiology of insulin resistance². Adding to the complexity, insulin also affects the activity of protein phosphatases such as PTB1B and SHP2^2,6, which dynamically alter the phosphorylation status of cellular proteins. Furthermore, multiple phosphorylation cascades overlap, possibly contributing to redundancy in insulin signaling pathways and amplification of cellular signals⁵. Collectively, the orchestrated action of insulin on the many phosphorylation sites may be essential for distributing and amplifying the insulin response within the cell and likely contributes to individual differences in insulin sensitivity and susceptibility to diabetes.

Previous studies have identified numerous insulin-regulated phosphosites using rodent cell culture models of muscle cells^7,8, induced pluripotent stem cells (IPS) from human donors differentiated into myoblasts⁹ and human skeletal muscle biopsies¹⁰ by unbiased mass spectrometry (MS)-based phosphoproteomics. Although these studies provide a deeper understanding of the metabolic and mitogenic pathways elicited by insulin, several limitations remain: (i) most studies have been informed by single time point analysis, thus not considering temporal propagation of insulin signaling; (ii) to the best of our knowledge, no phosphoproteomics study has investigated time-dependent insulin signaling in human primary myotubes that retain innate features of the donor but lack the complexity of the multicellular muscle tissue, and, (iii) analyzing the individual variation in insulin signaling of differentiated human skeletal muscle cells derived from different donors has received rather little attention in previous signaling studies.

Thus, the aim of this study is to explore the temporal development of insulin signaling and its donor variation using differentiated primary human myotubes derived from satellite cells of healthy donors. Our results identify signaling nodes conferring donor variability in insulin signal transduction and reveal the dynamics of temporal activation patterns for multiple pathways that are understudied in the context of insulin signaling such as mRNA splicing.

Results

Quantitative analysis of the insulin-stimulated phosphoproteome of differentiated human skeletal muscle cells

For a global analysis of insulin signal transduction in human skeletal muscle cells, we isolated satellite cells from muscle biopsies of five normal-weight male subjects as described in “Methods” section and Supplementary Table 1. Myoblasts were differentiated into myotubes for 6 days by serum deprivation and characterized for insulin-stimulated glycogen synthesis as previously described¹¹. After differentiation, the myotubes were serum starved for 6 h before being stimulated with 100 nM insulin for different time intervals (1 min, 2.5 min, 5 min, 15 min, 30 min, and 60 min). Moreover, we also incubated serum-starved myotubes in the absence of insulin for 2.5 min, 15 min, and 60 min as controls. Following incubation, the cells were lysed and snap-frozen, before being further processed for phosphoproteomic analysis. In addition, cell lysates were subjected to quantitative proteomics analysis as described in “Methods” section to account for potential donor-specific differences in total protein abundance.

Variation in the abundance of phosphopeptides between technical replicates for all five donors and seven time points was low with a median Pearson’s r of 0.98 for all 70 measurements, i.e. 35 pairs of technical replicates (Supplementary Figs. 1 and 2). Overall and across all time points, we quantified 20,470 peptides, including 13,196 phosphopeptides corresponding to 11,572 class I phosphosites (i.e. ≥75% localization probability), and mapping to 4415 phosphoproteins with the vast majority (>99%) being phosphorylated at Ser/Thr residues (Fig. 1a, Supplementary Fig. 3a). We also measured the total proteomes from each donor as described in “Methods” section. Briefly, myotubes at baseline, as well as 60 min after insulin stimulation, were lysed and analyzed, allowing quantitative determination of overall 25,593 peptides, corresponding to 3390 proteins. Stimulation with insulin for 60 min hardly had any effect on the protein abundance in total lysates of donor cells with only ten significantly differential proteins (adjusted p-value < 0.05; ⎸FC⎹ > 1.5; Supplementary Table 2; Supplementary Fig. 4a). We then identified 1367 peptide sequences that were quantitated in both analyses, total proteome derived from cell lysates and phosphoproteome derived from the TiO₂/Fe²⁺IMAC enrichment analysis. Likewise, 5322 phosphopeptides could be related to matching protein abundances of the total proteome. Normalization of individual phosphorylation site abundance to the corresponding peptide/protein abundance revealed that changes in phosphopeptide abundance are due to phosphorylation events and not alterations in protein abundance (Supplementary Fig. 4). As expected, the abundance of the non-phosphorylated peptidoforms correlated mainly with donor identity whereas differentially phosphorylated peptides correlated with insulin treatment (Supplementary Fig. 4e–g).

Fig. 1: Time-resolved phosphoproteomics analysis of insulin action in human skeletal myotubes.

a Study design and summary of quantified phosphopeptides. Total proteomes (25,593 quantified peptides) and phosphoproteomes (20,470 quantified peptides) were determined from tryptic peptides after RP₁₈ and TiO₂/IMAC chromatography, respectively. Myotubes from each donor were analyzed for insulin-stimulated glycogen synthesis, as detailed in “Methods” section. b Number of regulated and unregulated phosphopeptides (p-value < 0.05; ⎸FC⎹ > 1.5 vs. t₀) over time after insulin stimulation. c Donor and time-dependent coefficient of variation (CV) of peptide abundances (n = 11,612 CV values for each time point). In the box plots, median (horizontal line), minimum (lower whisker), maximum (upper whisker), 25th and 75th percentile (lower and upper hinge) values of the respective time-specific distributions of CVs are visualized. Median value 26.9% across all donors, time points. d Number of regulated phosphopeptides containing one (mono-), two (di-), three (tri-), and more phosphorylated sites across the time course of insulin stimulation. e–g Time point-specific volcano plots showing increasing and decreasing abundance of phosphopeptides after insulin stimulation relative to the basal state. Selected phosphosites are highlighted. Collectively, 822/749, 829/550, 914/548, 1023/604, 1057/543, and 1048/720 phosphopeptides were up/downregulated after 1 min, 2.5 min, 5 min, 15 min, 30 min, and 60 min (t-test p-values, not adjusted). Source data are provided as a Source Data file.

We observed only a few differences in the abundance of phosphopeptides in the absence of insulin within 60 min of incubation (i.e. p < 0.05; ⎸FC⎹ > 1.5 vs. t₀; Supplementary Fig. 5). Therefore, non-stimulated myotubes (t₀) from each individual donor served as respective basal control. Further analysis identified a varying number of phosphopeptides that were significantly differentially phosphorylated (i.e. p < 0.05; ⎸FC⎹ > 1.5 vs. t₀) over time in response to insulin stimulation (Fig. 1b).

Over the entire time course, a total of 2741 unique phosphopeptides were differentially phosphorylated for at least one time point, accounting for about 21% of all phosphopeptides analyzed. Interestingly, the overlap of differentially phosphorylated sites between the particular time points was rather moderate, reaching approx. 24–28% (Supplementary Fig. 3). In addition, we identified phosphosites that were exclusively regulated at one time point only, such as after one minute (262 sites) or 60 min (497 sites) of insulin stimulation, respectively (Supplementary Fig. 3). Most regulated phosphopeptides (approx. 75%) contained a single phosphorylated Ser/Thr residue, a sizable fraction of regulated peptides, however, was phosphorylated at multiple Ser/Thr residues (Fig. 1d).

Figure 1e–g illustrates the dynamics in the abundance of phosphopeptides for each time point relative to the basal state. Several phosphopeptides were not found in basal cells but readily detectable after insulin stimulation, such as peptides derived from the glucocorticoid receptor (GR) containing S134, and C-Jun-amino-terminal kinase-interacting protein 4 (JIP4) containing T595. Furthermore, at each time point, insulin stimulation resulted in both increases and decreases in the abundance of specific phosphopeptides when compared to the basal value, whereas the distribution of fold changes remained similar over the time course. After 1 min of insulin stimulation, 822 and 749 phosphopeptides were up- and downregulated relative to baseline, and after 60 min, 1048 and 720 phosphopeptides were up- and downregulated, respectively.

Donor and time-dependent clustering of insulin-regulated phosphopeptides

In further analyses, we considered only phosphopeptides that had been quantified at each time point for each donor, resulting in a set of 11,612 phosphopeptides with continuous quantitative information over the entire time course without missing values. We then investigated the individual response to insulin for the different donors over time using principal component analysis (PCA). When all quantified phosphopeptides were considered, the samples from individual donors showed the highest degree of similarity, mostly independent from the time of insulin treatment (Supplementary Fig. 7). In contrast, PCA of the most distinctively insulin-regulated phosphopeptides (p < 0.01 and FC > 3 or <¹/₃ vs. t₀; with continuous quantitative information for each time point and donor) revealed distinct clusters dependent of the duration of insulin stimulation: A cluster comprising phosphopeptides in the basal state (t₀), a cluster corresponding to 1 and 2.5 min of insulin stimulation (t₁ and t_2.5), a cluster after 30 min and 60 min (t₃₀ and t₆₀) of insulin stimulation and an intermediate cluster corresponding to 5 min and 15 min (t₅ and t₁₅) respectively (Fig. 2a). We therefore designated the clusters as reflecting early (t₁, t_2.5), intermediate (t₅, t₁₅), and late (t₃₀, t₆₀) phosphorylation events in insulin signaling. We next analyzed the overlap of significantly regulated insulin targets across the different clusters. As shown in Fig. 2b, a substantial fraction of phosphopeptides was unique for each of the clusters while phosphorylation of >400 insulin-stimulated sites persist over the entire time course.

Fig. 2: Donor and time-dependent clustering of peptides and proteins in human myotubes.

a Principal component analysis (PCA) of most regulated phosphopeptides for each donor and time point. b Time interval-specific and interval-shared phosphopeptides regulated by insulin. Intervals include “Early” (1 and 2.5 min), “Intermediate” (5 and 15 min), and “Late” (30 and 60 min) time points of insulin stimulation. Source data are provided as a Source Data file.

Targets of canonical insulin signaling in human skeletal myotubes

Next, we aimed at dissecting the time course of phosphorylation of the canonical insulin signaling pathway and selected targets of AKT, such as p70 ribosomal S6 kinase (p70S6K/RPS6KA/B) that is activated through the phosphatidylinositol 3-kinase (PI3K)-AKT-mammalian target of rapamycin complex 1 (mTORC1) pathway and p90 ribosomal S6 kinase (p90RSK/RPS6KB) that is activated through Ras-Raf-extracellular signal-regulated kinases (ERKs), MAPK1 and MAPK3 (Fig. 3a).

Fig. 3: Time course and donor-to-donor variability of insulin-stimulated phosphorylation for selected insulin targets in human myotubes.

a Schematic insulin signaling of established insulin-regulated kinases, substrates, and main regulatory phosphorylation sites. Quantitated phosphosites for AKT are labeled in red, and identified sites without quantification shown in gray. b Phosphorylation kinetics based on quantification of phosphopeptides for 28 selected targets. Abundance values for individual donors (symbols) and mean values (red lines) are shown. The individual abundance values were used to calculate the average Pearson’s correlation coefficient r as a measure of variability between donors. Targets with low variability (r < 0.8) are indicated by blue color in (a). c Representative Western blot analysis (n = 3 donors) of selected insulin-regulated phosphosites. Source data are provided as a Source Data file.

As illustrated in Fig. 3b, insulin stimulation resulted in rapid phosphorylation of both regulatory sites, T309 and S474 in AKT2, with maximal phosphorylation achieved after 5 min. Several AKT substrates including AKT1S1 (T246), FOXO3 (S253), GSK3B (S9), and MYO5A (S1652) were phosphorylated with similar kinetics whereas others showed a more delayed phosphorylation such as ACLY (S455), PDCD4 (S76), and TBC1D4 (T642). Phosphorylation of MAPK1 (T185) also peaked at 5 min and preceded that of MAPK3 (T202) and RPS6KA3 (S369), consistent with delayed phosphorylation of downstream targets such as BAD (S75), eIF4B (S422) and RPS6 (S236). Phosphorylation of mTOR (S2448), RICTOR (T1135), and p90RSK (S369) reached maximal levels after 30 min, whereas phosphorylation of ACLY (S455), BAD (S75), and eIF4B (S422) continued to increase further. Elevated phosphorylation levels were mostly maintained after 60 min of insulin exposure. Interestingly, insulin stimulation led to a rapid reduction in phosphorylation of S50, a major activating site in the tyrosine phosphatase PTPN1¹² (Fig. 3b). Western blot analysis was used to validate phosphorylation sites and kinetics of insulin target sites (Fig. 3c, Supplementary Fig. 8).

We then analyzed the response to insulin for each of the five donors separately, as well as their relation to each other by calculating the average Pearson’s correlation coefficient (Pearson’s r) between all pairs of donors in order to assess the individual variability of phosphorylation kinetics over time (Fig. 3b). In fact, variation in the phosphorylation kinetics for all five donors was moderate and amounted to a median Pearson’s r of 0.82 for all 28 selected phosphosites. The correlation was substantially higher for annotated AKT sites (median r = 0.914) with T246 in AKT1S1 showing the least inter-donor variability among the selected substrates (r = 0.982; Fig. 3b).

Donor variability of insulin targets

As Pearson’s r may reflect the propensity of a specific Ser/Thr residue to be directly phosphorylated by insulin-activated kinases, we calculated correlation coefficients for all of the 8301 monophosphorylated peptides in our dataset quantified for all donors at all time points. Pairwise correlation analysis (i.e. per phosphosite r values for all 10 different donor pairs were computed using their respective phosphorylation kinetics) revealed several phosphosites with low donor variability, i.e. mean r > 0.9 (Fig. 4a). AKT1S1-T246, GAB2-S210, and NACA-S1112 were the three phosphosites with the lowest variability in our dataset with r > 0.96. On the other hand, EVA1B-T158, TLE4-S292, and ADD3-S677 were the phosphosites with the highest donor variability with r < −0.218. The median of all 8,301 average r values was 0.31. In fact, we observed negative correlations to be less frequent than positive correlations in our 83,010 pairwise donor comparisons of phosphorylation kinetics. Interestingly, some AKT target sites like T246 in AKT1S1, S210 in GAB2, and S253 in FOXO3 had lower variability in their phosphorylation profile than the two regulatory phosphorylation sites in AKT2, T309, and S474. In addition to well-established insulin targets, we identified several other insulin target sites with very low donor variability in proteins like NACA, AFDN, and KHDRBS1 annotated for a variety of biological processes such as transcriptional regulation, mRNA processing and protein trafficking (Fig.4b).

Fig. 4: Donor variability of insulin-regulated phosphosites.

a Abundance values for monophosphorylated peptides from each of the five donors were used to calculate pairwise (empty circles) Pearson’s correlation coefficient r and mean r values (filled circles) as a measure for donor variability. b Time course of selected insulin targets with low donor variability and respective Gene Ontology annotations (GO:Biological process, BP). Symbols denote single donor values, red lines represent mean values. Source data are provided as a Source Data file.

We compared the phosphorylation of canonical and non-canonical insulin targets with phosphosites associated with cellular stress response¹³ and type 2 diabetes⁹. While the majority of stressors are linked to the phosphorylation of canonical insulin targets and mTOR signaling in L6 myotubes, human iPSC-derived myoblasts (iMyos) from T2D donors also showed alterations in insulin-stimulated phosphorylation of some non-canonical targets including ARHGEF12, NEK9, and PLEKHF2 as shown in Supplementary Fig. 9. In addition, we used data from a recent large-scale multi-ancestry meta-analysis¹⁴ to investigate the overlap of insulin target proteins with GWAS genes associated with T2D. Our analysis revealed 43 proteins, mostly non-canonical targets, that are differentially phosphorylated in response to insulin and coded by T2D candidate genes (Supplementary Data 1).

Time-resolved regulation of Ser/Thr kinases

In response to insulin stimulation, we found 49 protein kinases differentially phosphorylated, among those several reported to be involved in intracellular signaling pathways regulating cytoskeleton dynamics, vesicle transport, metabolism, apoptosis, mitosis, and transcriptional control (Supplementary Data 2). In order to detail temporal aspects of kinase signaling, we conducted kinase-substrate overrepresentation analysis (KSORA) as described in “Methods” section. Briefly, changes in phosphopeptide abundance over time (t_{1 …}t₆₀ versus t₀) were used to compute the respective kinase activation based on assignments of 13,664 annotated phosphosites and 403 unique human protein kinases in the PhosphoSitePlus database¹⁵. Our analysis revealed a distinct temporal pattern of kinase activation where for each time point, a different set of Ser/Thr kinases was identified to be activated (Fig. 5a). According to the kinase-substrate overrepresentation analysis, insulin stimulation led to activation of a total of 159 Ser/Thr kinases of which the most significant 19 are shown in Fig. 5a. At early time points (t₁…t_2.5), enriched kinases are annotated for regulating cell cycle and differentiation (e.g. PLK1, MARK2, CDK6). Intermediate time points (t₅…t₁₅) showed enrichment of well-known kinases of the insulin signaling pathway, such as AKTs upstream activating kinases PDPK1/PDK1, and a diverse set of kinases implicated as regulators of the cytoskeleton (e.g. OXSR1, MAPKAPK2), cell cycle and mitosis (AURKB, CDK1, CDK5), and metabolism (e.g. mTOR). Kinases enriched at late time points (t₃₀…t₆₀) include p70 and p90 ribosomal S6 kinases RPS6KA/RPS6KB, extracellular signal-regulated kinases MAPK1/MAPK3, MAPK14, and other kinases such PKCs, SGKs involved in metabolic actions, transcriptional control and skeletal muscle homeostasis.

Fig. 5: Time course of kinase-substrate overrepresentation after insulin stimulation of human myotubes.

Kinase-substrate overrepresentation analysis (KSORA) with 4030 differentially regulated phosphopeptides was performed as described in “Methods” section. a Heat map of the top 19 overrepresented kinases (one-sided Fisher’s exact test p-value ≤ 0.075, not adjusted) across all time points where colors represent z-scores of the respective −log10(p-values) of the KSORA. Asterisks denote p-values *<0.075; *<0.05; **<0.01; ***<0.001. b Corresponding numbers of substrates over time that are annotated for the most highly enriched kinases. Source data are provided as a Source Data file.

A varying number of substrates were regulated at different time points of insulin stimulation with AKT, as expected, representing the most dominant kinase among the panel (Fig. 5b) with up to 27 annotated substrates at t₃₀. The >500 known human protein kinases have been classified into 7 major groups, based on sequence comparisons of their catalytic domains¹⁶. Mapping of the enriched kinases to the phylogenetic kinome tree indicated that insulin stimulation in skeletal muscle cells results in group-specific activation of kinases, with AGC family members being the most predominant category (Supplementary Fig. 10).

Network propagation and subnetwork identification

We further aimed to structure the insulin-regulated phosphoproteome by applying the method of network propagation which has been successfully used previously in identifying critical nodes and functional modules of genes and proteins controlling complex traits¹⁷. Briefly, we utilized the “NetCore” algorithm that, in an iterative manner, superimposes protein–protein interaction (PPI)-based networks on nodes in the network consisting of kinases, phosphatases, and corresponding substrates while considering network topology features to eventually derive phosphorylation-dependent relationships of proteins¹⁸.

We used a filtered subset of 1552 significantly regulated phosphopeptides with continuous quantitative information for each time point (fold change versus t₀ for t₁–t₆₀) to initialize the PPI network and conducted network propagation utilizing 10,642 proteins and 139,266 interactions for each time point separately as described in “Methods” section. Each of the resulting networks was finally merged into a combined PPI network consisting of 220 nodes and 582 connections, representing first and second-degree neighbors/interaction partners that share similar temporal dynamics of phosphorylation in response to insulin stimulation (Fig. 6). Further functional classification identified a signaling module of 14 Ser/Thr protein kinases with 57 first-degree neighbors and 4 protein phosphatases connected to 47 first-degree neighbors. Both kinases and phosphatases are linked through a set of 15 shared first-degree neighbors that include adapter proteins (YWHAH, GRB2), regulators of transcription and splicing (e.g. SNW1, SNRNP70, SRSF1, ZC3H18) and actin cytoskeleton (CAPZA1, LIMA1, FLNA). Moreover, 5 kinases were identified as first-degree neighbors of phosphatases, including MAP3K3, MAPK1, MAPK3, NTRK1, and PRKD1 propagated by NetCore. Conversely, the two phosphatases PPP1CC and PPP1CA, the latter propagated by NetCore, were identified as first-degree neighbors of kinases. About half of the kinases inferred by NetCore were also identified by KSORA as insulin targets, including AKT2, MAPK1, MAPK3, MAPKAPK2, MYLK, and PDPK1.

Fig. 6: Donor- and time-dependent protein–protein interaction network of insulin-regulated kinases and phosphatases.

a Proteins with significantly regulated phosphopeptides in response to insulin (n = 1552; t-test p-values < 0.05, adjusted using the Benjamini–Hochberg method) were selected for network propagation and assembled into a combined PPI network using NetCore and ConsensusPathDB as described in “Methods” section. The resulting network of 220 nodes and 582 edges reflects annotated connections and similar temporal dynamics of phosphorylation in response to insulin stimulation. Indicated are protein kinases (green) and phosphatases (orange), and their respective first-degree protein–protein interaction neighbors (light colors). Diamonds denote proteins propagated from NetCore analysis. b Subnetwork of phosphoproteins with low donor variability and similar phosphorylation kinetics. Nodes represent proteins with the lowest inter-donor variability (mean of Pearson’s r between all donors >0.78), and edges indicate similar phosphorylation kinetics of phosphopeptides with low variability (mean of Pearson’s r between all donors >0.73). Shown are protein kinases (green) and phosphatases (orange), and their respective first-degree protein–protein interaction neighbors (light colors). RNA-binding proteins are marked in red. Dotted lines connect proteins with highly similar phosphorylation kinetics (Pearson’s r > 0.9). c Correlation matrix of Perason’s r values for edges (phosphorylation kinetics) between NetCore-connected phosphopeptides with low donor variability. Source data are provided as a Source Data file.

We further extracted a subnetwork of 18 phosphoproteins with the lowest donor variability between nodes (Pearson’s r > 0.78) and along edges (Pearson’s r > 0.73), reflecting high connectivity within the network interaction and similarity in phosphorylation kinetics (Fig. 6b, c). As expected, all nodes directly connected by an edge in this subnetwork exhibit similar phosphorylation kinetics, both on average and in donor-to-donor comparisons (Supplementary Fig. 11). Several proteins of the Reactome mTOR signaling pathway (AKT1S1, AKT2, RPS6, EIF4B) were significantly enriched (p = 7.39 × 10^-6) in this subnetwork, in addition to RNA-binding proteins involved in processing and splicing of mRNA (EIF4B, KHDRBS1, MATR3, RMB34, SNRNP70, THRAP3, ZFP36L1) or rRNA (LYAR), respectively.

Temporal regulation of biological pathways

For phosphopeptides significantly regulated at early (t₁, t_2.5), intermediate (t₅, t₁₅), and late (t₃₀, t₆₀) time points, we conducted functional enrichment analysis using the Reactome pathway database. Briefly, for each time interval, differentially regulated phosphopeptides were filtered (⎸FC⎹ > 1.3; p < 0.05) and mapped to a non-redundant set of proteins, which was used to conduct overrepresentation analysis as described in “Methods” section. Hierarchical clustering revealed a distinct temporal pattern of significantly enriched pathways related to growth factor signaling, cytoskeleton and muscle contraction, cell development, AKT signaling, protein trafficking, and mRNA processing and splicing (Fig. 7a and Supplementary Data 3).

Fig. 7: Time course and donor variability of overrepresented Reactome pathways after insulin stimulation of human myotubes.

a Reactome pathways overrepresentation analysis (RORA) was performed for three time intervals, early, intermediate, and late, corresponding to 1 and 2.5 min, 5 min and 15 min, and 30 and 60 min after insulin stimulation with 969 differentially phosphorylated proteins as described in “Methods” section. The heat map colors represent z-scores of the respective -log10(FDR) values of the RORAs to visualize differential phosphorylation within the respective pathway. b Numbers of differentially regulated phosphopeptides for selected Reactome pathways. The red and blue bars represent phosphopeptides having positive or negative log-fold changes versus baseline, respectively. c Donor variability of pathway enrichment. Pearson’s correlation coefficient r of abundance values were calculated for the differentially regulated phosphoproteins within each pathway, and the values displayed as a measure for variability between donors. Red points represent the mean values. Source data are provided as a Source Data file.

Interestingly, pathways associated with cell cycle control were enriched at very early time points (59 differentially phosphorylated proteins) with rapid phosphorylation of the mitotic protein kinase NEK9, proteins associated with the nuclear lamina (e.g. EMD, LMNA), and the nuclear pore complex (e.g. TPR, NUP188, RANBP2) or nuclear transport (e.g. SUMO1, NPM). Pathways associated with tyrosine kinase signaling (65 differentially phosphorylated proteins) were enriched at intermediate and late time points and include known docking proteins in signaling such as CRK, GAB2, IRS1, IRS2, SH2B2, SHC1, SOS1, mitogenic kinases like MAPK1, MAPK3 and MAPK14, and transcriptional regulators like MEF2C and YAP1. Conversely, the pathway for vesicle-mediated transport (74 differentially phosphorylated proteins) was enriched at late time points, including proteins involved in clathrin-dependent transport (e.g. AAK1, AP3B1), GEFs and GAPs of Rabs (e.g. DENND1A, RAB3GAP1, TBC1D2) and Arf GTPases (ARFGAP1, GBF1) and various steps in intra-organelle membrane traffic (e.g. KIF1B, KIF13B, SEC16A, SEC22B, SNX9, TRAPPC12, VAMP3).

In most pathways, insulin led to both increased or decreased phosphorylation of proteins compared to the basal state. We further investigated the dynamics of pathway activation by segregating phosphoproteins with either increasing or decreasing phosphorylation in response to insulin within each pathway. As shown in Fig. 7b, insulin increased the number of phosphorylated peptides over time for “mTOR signaling” (R-HSA-165159) without such changes in the pathway for “Muscle contraction” (R-HSA-397014). Conversely, insulin reduced the number of phosphopeptides with decreased phosphorylation over time compared to baseline in “Muscle contraction” but not in the “mTOR signaling” pathway. Other pathways, like “Membrane trafficking” (R-HSA-199991), displayed an intermediate pattern where the number of phosphorylated peptides increased over time and, in parallel, less peptides were found that had decreased phosphorylation compared to the basal state (Fig. 7b). Collectively, the data indicate complex temporal segregation of insulin signaling within the different biological pathways. We also analyzed donor-to-donor variability in pathway activation by conducting pairwise correlation analysis of phosphosites assigned to proteins of the enriched Reactome pathways (Fig. 7c). While activation of pathways related to mTOR and AKT signaling displayed the lowest donor variability (Pearsons’s r > 0.8), TGF beta signaling and signaling of non-receptor tyrosine kinases exhibited high donor variability (Pearsons’s r > 0.8) (Fig. 7c).

Temporal dynamics of relative phosphosite occupancy

We further investigated the temporal dynamics of relative phosphosite occupancy during insulin stimulation by tracking the phosphorylation of individual phosphopeptides over time. For the analysis, we selected a subset of 2741 phosphopeptides that have a significant regulation with respect to baseline (0 min) for at least one time point from the 11,612 phosphopeptides with continuous quantitative information for each of the 7 time points. We then assessed the phosphorylation of each phosphopeptide over the entire time course and aggregated the data in a flow diagram. As illustrated in Fig. 8a, insulin stimulation results in a complex pattern of phosphorylation and dephosphorylation events over time, where the former were more common than the latter (Supplementary Fig. 12). The individual rates of phosphorylation and dephosphorylation showed complementary, interval-specific patterns for each time point with highest dynamics early after insulin stimulation (Fig. 8b).

Fig. 8: Temporal dynamics of phosphopeptide regulation.

a Flow diagram of relative phosphosite occupancy for 2741 unique phosphopeptides. Each phosphopeptide is represented by a single line where the colors denote the initial time point of continuous differential phosphorylation in response to insulin (⎸FC⎹ > 1.3 vs. baseline; p < 0.05). Lines aggregate continuous and discontinuous differential phosphorylation of individual peptides over time. Numbers in boxes indicate the numbers of unregulated (top level) and differentially phosphorylated peptides for each time point. b Time-dependent relative abundances (normalized to maximum values) for phosphorylation and dephosphorylation of peptides from the strands in (a) using the same color code regarding the initial time point of differential phosphorylation. Strand-specific numbers of peptides: n₁ = 111, n_2.5 = 28, n₅ = 53, n₁₅ = 98, n₃₀ = 128, n₆₀ = 345 (phosphorylation); n₁ = 20, n_2.5 = 7, n₅ = 11, n₁₅ = 17, n₃₀ = 57, n₆₀ = 258 (dephosphorylation). In the box plots, median (horizontal line), minimum (lower whisker), maximum (upper whisker), 25th and 75th percentile (lower and upper hinge) values of the respective phosphopeptide abundances are visualized. c Phosphosite motif analysis for 206 targets phosphorylated at 5 min; of those, 49 sites continuously phosphorylated until 60 min; and 197 targets dephosphorylated at 5 min. The frequency of residues surrounding the phosphorylated residue is indicated by the size of the letters, as described in “Methods” section; d Time-resolved relative phosphosite occupancy of selected Ser/Thr residues in proteins significantly overrepresented in the Reactome mTOR signaling pathway (Supplementary Data 4). Source data are provided as a Source Data file.

Of the 976 phosphopeptides significantly regulated after 1 min of insulin stimulation, only 131 of those phosphopeptides (~13%) remained continuously differentially phosphorylated for 60 min. Notably, about a third of these phosphosites were in proteins associated with the cytoskeleton and cytoskeleton organization (GO:0005856, GO:0007010; 40 proteins, p-adj = 1.08e-08) and intracellular signal transduction (GO:0035556; 36 proteins, p-adj = 0.000201), including important upstream kinases for AKT, PDPK1 and mTOR (Supplementary Table 3). Interestingly, several of the upregulated phosphopeptides (10 out of 110) were derived from proteins implicated in insulin-stimulated GLUT4 traffic, such as BIN1, DENND4C, ESYT1, MYO5A, PHLDB1, RALGAPA2, RAPGEF1/C3G, STX4, STX16, and TXNIP, presumably also phosphorylated by kinases upstream of AKT (Supplementary Table 4). Conversely, of the 1133 significantly regulated phosphopeptides at 60 min, a large proportion of 603 peptides (53%) was derived from phosphorylation/dephosphorylation events that occurred at 60 min. Interestingly, these late phosphorylation sites targeted a large number of RNA-binding proteins (GO:0003723, 113 proteins, p-adj = 1.85e-26), with 398 of the 1133 differentially phosphorylated proteins annotated for nuclear localization (GO:0005634, p-adj = 6.22e-22).

Analysis of phosphorylated sequence motifs of all phosphopeptides in the subset revealed interval-specific differences during insulin stimulation. The basophilic consensus motif of AKT targets, (RxRxxS/T)⁵, was more common in phosphopeptides exhibiting continuous phosphosite occupancy, whereas peptides dephosphorylated in response to insulin appeared to be enriched for proline-directed consensus sequences (Fig. 8c). Moreover, proline-directed/acidic motifs were more enriched in phosphosites at very early time points (Supplementary Fig. 13). We further investigated the temporal coordination of relative phosphosite occupancy using Gene Ontology and Reactome databases with our subset of phosphopeptides with continuous quantitative information for each of the time points (Supplementary Fig. 14). Differential phosphorylation of proteins assigned to pathways was complex and interval specific. For instance, mTOR signaling was accompanied by early, intermediate, and late phosphorylation events occurring over the whole time course of insulin stimulation, in addition to continuous phosphorylation of mTOR S2248, and the majority of other proteins being phosphorylated between 5 and 15 min (Fig. 8d). Similar patterns of activation were found for other pathways like Ras signaling (Supplementary Fig. 15).

Insulin-stimulated phosphorylation of splicing-related proteins

A substantial number of insulin-regulated phosphopeptides were derived from proteins annotated for processing of mRNA at various stages. For instance, we found differentially regulated phosphopeptides mapped to 49 out of the 191 proteins within the Reactome pathway “mRNA Splicing” (R-HSA-72172.3). In accordance, recent phosphoproteomics studies also reported insulin-regulated phosphorylation of spliceosome-associated proteins in iPS-derived myoblasts but without information on the temporal coordination of the process⁹. We, therefore, investigated more closely the time-dependent phosphorylation of spliceosomal proteins in primary human myotubes.

As illustrated in Fig. 9a, mRNA splicing, i.e. the removal of introns from precursor mRNA, occurs through a series of highly coordinated steps carried out by the spliceosome, a very large RNA–protein complex that assembles anew on each individual pre-mRNA substrate to specific intermediate complexes¹⁹. While the spliceosome is composed of five small nuclear ribonucleoproteins (snRNPs), U1, U2, U4, U5, and U6 snRNPs¹⁹, and several spliceosome-associated proteins (SAPs), it dynamically assembles on the mRNA precursor in several, functionally distinct complexes (termed E, A, B, B*, and C) that altogether encompass the more than 120 different proteins involved in the different steps of mRNA splicing¹⁹. Our analysis revealed numerous spliceosomal proteins that are differentially phosphorylated in response to insulin, including SNRNP70, SF3A1, SNRNP200, and SART1, that constitute components of the U1, U2, U5, and U4/U6-U5 tri-snRNP, respectively and are required for the assembly of the pre-catalytic spliceosome complex A and B (Fig. 9a, b; Supplementary Fig. 16). SRSF protein kinase 2, a regulator of RNA-binding SR proteins, splicing factor SRSF1, SRSF2 and associated mRNA binding proteins SRRM1 and SRRM2, U2AF2, RBFOX2, and RMB10, all found to be involved in either constitutive or alternative splicing, were also phosphorylated or dephosphorylated after stimulation with insulin. Interestingly, the temporal pattern of phosphosite occupancy within spliceosomal proteins exhibited substantial variation during insulin stimulation. While some sites showed a rapid and transient increase in abundance (e.g. S207 in SNRNP200; t_max 1 min), others increased more steadily (e.g. S448 in SART1) or decreased rapidly (S154 in RBM10) in response to insulin.

Fig. 9: Phosphorylation of the spliceosome.

a Schematic overview of the splicing process steps consisting of intron/exon recognition, ligation (complexes A, B), catalytic activation (complex C), and cleavage and exon ligation step. b Phosphorylation kinetics of spliceosomal phosphosites in 12 selected spliceosomal proteins, and their assignment to the splicing complexes indicated by colors. Data are mean values ± SEM (n = 5). c Distribution of different mRNA splicing events in myotubes from 3 donors after 1 h stimulation with insulin. d Insulin-regulated alternative splicing of transcripts. The percent spliced in index (PSI; n = 3 ± SEM) refers to the ratio between reads including or excluding exons (denoted by asterisks) as described in “Methods” section. A total of 153 transcripts were found to be differentially spliced in response to insulin (Supplementary Data 5). Source data are provided as a Source Data file.

We further investigated a possible association of insulin-stimulated phosphorylation of spliceosomal proteins with insulin-regulated alternative splicing using the same cells from our donors that were used for phosphoproteome analysis. After starvation for 6 h, fully differentiated primary myotubes from 3 different donors were incubated with and without 100 nM insulin for 1 h, harvested, and the resulting RNA was subjected to Next Generation Sequencing as described in “Methods” section. Differential expression at the whole transcriptome level was only moderately changed, affecting 21 genes (20/1 up/downregulated), some of which are implicated in different metabolic and biological responses such as regulation of insulin sensitivity (FOXC2, NR4A1, EGR1), cell growth (FOS, GADD45B) and mitochondrial biogenesis (PPRC1) (Supplementary Table 5). We applied multivariate analysis of transcript splicing (MATS) to identify transcripts that were differentially spliced in response to insulin. We found 153 transcripts for 131 coding genes and 9 long non-coding RNAs (lncRNAs) that were differentially spliced in response to insulin (Supplementary Data 5), and, according to principal component analysis, independent of the donor (Supplementary Fig. 17b). The majority of detected splicing events included skipped exons (SE), followed by alternative 3′ splicing (A3SS), mutually exclusive exons (MXE), retained introns (RI) and 5′ splicing (A5SS; Fig. 9c). The splicing of several mRNAs (e.g. UDP-glucose ceramide glucosyltransferase, UGCG) is predicted to result in nonsense-mediated decay (NMD; Fig. 9d; Supplementary Data 5).

Discussion

In the present study, we conducted a time-resolved analysis of the global insulin-regulated phosphoproteome of human primary skeletal muscle cells. Our results identify a complex set of temporally regulated networks and their donor-to-donor variation and provide novel insights into the impact of insulin on mRNA processing and splicing in human myotubes.

Insulin action in skeletal muscle is critical for whole-body energy substrate metabolism and homeostasis, and impaired insulin signaling in muscle has been associated with early onset and development of type 2 diabetes³. Insulin induces rapid and reversible phosphorylation of target proteins, but the complex network of its signaling cascade in muscle is not sufficiently understood. Previous studies have investigated insulin-regulated muscle phosphoproteomes at a fixed time point after hormone stimulation using rat L6 myotubes (15 min; 100 nM insulin), human immortalized muscle precursor cells (10 min; 100 nM insulin), and biopsies from human skeletal muscle tissue (2 h; hyperinsulinemic-euglycemic clamp)^8,9,10. Consistent with these studies, our phosphoproteomics analysis of human primary myotubes revealed a large number of insulin-regulated phosphopeptides, reflecting the phosphorylation status of >11,500 unique Ser/Thr phosphosites over the course of 1 h in the presence of the hormone. However, a few sites (e.g. TBC1D4 S570), were phosphorylated in both starving and insulin-stimulated conditions, presumably reflecting activation of AKT and AMPK, two distinct metabolically controlled kinases²⁰. Importantly, our data demonstrate strong temporal dynamics in insulin-regulated phosphorylation without corresponding changes in protein abundance. While some sites were phosphorylated transiently only at very early time points (e.g. SNRNP200 S207) others showed multimodal phosphorylation patterns (e.g. NACA-S1112) or had phosphorylation levels reaching a plateau (e.g. AKT1S1-T246). In fact, a substantial fraction of phosphopeptides, approximately 30%, is differentially regulated only at specific intervals during the course of insulin stimulation, and approx. 1400 regulated phosphosites would not have been observed when analyzing only phosphopeptides after 15 min of insulin stimulation. Thus, time-resolved analysis of the phosphoproteome provides enhanced sensitivity and resolution of phosphorylation events and may aid in identifying novel aspects of insulin signaling in muscle cells. For instance, while recent genome-wide association studies (GWAS) have identified numerous loci comprising >600 human protein-coding candidate genes for insulin resistance and T2D¹⁴, little is known about their potential role as targets for insulin-regulated phosphorylation. We explored the overlap of our phosphorylation data and found 43, mostly non-canonical target proteins, which may be of interest for further studies on insulin action in skeletal muscle cells²⁰.

As expected, substrates for AKT kinase are overrepresented in the phosphoproteome after insulin stimulation. In addition, several other kinases assigned to multiple Ser/Thr kinase phyla appear to be regulated at distinct time points, either preceding or following the activation of AKT. Among the kinases whose substrates are enriched relatively early in the time course (<5 min) are PDPK1, a well-known upstream activating kinase of AKT²¹, as well as cell cycle regulating kinases CDK1/5/6 and PLK1, whereas substrates for p70 and p90 ribosomal S6 kinases are enriched at later time points. Thus, consistent with a previous kinase-substrate prediction analysis for AKT, mTOR, and PKA in 3T3-L1 adipocytes²², we found interval-specific phosphorylation patterns of kinases and respective substrates in primary human myotubes. Interestingly, our enrichment analysis reveals a distinct temporal pattern of kinase activation with phase lengths in the order of minutes, resulting in little to moderate overlap between the different time points after insulin stimulation. This transient kinase activation pattern in the presence of insulin implies persistent activity of protein phosphatases, continuously reverting phosphorylation by the kinases, as well as kinase-driven feedback loops. While previous studies found phosphorylation-dependent feedback regulation in insulin signaling at multiple levels, including insulin receptor, IRS1, and AKT signaling via S6K1-mediated phosphorylation of PDK1, relatively little is known about the complex regulatory network of protein phosphatases in insulin signaling^12,23,24,25.

Following a kinase/substrate annotation-independent approach, we analyzed potential interdependencies of insulin signaling components by using the propagation algorithm “NetCore” that utilizes PPI data to segregate and identify functional networks within large and dynamic datasets and considers node coreness as a topological feature, i.e. whether a node belongs to a densely connected part of the network or its periphery¹⁸. This unsupervised method identified a signaling module encompassing 220 members, including 15 protein kinases and 4 phosphatases and their respective first-degree neighbors, some of which share similar temporal dynamics of phosphorylation in response to insulin stimulation. Module-associated phosphatases include PPP1CA and PPP1CC, catalytic subunits of PP1, a regulator of glycogen synthase²⁶, substrate specificity conferring regulatory subunits like the myosin targeting subunit PPP1R12B, and Calcineurin (PPP3CC), a Ca²⁺/Calmodulin-regulated Ser/Thr protein phosphatase shown to be involved in regulating insulin signaling proteins and insulin-stimulated glucose uptake in skeletal muscle²⁷. Moreover, the module includes CDC14B, which preferably dephosphorylates targets of proline-directed Ser/Thr kinases and has been recently implicated in regulating ciliogenesis^28,29. Interestingly, 52 nodes, about a quarter of the network, were first-degree neighbors of four phosphatases, whereas 72 nodes were exclusively connected to 15 kinases, emphasizing the different levels of integration of kinases and phosphatases in the signaling network. However, incomplete kinase/substrate and phosphatase/substrate annotation (“dark proteome”) and overlapping substrate specificities of kinases represent major caveats in understanding the complexity of insulin signaling³⁰. In fact, despite using large data repositories, the majority of our quantified phosphosites could not be unequivocally associated with a specific kinase.

Using human primary cells from matched healthy donors allowed us to address individual variability in insulin signaling, albeit with the caveat of analyzing a relatively low number of donor samples. Analysis of donor variability within the PPI-based network and correlation of phosphorylation kinetics revealed a subnetwork of proteins with similar phosphorylation patterns enriched in proteins involved in RNA processing, indicating insulin’s regulatory effect on transcription at multiple levels. While some connections between proteins were consistent with PPI-based associations, several highly correlated phosphosites may reveal novel associations and need further exploration. A comparison of site-specific kinetics revealed variability across donors in human skeletal myotubes. Donor variability was relatively low for phosphosites belonging to the canonical insulin signaling pathway with Pearson’s r > 0.9, such as AKT1S1 (T246), FOXO3 (S253), and TSC2 (S939). Consequently, donor variability of Reactome pathway regulation was low for AKT-driven processes, such as “mTOR signaling” (R-HSA-165159) and GLUT4 translocation pathway (R-HSA-1445148). This result is consistent with previous studies that have shown that co-occurring phosphorylation and correlated phosphorylation levels between different phosphosites indicate their functional association^31,32. Proteins that are part of a protein–protein complex participate in the same biological processes and may become co-phosphorylated, presumably by the same kinase^31,32, thus reducing variability in phosphorylation kinetics. Low donor variability represents another aspect of phospho-kinetic similarity that reflects highly efficient target sequences (e.g. in AKT1S1). Nodes with phosphorylation kinetics that are robust to individual differences of healthy donors could be particularly important switching points in biological signaling networks. Hence, the deviating kinetics of such phosphosites may identify novel, biologically relevant insulin targets and might be used to exploit future precision medicine approaches. Conversely, phosphosites exhibiting substantial donor variability are less likely to contribute to altered systemic insulin action or serve as informative biomarkers for insulin sensitivity³³. Interestingly, several substrates with very low donor variability constitute RNA-binding proteins that were not previously assigned to insulin signaling pathways, indicative of insulin’s possible role in regulating RNA processing and splicing.

We investigated the dynamics of relative phosphosite occupancy and pathway regulation for a subset of phosphopeptides with continuous quantitative information for each time point. Interestingly, while many phosphorylations were transient, even in the presence of insulin, a fraction of sites were differentially phosphorylated for longer periods of time. More than 100 sites were phosphorylated within 1 min in response to insulin stimulation and maintained their phosphorylation status for up to 1 h. Of those, many proteins were associated with the cytoskeleton and functions relayed by small GTP-binding proteins, indicative of insulin’s fundamental role in regulating cell structure and subcellular organization. Conversely, more than 340 of the sites phosphorylated after 1 min apparently became dephosphorylated thereafter, most of them within the following 5 min during the time course.

Kinetic analysis of the individual rates and consensus sequences of phosphorylated and dephosphorylated peptides at each interval substantiate time-dependent preferences for kinases and/or phosphatases, consistent with phased regulation of the enzymes. Sites phosphorylated early and maintained for longer intervals become enriched for more basophilic motifs, particularly the AKT target motif RxRxxS/T and RxxS/T, the latter being also a potential target for many AGC family kinases from such as PKA, PKC, RSK, ROCK, and multifunctional calmodulin-dependent protein kinase II³⁴. Conversely, early dephosphorylated sites are enriched for Proline at position -1, preferred sites for Pro-directed kinases such as conventional and atypical MAPKs and Cyclin-Dependent Kinases³⁴, indicating consecutive deactivation of these kinases over time in the presence of insulin and/or activation of phosphatases targeting proline-containing motifs, such as CDC14B, identified by our Netcore analysis.

Consistent with the distinct temporal pattern of kinase activation, we found a similar time-dependent relationship in the regulation of biological pathways as annotated in the Reactome database. Transient regulation of pathways such as cell cycle control occurred shortly after stimulation with insulin whereas pathways related to e.g. mRNA processing and splicing became differentially regulated after prolonged stimulation with the hormone.

Interestingly, pathways such as mTOR signaling appear to be activated in a phased fashion where early onset insulin phosphosites are dephosphorylated over time to be replaced with phosphosites in other proteins annotated for this pathway, while only mTOR phosphosites are continuously phosphorylated during the time course. Similar patterns were found for other regulated pathways, indicating that pathway components, some of which are shared among different pathways, contribute in a time-synchronized manner to multiple signaling pathways to form the circuitry connecting insulin effector sites. Thus, activation of insulin-regulated biological pathways follows a complex time-dependent pattern that requires consideration when evaluating the efficacy of insulin signaling and phosphorylation of downstream effectors. It could be speculated that insulin resistance may have a different impact on targets regulated at the very early or late stages of insulin signaling, respectively, as each time point displays a different pattern of phosphorylated targets.

In skeletal muscle, insulin regulates both transcription and translation of proteins where the former has been attributed to FoxO transcription factors, targets of AKT, that control expression of mitochondrial genes and may contribute to the decreased muscle strength and mass in diabetes³⁵. In muscle cells, insulin has also been reported to affect alternative splicing of the insulin receptor and PKCβ, and growing evidence indicates that alternative mRNA splicing can regulate energy metabolism and may affect cellular insulin sensitivity and susceptibility for diabetes^36,37,38.

Among the most enriched pathways in our study with differential phosphorylation in response to insulin were sets of proteins involved in mRNA processing and splicing. The spliceosome, a key component of the splicing machinery, constitutes a very large, macromolecular RNA–protein complex that assembles de novo on nascent pre-mRNA molecules for subsequent removal of intronic sequences¹⁹. We identified several spliceosomal proteins with distinct temporal patterns of differential phosphorylation upon insulin stimulation. The abundance of the respective phosphopeptides exhibited site-specific kinetics, indicating that insulin stimulation results in increases or decreases of phosphorylation of targets at different time points. The majority of these phosphopeptides were derived from proteins involved in the early steps of spliceosome formation, including assembly of the pre-catalytic spliceosome complexes A and B that precedes the removal of introns and the ligation of neighboring exons to produce mature mRNA for protein biosynthesis on the ribosome.

To further investigate the link between insulin and alternative splicing, we conducted an RNASeq analysis of human primary myotubes with and without insulin stimulation. Consistent with a previous transcriptional analysis of skeletal muscle in mice, short-term (<1 h) stimulation with insulin has only a relatively minor effect on differential gene expression³⁹. Early responding genes were associated with cellular metabolism and transcriptional control, including genes (e.g. NR4A1, EGR1) previously implicated in regulating insulin sensitivity in skeletal muscle cells^40,41. In addition, insulin stimulation was associated with the altered abundance of alternatively spliced mRNAs, including transcripts with retained introns and NMD variants that are predicted to reduce the half-time of the respective mRNAs⁴².

SRSF protein kinase 2 (SRPK2), a target for mTORC1-activated S6K1 (S494) and subsequently CK1 (S497) has been demonstrated to undergo nuclear translocation in its phospho-form to phosphorylate splicing factors such as members of the SR family proteins that promote exon-inclusion or exon-skipping during alternative splicing⁴³. Interestingly, in HEK293E cells, the knockdown of SRPK2 led to increased intron retention of lipogenic genes, triggering nonsense-mediated mRNA decay (NMD), whereas in MCF-7 breast cancer cells, stimulation with IGF1 for 6 h had the opposite effect, reducing intron retention and increasing the stability of transcripts^44,45. In line with these data, we found that short-term simulation of primary human myotubes with insulin was associated with phosphorylation and activation of SRPK2, phosphorylation of SR proteins (e.g. SRSF1 and SRSF2) and reduced intron retention of transcripts like EPS8, a signaling protein involved in regulating actin dynamics⁴⁶ or exon skipping and NMD of UGCG, a ceramide glucosyltransferase and regulator of glucose metabolism⁴⁷.

A recent study of iPSC-derived human myotubes from donors with T2D reported diabetes-associated alterations in the phosphorylation of proteins involved in mRNA splicing⁹. While several components of the spliceosome and pre-mRNA binding proteins, SRRM2, SF3B2, SRSF5, SRSF7, and HNRNPM displayed reduced phosphorylation in cells from T2D donors, phosphorylation of other splicing-associated proteins including SNRNP70, SF1, and HNRMPAB were increased compared to controls⁹. However, in this study, only a few phosphosites in these spliceosomal proteins were reported to be regulated by insulin, SF3B1 (S190), and SNRNP70 (S410), whereas our analysis revealed many additional insulin-regulated phosphosites in the spliceosome over the entire time course. These differences may relate to the different cell types, i.e. primary cells with native epigenetic marks vs. stem cells, quantitation of different phosphopeptides, and/or differences in temporal resolution, the latter being highly likely as evidenced by the intermittent phosphorylation patterns of many sites. Nevertheless, our finding that phosphorylation of spliceosomal proteins, which has been shown to be reduced in T2D, was, in fact, increased in response to insulin, albeit at different target sites, in primary myotubes supports the hypothesis that insulin resistance may generally impair the functionality of the splicing machinery. Thus, while our data indicate that in human primary myotubes, insulin acutely regulates transcription and alternative splicing through phosphorylation of spliceosomal proteins, further work is necessary to investigate the scale and scope of insulin’s acute effects.

Our study has several limitations, including (i) only partial overlap of the phosphoproteome and global proteome, preventing the calculation of absolute phosphosite occupancies, (ii) limited sample size due to the exploratory design of the study which needs follow-up in a cross-sectional clinical study, covering a broader spectrum of phenotypes, (iii) rudimentary assignment of substrates with kinases and phosphatases, and (iv) insufficient information on the functional relevance of insulin-regulated phosphosites and splicing events, which require further studies. In conclusion, our analysis revealed time-dependent phasic phosphorylation patterns in insulin-induced signal propagation, demonstrating that regulation of specific kinases, phosphatases, and corresponding downstream pathways are timely synchronized through cellular feedback loops. Our data highlight the importance of analyzing the kinetics in insulin signaling and provide a resource for assessing donor-specific and temporal variation, which might aid the identification of critical nodes in the circuitry of insulin signaling. Finally, our study links insulin-regulated phosphorylations of the splicing machinery with alternative splicing of mRNA, suggesting hormonal control of the structure and function of the transcriptome.

Methods

Skeletal muscle satellite cell isolation

All study participants provided informed written consent regarding the sampling of biopsies and publication of data. The study design and conduct complied with all relevant regulations regarding the use of human study participants and was conducted in accordance with the criteria set by the Declaration of Helsinki. Ethical approval was granted by the Regional Committee for Medical and Health Research Ethics South East, Oslo, Norway (reference number: 2011/2207). The donors were healthy men, 24 ( ± 1.04) years old, with a body mass index of 24.57 ( ± 0.74) kg/m² (Supplementary Table 1). Biopsies of the vastus lateralis muscle were obtained under local anesthesia (lidocaine) with a Bergström-needle. The biopsy material was released from fat tissue and kept in Ham’s F10 medium (Thermo Fisher Scientific) at 4 °C. The tissue was minced, and the satellite cells were isolated by three subsequent trypsin (0.125%) digestions at room temperature. The supernatants from this digestion were combined and spun down at 500 × g for 7 min to collect the cells. The cells were resuspended in PromoCell skeletal muscle cell growth medium (PromoCell GmbH) supplemented with PromoCell SupplementMix (PromoCell GmbH), 25 IU penicillin (Thermo Fisher Scientific), 25 µg/ml streptomycin (Thermo Fisher Scientific) and 1.25 µg/ml amphotericin B (Thermo Fisher Scientific). The resuspended cells were transferred to a collagen-coated (0.01%) 25 cm² cell culture flask. Cells were kept at 37 °C, 5% CO₂, and medium was changed every 24–48 h until confluent.

Culture of human myotubes

The cells were cultured on cell culture dishes (145 × 20 mm) in DMEM-Glutamax^TM (5.5 mM glucose) (Thermo Fisher Scientific) supplemented with 10% FBS (Gibco), 25 IU penicillin, 25 µg/ml streptomycin, 1.25 µg/ml amphotericin B, 50 ng/ml gentamicin (Sigma-Aldrich), 0.05% BSA, 10 ng/ml hEGF (Thermo Fisher Scientific), 0.39 µg/ml dexamethasone (Sigma-Aldrich) and 25 mM HEPES (Sigma-Aldrich). When the cells had grown to approximately 80% confluence, the growth medium was replaced by differentiation medium (DMEM-Glutamax^TM (5.5 mM glucose) supplemented with 10% horse serum (ATCC), 25 IU penicillin and 25 µg/ml streptomycin). The cells were incubated at 37 °C in a humidified 5% CO₂ atmosphere, and the medium was changed every 2–3 days. All experiments were performed after 6 days Insulin-stimulated glycogen synthesis was performed as previously described¹¹.

Insulin stimulation for phosphoproteomics analysis

On day 6 of differentiation, the medium was aspirated, and the cells were washed twice with PBS before given starvation medium (DMEM without phenol red (5.5 mM glucose) supplemented with 4 mM Glutamax, 25 IU penicillin and 25 µg/ml streptomycin). The cells were serum starved for six hours prior to stimulation with or without 100 nM insulin for 0, 1, 2.5, 5, 15, 30, and 60 min (t₀, t₁, t_2.5, …, t₆₀) before being lysed and harvested in a buffer containing Tris, DTT, protease and phosphatase inhibitor, and snap-frozen in liquid N₂. Protein concentration was determined using the Pierce^TM BCA assay kit.

Sample preparation and mass spectrometry (MS)

For phosphoproteomic profiling, differentiated skeletal muscle cells were solubilized in denaturing SDS buffer (100 mM Tris-HCl, 4% SDS, and 100 mM DTT, supplemented with Complete^TM (Roche) and PhosSTOP (Merck)) and homogenized by passing 10 times through a syringe/26 gauge needle, followed by sonication. After centrifugation at 75,000 × g for 30 min at 4 °C, supernatants were transferred to fresh reaction tubes and protein concentrations were determined by direct photometric measurements (Nanodrop, Thermo Fisher Scientific). Protease digestion was conducted with the filter-aided sample preparation (FASP) method⁴⁸. Briefly, 1 mg of total protein lysate was incubated for 10 min at 96 °C and subsequently diluted 1:10 (v/v) with urea buffer (UB: 8 M urea, 100 mM Tris-HCl pH 8). After concentration using filter devices (Amicon 30 kDa, Merck), samples were alkylated (50 mM iodoacetamide) for 15 min at RT. After three times washing with UB, protein digestion was performed utilizing LysC/Trypsin mix (Promega) in a 1:25 (w/w) enzyme/protein ratio overnight at 37 °C.

Peptides were collected by centrifugation through the filter device. Subsequently, peptides were acidified with trifluoroacetic acid (TFA), a final concentration of 0.1% (v/v), and purified by C18 solid phase extraction (Strata C18-E, 1 g/6 ml, Phenomenex) according to the manufacturer’s instructions. Eluates were lyophilized and subjected to Sequential enrichment from Metal Oxide Affinity Chromatography (SMOAC)⁴⁹. For consecutive enrichment of phosphopeptides, we applied the peptides to HiSelect TiO₂ phosphopeptide enrichment columns (Thermo Fisher Scientific PN#A32993) and used the corresponding flow-through and wash fractions to load HiSelect Fe-NTA phosphopeptide enrichment columns (Thermo Fisher Scientific PN#32992). All steps of SMOAC enrichment were carried out according to the manufacturer’s instructions.

For MS analysis, lyophilized phosphopeptides were then reconstituted in 1% TFA (v/v), including iRT peptides (indexed retention time peptides, Biognosys) and separated by liquid chromatography (Ultimate 3000, Thermo Fisher Scientific). Accordingly, peptides were trapped on an Acclaim PepMap C18-LC-column (ID: 75 μm, 2 cm length; PN#164535 Thermo Fisher Scientific) and subsequently separated on a 50 cm µPAC^TM (PN#552503518050, Pharmafluidics) column connected to an EASYspray ion source of an Orbitrap Fusion Lumos mass spectrometer (Thermo Fisher Scientific).

Peptides were separated using a 100 min linear gradient from buffer A (0.1% formic acid) to 1–34% buffer B (80% ACN, 0.1% formic acid) at a flow rate of 300 nl/min followed by a 20 min linear gradient increasing buffer B to 50% and a 1 min linear gradient increasing buffer B to 90%. Column temperature was set to 40 °C.

MS data were acquired in duplicate on an Orbitrap Fusion Lumos instrument (Thermo Fisher Scientific) operated in data-dependent mode. MS spectra were obtained at 120,000 resolution (3 s cycle time), m/z range of 350–1600, and a target value of 4⁵ ions, with a maximum injection time of 50 ms. For fragmentation precursor selection, filters were set to charge state between 2 and 7, dynamic exclusion of 30 s, and intensity threshold of 2.5⁴. Fragmentation of precursors was done with an isolation window (m/z) 1.2, HCD (Higher-energy collisional dissociation) energy of 32%, Orbitrap resolution of 30,000, and an ion target value of 1.0⁵ with a maximum injection time of 120 ms.

Analysis of mass spectrometry raw data

Mass spectrometry raw files were analyzed with Proteome Discoverer^TM 3.0 software (Thermo Fisher Scientific). ‘SpectrumRC’ node was used with the FASTA database (UniProtKB database, reviewed SwissProt, Homo Sapiens TaxID=9606, v2022-01-30 for phosphoproteome and v2022-03-30 for total proteome analysis) to recalibrate spectra. For quantification purposes, the ‘Minora Feature Detector’ node was used with standard settings (minimum trace length 5, max. delta RT of isotope pattern multiplets of 0.2 min, and for feature to ID linking, use only high confident PSMs). HTsequest search was done against UniProtKB database (reviewed SwissProt, Homo Sapiens TaxID=9606, v2022-01-30 for phosphoproteome and v2022-03-30 for total proteome analysis), the enzyme was set to trypsin with maximum of 2 missed cleavage sites allowed. For HCD fragmentation, b and y Ions were selected with a fragment mass tolerance of 0.02 Da. Carbamidomethylation of cysteine was set as a fixed modification. N-terminal acetylation, N-terminal methionine loss, N-terminal methionine loss, and acetylation, methionine oxidation, as well as phosphorylation (+79.699 Da) of serine, threonine, or tyrosine were allowed as variable modifications. The percolator node (max delta Cn: 0.01) and the ptmRS node included in Proteome Discoverer^TM were applied. Label-free quantification of phosphopeptides using nested design was performed with RT alignment of max 10 min, precursor mass tolerance of 10 ppm, and a minimum S/N threshold of 5 for linked and mapped features. Precursor abundance quantification was based on intensity and a minimum occurrence of 20% in the replicate features. Sample loading normalization of phosphopeptide abundances was performed by the ‘Precursor Ions Quantifier’ node based on total phosphopeptide content. Missing values were imputed using replicate-based resampling.

Sample preparation for RNA sequencing

Human myotubes were cultured as described above. RNA from basal and insulin-stimulated cells (100 nM, 1 h) was isolated using the RNeasy® Plus Mini Kit (Qiagen PN#74134) according to the Manufacturer’s instructions. Whole transcriptome sequencing (2 technical replicates) was conducted at the Genomics & Transcriptomics Laboratory (https://www.gtl.hhu.de/en/). Prior to library preparation, total RNA samples were fluorometrically quantified via Qubit RNA HS Assay (Thermo Fisher Scientific), and quality was measured by capillary electrophoresis using the Fragment Analyzer and the ‘Total RNA Standard Sensitivity Assay’ (Agilent Technologies). All samples in this study showed high-quality RNA Quality Numbers (RQN above 9.8). The library preparation was performed according to the manufacturer’s protocol using the ‘Illumina Stranded Total RNA Prep, Ligation with Ribo-Zero Plus Library Prep Kit’ (Illumina Inc.; Document # 1000000124514 v02). Briefly, 500 ng total RNA was used for ribosomal depletion, fragmentation, synthesis of cDNA, adapter ligation, and library amplification. Bead-purified libraries were normalized and subsequently sequenced on the NextSeq 2000 system (Illumina Inc.) with a paired-end read setup of 2 × 151 bp and a sequencing depth of around 30 million reads per sample. The bcl2fastq2 tool was used to convert the bcl files to fastq files as well for adapter trimming and demultiplexing. The median/mean CVs of the read counts (26,848 transcripts) from the 3 donors were approx. 27%/33% for both basal and insulin-stimulated states, indicating moderate variation of transcript abundance in the samples.

Bioinformatics

General

All bioinformatic analyses, unless explicitly described otherwise below, were performed using a custom data analysis pipeline that we implemented in R v4.2.1⁵⁰. For data analysis using this pipeline, quantitative phosphoproteomics data were imported from ProteomeDiscoverer into R. The normalized abundances of peptides, time point-specific abundance ratios of peptides against basal time point t₀, abundance ratio p-values, adjusted abundance ratio p-values, peptide and protein identifications, and phosphosite localizations generated by Proteome Discoverer were adopted for our bioinformatic analyses.

Regulated phosphopeptides

As a basic analysis step, we identified per time point those phosphopeptides that showed significant differential phosphorylation compared with the basal time point t₀. According to our criteria, these differential or regulated phosphopeptides had to have an abundance ratio against t₀ of >1.5 (upregulated phosphopeptides) or <2/3 (downregulated phosphopeptides), as well as an abundance ratio p-value of <0.05. To retain a sufficiently large number of differential phosphopeptides in all time points over time and for all donors as input for the downstream cross-time point analyses, the unadjusted p-value was chosen for the threshold at this point. The 2741 differential phosphopeptides identified using these criteria were visualized with volcano plots, where the dashed lines represent the thresholds p < 0.05 and fold change >1.5 or <2/3 vs. t₀. In contrast, if the adjusted p-value had been chosen here, only 1521 differential phosphopeptides would have been included in the downstream analyses.

Principal component analysis and total proteome

For principal component analysis (PCA) and all subsequent time course overlapping analyses, a filtered dataset was used with 11,517 phosphopeptides quantified at all time points t₁–t₆₀ in all subjects (i.e., no missing (‘NA’) values; for comparison: 11,612 phosphopeptides were counted as quantified if they were quantified at all time points in at least one donor only). Moreover, for PCA computation and visualization based on regulated phosphopeptides, the most differential 714 phosphopeptides were filtered using the more stringent criteria of p-value < 0.01 and ratio against t₀ > 3 or <1/3 (for comparison: 809 phosphopeptides were quantified at all time points in at least one donor only) to obtain a robust group separation. For comparison, PCA with all 11,517 phosphopeptides without NA values has also been performed.

To account for potential differences in protein abundance, the 714 most differential phosphopeptides of the phosphoproteome were mapped to the 167 corresponding non-phosphorylated peptides quantified in the total proteome with the same amino acid sequence. Conversely, all 25,803 non-phosphorylated peptides quantified in the total proteome were mapped to the 125 corresponding phosphopeptides quantified in the analysis of the phosphoproteome and included in the set of the 714 most differential phosphopeptides. The respective mapped peptides in the phosphoproteome and total proteome were used to perform PCAs in the phosphoproteome and total proteome, respectively. For comparison, also PCA with all 25,803 total proteome peptides was performed.

Analysis of phosphorylation kinetics

To investigate the inter-donor variability in phosphorylation kinetics of phosphosites, only the abundance values of the 8301 monophosphorylated phosphopeptides quantified at all time points t₀–t₆₀ and in all donors were considered. Pearson’s correlation coefficients (Pearson’s r) were calculated for each of these phosphopeptides for all 10 different donor pairs. In the method used, the seven time point-specific abundance values of a particular donor were combined into a temporal profile for each phosphopeptide. Mathematically, the profile ${a}_{i}$ for the Donor $A$ and phosphopeptide $i$ is a seven-dimensional vector of abundance values ${a}_{{ij}}$ where $j\in \{{t}_{0},{t}_{1},{t}_{2.5},{t}_{5},{t}_{15},{t}_{30},{t}_{60}\}$ represents one of the seven time points. Each of the five donor-specific temporal profiles ${a}_{i}$, ${b}_{i}$, ${c}_{i}$, ${d}_{i}$ and ${e}_{i}$ represents the realization of a variable that is independent of the other four variables, since sample-matched values, by definition, are only found within a profile and not between two profiles. To assess donor variability, we calculated Pearson’s correlation coefficient $r$ for each of the 10 pairwise combinations $({a}_{i}\,{,{b}}_{i})$, $({a}_{i}\,,\,{c}_{i})$, $({a}_{i}\,,\,{d}_{i})$, $({a}_{i}\,,\,{e}_{i})$, $({b}_{i}\,,\,{c}_{i})$, $({b}_{i}\,,\,{d}_{i})$, $({b}_{i}\,,\,{e}_{i})$, $({c}_{i}\,,\,{d}_{i})$, $({c}_{i}\,,\,{e}_{i})$, $({d}_{i}\,,\,{e}_{i})$ of the five independent profiles, as well as the average correlation coefficient across the 10 pair combinations. Pearson’s correlation coefficient, which identifies linear correlation relationships, allowed us to validate kinetics with very low inter-donor variability using 8301 plots of donor-specific curves of phosphopeptide abundance values in the log2 scale. We used the average of the respective Pearson’s r of the 10 donor pairings as a measure of the inter-donor variability of the phosphorylation kinetics of the respective phosphosites.

Kinase-substrate overrepresentation analysis

To analyze kinase activity over time after insulin stimulation, we performed a kinase-substrate overrepresentation analysis (KSORA). In this analysis, an overrepresentation analysis (ORA) was used to test whether kinases are statistically significantly overrepresented at a given time point t₁–t₆₀ because their respective phosphosites are located in phosphopeptides that are differential at that time point. We downloaded the necessary information on which phosphosite is a substrate of which kinase from the PhosphoSitePlus database¹⁵ (database version as of 09/09/2022). This included 403 unique human kinases with their annotated substrate phosphosites (13,664 kinase phosphosite annotations, i.e. an average of 33.9 phosphosites per kinase and 9540 unique phosphosites). For KSORA, we used a total of 319 phosphosites quantified in our dataset (out of a total of 2996 phosphosites of all 4030 regulated phosphopeptides, which did not need to be quantified at all time points) for which we found a matching kinase-substrate annotation in PhosphoSitePlus. This made it possible to check which of these phosphosites were located in a differential phosphopeptide at which time point t₁–t₆₀. With this information, a Fisher’s exact test-based p-value was calculated for each of the corresponding kinases using the function enrichmentTest of the R package ClueR v1.4⁵¹, expressing their overrepresentation. These p-values were adjusted using the Benjamini–Hochberg method to obtain associated FDR values. Finally, KSORA results were visualized in a heat map of 19 kinases with the lowest p-values (p < 0.075), in which the colors represent the associated −log10(FDR) values in the z-score scale of each kinase across the six time points. To further clarify the KSORA results, the respective numbers of substrates over time were visualized with bar plots for the 10 most overrepresented kinases. Finally, the top 35 kinases with the lowest p-values (p < 0.15) were placed in a phylogenetic kinase family tree. The software tool Coral was used for this purpose⁵².

Network propagation

The PPI network provided by the ConsensusPathDB, version 35⁵³ was used as a scaffold for mapping the phosphoprotein data. This network has been integrated from 18 different public resources such as BIND, Biogrid, Reactome, and INTACT, among others, and showed very good performance in a comprehensive network propagation benchmark study⁵⁴. Individual PPIs have been assigned a confidence value, c, (range [0,1]) that incorporates topological as well as annotation-based criteria⁵⁵. For this study, we used the high confidence PPIs (c ≥ 0.9) resulting in a network of 10,642 different proteins and 139,266 interactions.

We used a subset of 1552 significantly regulated phosphopeptides with continuous quantitative information for each time point (no missing values in t₁–t₆₀; however, missing values in t₀ allowed; adjusted p-value < 0.05 in at least one of t₁–t₆₀), which were mapped to 905 unique HNGC accessions. For each time point t₁–t₆₀, a separate network propagation was performed: proteins (network nodes) were scored according to the significance of their deviation from t₀. Let p_ij be the adjusted p-value of the statistical test for an abundance of phosphopeptide i in time point t_j compared with time point t₀, then the score for the corresponding protein is -log(p_ij). If more than one phosphopeptide is associated with the protein, the largest score is used. Thus, for each time point, the network initialization reflects the significance of the phosphopeptide abundance differences from the control time point.

Network propagation was done with the NetCore tool¹⁸ with a two-step procedure: In the first step, initial node weights were spread through the network with a modified random walk with a restart procedure to obtain a re-ranking of the initial node weights. The reweighted nodes were then compared against random permutations of the network, and a p-value was assigned to each node. In the second step, significantly reweighted nodes were connected with seed nodes (top 100 nodes from the initialization) in order to compute subgraphs that represent parts of the PPI network that are mostly affected by the data.

Network dynamics analysis

Following the network propagation for each time point, the resulting p proteins and their interactions were merged into a single non-redundant network, which was subsequently integrated as nodes and edges into Cytoscape V3.9.1⁵⁶. Additional databases and manually curated information regarding protein function (Kinase or Phosphatase) as well as derivation of the nodes (through input or propagation) were added and mapped to the network as color and shape, respectively.

Pathway overrepresentation analysis

For the pathway overrepresentation analysis pathways from the Reactome pathway database⁵⁷ were used. Our data analysis pipeline used the R package ReactomePA v1.4.0⁵⁸ to perform the pathway overrepresentation analysis based on differential proteins found in early (t₁, t_2.5), intermediate (t₅, t₁₅), and late (t₃₀, t₆₀) time points. Thus, a total of 734 unique differential phosphopeptides from the early time points, 816 from the middle time points, and 1080 from the late time points were mapped to the respective genes (Entrez IDs) of their respective unique proteins. The resulting three sets of genes were passed as input to the enrichPathway function of ReactomePA using Fisher’s exact test to calculate p-values and adjust them (FDR computation via the method of Benjamini and Hochberg) in order to identify overrepresented Reactome pathways. Then, 60 statistically significantly overrepresented pathways with the lowest p-values (where all respective FDRs <0.00175 after p-value adjustment) were taken from each of the three time periods early, intermediate, and late. Due to overlaps, this was a total of 90 distinct overrepresented pathways. Of these, 30 pathways covering the most diverse biological functions were selected and visualized in a heat map (function heatmap.2 of the R package gplots v3.1.3⁵⁹), where colors represent z-scores of −log10(FDR) values, to show the varying overrepresentation of their respective phosphoproteins in the early, intermediate and late time points after insulin stimulation. For six selected Reactome pathways, the respective numbers of associated up- and downregulated differential phosphopeptides were visualized with bar plots across the three time periods using the R packages ggplot2 v3.3.6⁶⁰ and ggpubr v0.4.0⁶¹. Moreover, for each of the 30 pathways shown in the heat map, the mean values of their respective average inter-donor Pearson’s r of the respective associated differential phosphopeptides were computed as a measure for inter-donor pathway variability. After sorting the 30 pathways by these mean values, the pathway-specific distributions of average inter-donor Pearson’s r were visualized along with their respective mean values in a strip chart.

Splicing analysis

The general pipeline was generated with Python (v.3.7) and R (v4.2.3) in snakemake⁶². RNA reads were mapped to the human genome (hg38, GRCh38.primary_assembly.genome.fa) with STAR v2.7.10b⁶³ using the GENCODE annotation (gencode.v44.primary_assembly.annotation.gtf)⁶⁴. Subsequently, the genome alignments were analyzed with rMATS-turbo (v4.2.0)⁶⁵ and additionally counted with kallisto⁶⁶. Further downstream analysis was done following the pipeline and parameters according to Love⁶⁷ with DRIMSeq⁶⁸.

Temporal dynamics of relative phosphosite occupancy

We used the abundance value ratios of the phosphopeptides relative to time zero (t₀) to determine the relative occupancy of the respective site. The information on the number and identity of phosphopeptides that are differentially regulated compared to t₀ for each time point (t₁–t₆₀) was plotted as a Sankey flow diagram to visualize the temporal dynamics of the underlying phosphorylations and dephosphorylations. The R packages ggplot2 v3.3.6⁶⁰, ggsankey v0.0.99999⁶⁹, and dplyr v1.0.9⁷⁰ were used for this purpose. We used a self-developed R code to group the regulated phosphopeptides into temporal strands, allowing assessment of kinetics and duration of regulated phosphorylation events, as well as visualization of strand-specific properties and distributions of relative abundances (each normalized by the maximum abundance) of the respective up- or downregulated differential phosphopeptides using box plots.

Strand- and time-point-specific sequence motifs of the phosphorylation sites were visualized with the function ggseqlogo of the R package ggseqlogo v0.1⁷¹. For the visualization of the sequence motifs using ggseqlogo, the log10(bits+1) height method was used, with the calculation of bits and their log10 scaling performed with self-developed R code. The data structure was also used to automatically perform time point-specific pathway ORAs using the R packages ReactomePA v1.40.0⁵⁸ and topGO v2.48.0⁷² for the respective temporal strands and compare manually selected ORA results using multiple databases with the CPDB web application⁵³. Based on the ORA results, the mapping to temporal strands was manually visualized for some selected Reactome and GO pathways in diagrams depicting differential phosphopeptides of the respective pathway and their respective phosphosites as color-coded containers. In these diagrams, the presence of a container represents the quantification of the corresponding differential phosphopeptide at that particular time point, and its color code corresponds to the strand-specific color code in the Sankey flow diagram.

Authors’ relationships and activities

The authors declare no relationships or activities that might bias, or be perceived to bias, their work.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.