Metabolic disease is caused by a combination of genetic and environmental factors, yet few studies have examined how these factors influence signal transduction, a key mediator of metabolism. Using mass spectrometry-based phosphoproteomics, we quantified 23,126 phosphosites in skeletal muscle of five genetically distinct mouse strains in two dietary environments, with and without acute in vivo insulin stimulation. Almost half of the insulin-regulated phosphoproteome was modified by genetic background on an ordinary diet, and high-fat high-sugar feeding affected insulin signalling in a strain-dependent manner. Our data revealed coregulated subnetworks within the insulin signalling pathway, expanding our understanding of the pathway’s organisation. Furthermore, associating diverse signalling responses with insulin-stimulated glucose uptake uncovered regulators of muscle insulin responsiveness, including the regulatory phosphosite S469 on Pfkfb2, a key activator of glycolysis. Finally, we confirmed the role of glycolysis in modulating insulin action in insulin resistance. Our results underscore the significance of genetics in shaping global signalling responses and their adaptability to environmental changes, emphasizing the utility of studying biological diversity with phosphoproteomics to discover key regulatory mechanisms of complex traits.
This fundamental study provides a unique tool for assessing the range of phosphorylation in insulin reactions due to genetic variation and dietary influence through the utilization of genetically distinct mouse strains. The discoveries of this study hold substantial importance, as they shed light on the interplay between genetic attributes and environmental conditions in shaping the insulin-signaling network within skeletal muscle, a crucial regulator of metabolism. The supporting evidence presented is compelling, and the work is anticipated to captivate a wide audience within the metabolism discipline due to its extensive appeal and by providing inspiration for further hypothesis-driven research.
Protein post-translational modifications such as phosphorylation enable cells to rapidly respond to environmental changes by modifying protein function at low metabolic cost1. As a result of this high metabolic efficiency, phosphorylation is involved in nearly all biological processes and is dysregulated in numerous complex diseases2. Advances in mass spectrometry-based phosphoproteomics – the unbiased identification and quantification of protein phosphorylation – have led to the discovery of more than 100,000 phosphosites, revealing that the phosphoproteome comprises vast, interconnected phosphorylation networks3–6, rather than the textbook view of isolated, linear kinase cascades.
The insulin signalling network is among the most studied phosphorylation networks. Post-prandial increases in blood glucose stimulate pancreatic insulin secretion, coordinating a metabolic switch in target tissues like skeletal muscle and adipose7. Insulin increases glucose uptake into these tissues by promoting translocation of the glucose transporter GLUT4 to the plasma membrane, and serves other functions like enhancing protein synthesis, downregulating lipid catabolism, and altering gene transcription7. To coordinate these functions, insulin triggers a phosphorylation cascade primarily involving activation of the Ser/Thr kinase Akt, regulation of downstream kinases including mTORC1 and GSK3, and modulation of parallel signalling arms7–9. Seminal phosphoproteomics studies demonstrated that this cascade regulates over a thousand phosphosites, with many still uncharacterised in insulin action10–12. Insulin resistance – the failure of insulin to promote glucose uptake in its target tissues – is triggered by genetic and environmental factors such as family history of metabolic disease and high-calorie diets13. Although insulin resistance is a major precursor of metabolic disease including type 2 diabetes, its mechanistic basis remains unresolved13, 14.
Interactions between genetics and environment significantly regulate biomolecular processes, including insulin resistance15–18. As signalling pathways connect the extracellular environment to intracellular proteins, they are likely a major conduit of gene-by-environment interactions. Yet, how global phosphorylation signalling networks are regulated across different genetic backgrounds is relatively unexplored. Recent phosphoproteomics studies in yeast19 and mice20 identified genetic variants affecting multiple phosphosites, but did not analyse the phosphoproteome’s response to acute perturbation, which is crucial to its role as a signal transduction system. We have also shown marked variation in acute signalling responses to exercise or insulin across individuals21. However, this study did not systematically assess the relative contributions of genetics and the environment21. Understanding how these variables intersect with signal transduction is fundamental to our basic knowledge of signalling and the advancement of personalised medicine, which advocates individualised treatment regimens based on genetic risk factors and gene-by-environment interactions15, 22.
Inbred mice allow precise control of genetics and environment unachievable in human studies, enabling examination of how these factors interact to influence biomolecular systems23. Here, we performed phosphoproteomics on insulin-stimulated or control skeletal muscle from five genetically distinct inbred mouse strains fed either an ordinary chow diet or a high-fat, high-sugar “western style” diet. Strikingly, we found that genetic background influenced both the phosphoproteomic insulin response of chow-fed mice, as well as how these responses were modified by high-fat, high-sugar feeding. These signalling changes were reflected in altered activity profiles of multiple kinases and provided insight into the functional organisation of the insulin signalling network by revealing subnetworks of coregulated phosphosites. A major challenge in phosphoproteomics studies is pinpointing important regulatory events among the many responding to a stimulus. We reasoned that associating changes in protein phosphorylation across the gene-by-environment landscape with phenotypic change – in this case insulin-stimulated glucose uptake – would dissect mechanistic targets with greater fidelity. This approach generated known as well as candidate regulators of insulin-stimulated glucose uptake, leading us to demonstrate that glycolytic upregulation reverses insulin resistance. Our work represents the first global portrait of insulin signalling network plasticity in response to genetic and environmental variation, which will serve as an important resource in future studies of insulin action and resistance.
Phosphoproteomics of insulin signalling in mouse skeletal muscle
To study how protein phosphorylation networks are affected by genetics and environment, we examined insulin signalling in five genetically distinct inbred mouse strains including four lab strains with diverse metabolic phenotypes (C57Bl6J, NOD, BXH9, and BXD3417), and the wild-derived CAST strain (Fig. 1a). Mice underwent a six-week diet regimen of standard lab diet (CHOW) or a high-fat high-sucrose diet (HFD), which is commonly used to induce insulin resistance17, 24. Consistent with their diverse genetics, these strains differed in morphometric parameters (body weight, adiposity) and metabolic traits (fasting blood glucose, fasting blood insulin, glucose tolerance) both on the CHOW diet and in their response to HFD-feeding (Fig. S1).
We focused on skeletal muscle, as it is the site of greatest post-prandial glucose uptake and the most significant contributor to impaired glucose disposal in type 2 diabetes25. Specifically, we chose to examine the soleus muscle, because its largely oxidative fibre composition resembles human muscle tissue more than other murine muscles26. Mice were injected retro-orbitally with saline control or insulin for 10 minutes, and the soleus was collected for phosphoproteomic analysis (Fig. 1a). A tritiated 2-deoxyglucose (3H-2DG) tracer was co-injected to measure soleus glucose uptake.
Using the EasyPhos workflow and data-independent acquisition (DIA) mass spectrometry27, 28, we quantified 28,809 phosphopeptides across 95 biological samples, corresponding to 23,126 unique high-confidence phosphosites (Class I; localization score > 0.75) on 3,507 proteins (Fig. 1b, Table S1). On average we quantified 15,395 phosphopeptides in each sample, with minimal variation (Fig. S2a). Samples from animals of the same strain and diet were highly correlated and generally clustered together, implying the data are highly reproducible (Fig. S2b-d).
To explore the soleus insulin signalling network, we examined phosphopeptides altered by insulin stimulation in at least one of the ten strain-diet combinations. We filtered our data to 10,432 phosphopeptides highly quantified across most conditions (see Methods) and identified those with significant differences between unstimulated and insulin-stimulated samples (three-way ANOVA insulin main-effect q-value < 0.05) that was of sufficient magnitude in at least one strain-diet combination (insulin/unstimulated fold change > 1.5). This resulted in 441 insulin-regulated phosphopeptides on 232 proteins, which is noticeably more than recent studies of 10-minute insulin signalling in patient-derived myoblasts (174 phosphopeptides29 and 242 phosphopeptides30) or mouse adipose tissue (319 phosphopeptides31) (Fig. 1b, c, Table S1).
Our analysis recovered many well-studied insulin-regulated phosphosites, including Akt substrates such as T247 on Akt1s1 (PRAS40), S939 on Tsc2, and S9 on Gsk3β (Fig. 1d-f), as well as targets of downstream kinases including the Gsk3β substrate S641 on Gys1 (Fig. 1g). Enrichment of Gene Ontology (GO) biological processes recapitulated canonical insulin signalling axes including “insulin receptor signaling pathway”, “phosphatidylinositol 3-kinase signaling”, “protein kinase B signaling”, and “TORC1 signaling”, and multiple pathways related to glucose metabolism, fatty acid metabolism, autophagy, and protein translation, reflecting known targets of insulin action (Fig. S3a). Furthermore, insulin-regulated phosphosites were enriched for insulin regulation in our previous human skeletal muscle phosphoproteome (fold enrichment = 4.22, p = 9.50×10-24, one-sided Fisher’s exact test, Fig. S3b)21. Despite this, only half of all insulin-regulated phosphopeptides (228/441) were previously annotated as insulin-regulated in the PhosphositePlus database32 (Fig. S3c), highlighting the potential of our data to discover novel aspects of insulin signalling while recapitulating known components. Overall, our phosphoproteomics data provide a comprehensive and high-quality atlas of insulin signalling in mouse skeletal muscle.
Genetics and diet modulate insulin signalling
The influence of genetic and environmental variation on global insulin signalling responses is largely unknown. We therefore developed a pipeline to address this question using our phosphoproteomics data (Fig. 2a). First, we converted the intensity values of each insulin-regulated phosphopeptide to “insulin response” values, by normalising insulin-stimulated data to the unstimulated median of the corresponding Strain-Diet combination. Since protein expression should not change within a 10 min insulin stimulation31, this allowed us to parse out protein abundance differences across strains and diets and focus solely on acute signalling processes. We then assessed the impact of genetics in CHOW-fed mice by identifying phosphopeptides with differing insulin responses in one or more strains compared to C57Bl6J (“Strain effect”). Lastly, we explored the effects of HFD-feeding on signalling as two types of “Diet effects”, either as a “Uniform diet effect” – where HFD-feeding affects each strain similarly – or a “StrainxDiet effect” – where its impact depends on the mouse strain. Analyses of “Strain effects” and “Diet effects” were performed separately, so a phosphopeptide could have both a Strain effect and a StrainxDiet/Uniform Diet effect.
Almost half of the 441 insulin-regulated phosphopeptides displayed a Strain effect (Fig. 2b, Table S1). These included phosphopeptides where C57Bl6J had a stronger insulin response than other strains, such as S15 on the RNA methyltransferase Rnmt (Fig. 2c), and phosphopeptides where C57Bl6J had a weaker insulin response, such as S48 on the vesicle fusion regulator Vamp3 (Fig. 2d). Vamp3 S48 is predicted to be highly functional (functional score = 0.75033), its phosphorylation correlates with glucose uptake in insulin-stimulated and/or exercised human skeletal muscle21, and Vamp3 overexpression rescues GLUT4 translocation in insulin resistance34, suggesting that this site may represent a genetically variable control point of GLUT4 trafficking. In general, insulin responses were weaker in the four remaining strains compared to C57Bl6J, though the extent of this trend was strain-dependent (Fig. 2e). In all, the strain-affected phosphopeptides reveal a unique fingerprint of insulin signalling within each strain (Fig. 2f), highlighting the complex and widespread effects of genetic variation on signalling networks.
We next examined the impact of HFD-feeding in insulin signalling. StrainxDiet effects were more prevalent than Uniform diet effects (110 vs 10 phosphopeptides, Fig. 2g, Table S1), suggesting that the molecular impact of dietary perturbation was strongly modulated by genetic background. StrainxDiet effects impacted known regulatory phosphosites such as the inhibitory site S78 on Map2k4, whose insulin response was attenuated by HFD-feeding only in C57Bl6J and CAST (Fig. 2h). Map2k4 activates p38 and JNK kinases which have been implicated as drivers of HFD-induced insulin resistance8, and based on the behaviour of S78, the orchestration of this detrimental signalling axis in HFD-feeding may depend on genetic background.
HFD-feeding exerted complex effects on signalling, with the balance between suppressed and enhanced insulin responses varying across strains (Fig. 2i). Furthermore, insulin-regulated phosphopeptides were largely altered by HFD-feeding in only a single strain, and when multiple strains were affected, they often changed in opposite directions (Fig. 2j-k). For instance, multiple insulin responses were strengthened in BXH9 but weakened in C57Bl6J or CAST (Fig. 2k), such as S500 on the translation regulator Larp4b (Fig. S3d). Principal component analysis supported the highly divergent impact of dietary perturbation, as HFD-feeding displaced each strain in a distinct direction in principal component space (Fig. 2l). StrainxDiet and Strain effects were driven by a mixture of changes to insulin-stimulated phosphorylation, unstimulated phosphorylation, or both, highlighting the complexity of these signalling alterations (Fig. S3e-f). This analysis demonstrates the pervasive role of genetics in shaping signalling networks, as genetic background profoundly modulated the effect of HFD-feeding on insulin signalling.
Exploring genetic and dietary modulation of the insulin signalling network
To understand the insulin signalling circuitry and functional pathways modulated by genetics and diet, we curated a network of 160 insulin-regulated phosphosites, comprising sites from a knowledge pathway-derived list of canonical insulin signalling proteins21 and substrates of insulin-regulated kinases (Fig. 3, see Methods). Strain and diet affected multiple highly studied signalling sites, including Tsc2 S939 (Diet effect), Gsk3α S21 (Strain and Diet effects), and Tbc1d4 T649 (Strain and Diet effects), while other sites such as Gsk3β S9 and Akt1s1 T247 were unaffected. Interestingly, strain and diet affected both canonical and non-canonical insulin signalling proteins to a similar extent (Fig. S3g-h). Non-canonical phosphosites could shed light on underappreciated outcomes of insulin action either altered or unaffected by genetics and the environment, such as the p70S6K substrate S47 on Dnajc2 (Strain and Diet effects), which drives cellular senescence35, and S315 on Pcyt1a (no Strain or Diet effect), which inhibits phosphatidylcholine biosynthesis36 (Fig. 3).
No functional pathways were overrepresented within strain or diet-affected phosphosites relative to all insulin-regulated sites, implying that genetics and environment modulate most or all outcomes of insulin. For instance, strain and diet affected regulatory phosphosites controlling vesicle trafficking (S348, T575, S595, and T649 on the GLUT4 trafficking regulator Tbc1d4); glucose metabolism (S469 and S486 on Pfkfb2, which promote glycolysis); mitochondrial structure and function (S129 on Mff, which promotes mitochondrial fission); autophagy (S555 on the master autophagy regulator Ulk1); gene transcription (the inhibitory site S108 on the transcription factor Tfeb); and mRNA translation (S236 on ribosomal proteins S6, S422 on Eif4b, and S64 and T69 on Eif4ebp1, which promote translation). Interestingly, Strain and Diet effects overlapped significantly (fold enrichment = 1.50, p = 4.40×10-9, two-sided Fisher’s exact test, Fig. S3i), implying some phosphosites may be more amenable to regulation overall. As a notable exception, all six insulin-regulated phosphosites on Plin1 had Diet effects while only one had a Strain effect (Fig. 3). Plin1 coats and regulates lipid droplets, hence this enrichment of Diet effects may represent a signalling response to increased intramuscular lipids in the HFD condition. Overall, genetics and environment triggered widespread alterations in insulin signalling impinging on diverse cellular pathways.
Genetics and diet rewire insulin-regulated kinase signalling
The extensive signalling changes caused by genetics and diet may result from altered kinase regulation. We tested this hypothesis using a kinase substrate enrichment analysis (KSEA)37 on phosphopeptide insulin responses. KSEA accurately captured the activation of canonical insulin-regulated kinases (Akt, mTOR, p70S6K, and p90RSK) and the deactivation of GSK3, confirming the validity of the approach (Fig. 4a). Focussing on CHOW-fed mice, we identified seven kinases differentially enriched across mouse strains (ANOVA adjusted p-value < 0.05, Fig. 4b). For example, insulin activated SGK and deactivated GSK3 more in C57Bl6J and NOD than in other strains (Fig. 4b). Extending this analysis to all mice, we identified kinases with Uniform diet or StrainxDiet effects (Fig. 4c). Akin to our analysis of individual phosphosites (Fig. 2), StrainxDiet effects were more prevalent than Uniform diet effects (five kinases compared to one), indicating that genetic background strongly influences the impact of HFD-feeding on kinase signalling. These results suggest that the observed phosphosite signalling changes could be partly due to altered insulin regulation of multiple kinases.
Biological variation reveals functional organisation of the insulin signalling network
KSEA predicted changes in overall kinase activity, but we questioned if substrates of the same kinase could be differentially regulated by genetic and environmental variation. As a case study we examined substrates of Akt – a master regulator of insulin signal transduction – to assess the similarity of their insulin responses across strains and diets. Strikingly, we observed a range of both positive and negative correlations (Fig. 5a). For instance, while Tsc2 S939 and Akt1s1 T247 both activate mTORC1, their insulin responses correlated poorly (r = 0.202, p = 0.168, Fig. 5b). Supporting these findings, similar heterogeneity in Akt substrate insulin responses has previously been observed in skeletal muscle from humans with differing metabolic health38. Hierarchical clustering revealed four distinct groups of positively correlated Akt substrates (Fig. 5a), suggesting these groups may coordinate distinct functional outputs of Akt signalling.
We next explored if the genetic and environmental variation in our study could reveal organisational principles of the entire insulin signalling pathway. By performing weighted gene correlation network analysis (WGCNA39, 40), we identified eight subnetworks of coregulated insulin-responsive phosphopeptides (Fig. 5c) varying in size from 16 to 120 phosphopeptides, with 91 assigned to no subnetwork (Fig. 5d, Table S2). Examining the subnetwork “eigenpeptides” – a weighted average of the constituent phosphopeptides39, 40 – revealed that the subnetworks captured distinct effects of genetics and diet on insulin signalling (Fig. 5e). For example, HFD-feeding attenuated the insulin response of subnetwork I in CAST and C57Bl6J strains, while subnetwork II was affected by HFD-feeding only in CAST and NOD (Fig. 5e). This suggests that the subnetworks may be sensitive to distinct cellular information.
Next, we characterised the regulatory and functional nature of these subnetworks. Canonical insulin-regulated kinases such as Akt and mTOR were enriched across multiple subnetworks (Fig. 5f), confirming that genetic and environmental variation can reveal uncoupling of substrates targeted by the same kinase (Fig. 5a-b). Nevertheless, visualising these subnetworks within our curated insulin signalling pathway (Fig. 3) revealed cases where signal flowed through a single subnetwork, such as from Erk2 (Y185) to its target kinase Rsk2 (T365 and S369) and Rsk2 substrates (Gab2 S211 and Nos1 S847) within subnetwork III (Fig. S4). Within multiply phosphorylated proteins, phosphosites either belonged to the same subnetwork (e.g. Plin1) or diverse subnetworks (e.g. Tbc1d4 and the transcription factor Nfatc2), suggesting the latter may serve as hubs integrating diverse cellular information (Fig. S4).
Some subnetworks were enriched in specific cellular compartments (GO ontologies), implying that common localisation may facilitate coregulation of phosphosites (Fig. 5g). Coregulation may partition functional outcomes of insulin action, as certain biological processes were enriched only in select subnetworks (Fig. 5h). These included known insulin targets like “negative regulation of lipid catabolic process” in subnetwork I and “positive regulation of glycogen biosynthetic process” in I, VI, and VIII (Fig. 5h). Lastly, we leveraged our previous phosphoproteomic time course of insulin signalling to interrogate subnetwork dynamics10 and found that phosphopeptide insulin response dynamics varied across subnetworks (Fig. 5i). This reveals distinct temporal regulation as another feature underlying the substructure of the insulin signalling network. Overall, genetic and environmental diversity illuminated the complex coregulation structure of insulin signalling, featuring subnetworks that evade known network circuitry and present unique functional signatures.
Leveraging biological variation to identify drivers of insulin responsiveness
We have so far described the marked influence of genetic background and HFD-feeding on skeletal muscle insulin signalling, evident at the level of individual phosphosites, protein kinases, and co-regulated network modules. We hypothesised that by associating this signalling diversity with a phenotypic output of insulin, such as enhanced glucose uptake, we would filter out mechanistically irrelevant phosphosites and hence prioritise molecular regulators of the phenotype. To test this hypothesis, we measured in vivo glucose uptake with 3H-2DG tracer in the same muscle samples used for phosphoproteomics. Insulin-stimulated glucose uptake differed by more than twofold across strains (two-way ANOVA strain effect p = 4.78×10-7) and was almost uniformly decreased by HFD-feeding (two-way ANOVA diet effect p = 1.83×10-5) (Fig. 6a). This highlights that genetic background and dietary status are key determinants of insulin responsiveness.
To prioritise signalling nodes responsible for differences in insulin responsiveness, we correlated all insulin-regulated phosphopeptides with glucose uptake in insulin-stimulated muscles, resulting in 37 significantly correlated phosphopeptides (r > 0.35 or < -0.35, q-value < 0.1, Fig. 6b). The most significantly correlated phosphopeptide contained T1174 and S1176 on the nitric oxide synthase Nos3. The latter serves as a positive control for our analysis, as this site is known to be phosphorylated in response to insulin to promote in vivo glucose uptake by vasodilation41–43 (Fig. 6c). Other correlated phosphopeptides that could modulate insulin responsiveness include S1751 on Afdn, a phosphosite implicated in insulin action44, and S196 on the Prkag2 subunit of AMPK, a major metabolic signalling hub promoting glucose uptake45 (Fig. 6b). These examples suggest that our analysis prioritised regulators of glucose uptake.
While the above analysis identified phosphosites associated with glucose uptake through their absolute abundance, we hypothesized that for some phosphosites, the magnitude of their response to insulin may be a stronger determinant of insulin action. We found that the insulin response values of 13 phosphopeptides correlated with insulin-stimulated glucose uptake (r > 0.35 or < -0.35, q-value < 0.1, Fig. 6d). These were largely distinct from the 37 phosphopeptides identified in our first analysis, indicating that the two approaches captured complementary information. Several of these phosphopeptides could regulate insulin-stimulated glucose uptake, such as the regulatory site S469 on the enzyme Pfkfb2 which activates glycolysis, a major pathway for glucose consumption46–48 (Fig. 6e), and S177 on Rcsd1, which could affect GLUT4 vesicle transport via actin cytoskeleton remodelling49 (Fig. 6d).
In addition to individual phosphosites, the status of larger signalling network components could also influence insulin responsiveness. Kinase enrichment scores affected by strain or diet did not correlate with glucose uptake (Table S3), suggesting insulin action is not dominated by the net activity of specific kinases. Interestingly, two WGCNA-derived insulin signalling subnetworks correlated with glucose uptake: subnetwork V (r = 0.436, p = 0.00173) and subnetwork I (r = 0.332, p = 0.0197, Fig. 6f). Subnetwork V could modulate glucose uptake through actin cytoskeleton remodelling via Rscd1 S177, through glucose metabolism promotion via Gys1 S641 (Table S2), and by influencing GLUT4 vesicle trafficking due to its enrichment at “cytoplasmic vesicle membranes” (Fig. 5g). It was also enriched in substrates of GSK3, which has been implicated in insulin resistance in skeletal muscle50–52 and adipose tissue31. Subnetwork I, the largest cluster containing 27% of insulin-regulated phosphopeptides, was enriched in multiple kinases and biological processes (Fig. 5f, h), suggesting it may be a central regulatory hub for various outcomes of insulin action including glucose uptake. Overall, compartmentalisation of insulin-responsive phosphosites into subnetworks may enable independent control of insulin’s functional outputs, since only two subnetworks correlated with insulin-stimulated glucose uptake.
Upregulating glycolysis reverses insulin resistance
We next aimed to validate our approach for identifying regulatory mechanisms of insulin-stimulated glucose uptake. S469 on Pfkfb2 correlated highly with glucose uptake follow insulin stimulation (Fig. 6e). Phosphorylation of this site leads to increased production of F2,6BP, a potent glycolytic agonist, suggesting that activating glycolysis may play a key role in muscle insulin responsiveness. This is consistent with our previous findings that glycolytic enzyme abundance was strongly associated with ex vivo insulin-stimulated glucose uptake in muscle from inbred mice17, and that decreasing glycolytic flux caused insulin resistance in vitro53. To further establish glycolysis as a regulator of insulin responsiveness in skeletal muscle, we decided to investigate whether upregulating glycolysis through F2,6BP production can restore insulin-stimulated glucose uptake in insulin resistance. We endeavoured to stimulate glycolytic flux in muscle cells independently of regulatory input in case this is compromised by insulin resistance. While Pfkfb2 requires activating phosphorylation by Akt to produce F2,6BP substantially, its paralog Pfkfb3 has high basal production rates and lacks an Akt motif at the corresponding phosphosites48. We therefore rationalised that overexpressing Pfkfb3 would robustly increase F2,6BP production and enhance glycolysis regardless of insulin stimulation and Akt signalling (Fig. 6g). To avoid systemic effects of Pfkfb3 overexpression we studied cultured L6-GLUT4-HA myotubes, which display robust insulin regulation of GLUT4 trafficking and develop insulin resistance upon palmitate treatment, mimicking lipotoxicity, a trigger of in vivo insulin resistance54.
As anticipated, Pfkfb3 overexpression increased glycolytic capacity in L6-GLUT4-HA myotubes as measured by extracellular acidification rate (Fig. S5a-c). Pfkfb3 overexpression also restored insulin-stimulated glucose uptake to normal levels in palmitate-treated cells (Fig. 6h). This effect was only observed in cells treated with palmitate and insulin, suggesting it specifically modulated insulin action rather than non-specifically increasing glucose uptake through enhanced glucose consumption. These results further establish glycolytic flux as a major determinant of the glucose uptake arm of muscle insulin action and highlight the power of studying phosphoproteomics across the gene-by-environment landscape to identify causal drivers of complex phenotypes. We anticipate that our catalogue of glucose uptake-correlated phosphosites will provide a rich starting point for future investigations into mechanisms of insulin action and resistance.
The environment shapes the flow of information from genotype to phenotype. Many studies have interrogated the role of intermediate molecular layers such as the transcriptome or proteome, however few studies have examined how protein post-translation modifications participate in this information transfer. Here we have approached this problem by leveraging diverse inbred mouse strains and phosphoproteomics to examine the insulin signalling network across a landscape of genetic and dietary variation. Genetic background significantly altered the insulin signalling network both independently and in concert with dietary status, affecting myriad phosphosites and multiple kinases. We exploited this variation in signalling responses in two ways – to study the partitioning of the Akt and insulin signalling pathways into distinct subnetworks of coregulated phosphosites; and to identify potential regulators of insulin responsiveness by associating phosphorylation with insulin-stimulated glucose uptake. Finally, validation studies in L6 myotubes confirmed the major role of accelerated glycolysis as a key regulator of insulin responsiveness.
Genetic and diet-driven signalling changes did not transmit linearly through our current model of the insulin signalling network, illustrating that this model remains incomplete. Notably, substrates of kinases such as Akt clustered into distinct groups based on differing insulin responses. Hence, it is an oversimplification to model signalling pathways as networks of individual kinases since substrates of the same kinase display independent regulation. This could arise from localisation of a kinase to distinct substrate pools55–57; interactors targeting a kinase to different substrates55; substrate phosphorylation by alternate kinases58; the dephosphorylation of specific substrates by phosphatases; kinase post-translational modifications altering substrate specificity; and distinct substrate phosphorylation kinetics10. As our knowledge of the repertoire of kinase substrates continues to deepen59, future research should explore how the above mechanisms contribute to finer regulation of these substrates. Genetic and environmental variation also exposed a coregulation subnetwork structure within the insulin signalling network. The enrichment of subnetworks in distinct biological processes, and the selective association of two subnetworks with glucose uptake, suggests that this coregulation structure may direct independent control of distinct outcomes of insulin action. This exciting possibility necessitates further investigation, including replication in independent cohorts, spatiotemporal characterisation of subnetwork dynamics, and association of additional insulin-regulated phenotypes with subnetwork profiles.
Muscle insulin signalling responses vary across individuals21, 38, and our results suggest that baseline genetic differences and an individual’s environment both alter signalling, with the environment’s influence depending strongly on genetic background. Signalling pathways are popular therapeutic targets due to their importance in human health and the relative ease of pharmaceutical interventions60. Our results advocate for a personalised approach to such therapies, implying that the efficacy of treatments aiming to correct pathological signalling responses will depend on an individual’s genetic background. In cancer, where signalling networks are dysregulated heterogeneously, modelling patient-specific networks has already shown promise for identifying personalised drug targets61, 62. Personalised medicine approaches will also be aided by a comprehensive understanding of how genetics shape signalling networks and potentiate their modulation by the environment. Recent studies have made the first step, revealing that the ground-state phosphoproteome can be altered by mutations affecting network components such as kinases, phosphatases, and phosphoproteins, as well as the molecular milieu the network is exposed to including extracellular signalling ligands19, 20. An important corollary of such genetic factors is that multiple genetic backgrounds should be studied when establishing generalizable signalling responses. Our data indicate that insulin responses in C57BL6J – the most commonly studied mouse strain – are not necessarily generalizable, indicated by phosphosites such as S15 on Rnmt that were insulin regulated almost exclusively in C57Bl6J.
A major challenge in studying signal transduction with omics technologies is that hundreds to thousands of molecular events typically respond to a cellular signal, making it difficult to pinpoint the most crucial regulatory nodes. To tackle this challenge, we previously demonstrated that associating phosphoproteomics with physiological phenotypes across diverse individuals enriches for phosphosites more likely to modulate biological responses21. Here we have elaborated on this approach, revealing that associating phosphorylation with phenotype across a genetic and environmental landscape can identify regulators of specific biological processes, such as insulin-stimulated glucose uptake. Our results recapitulated known glucose uptake regulators and led to further validation of glycolytic flux as a modulator of insulin responsiveness. We have previously demonstrated that reduced glycolytic flux impairs GLUT4 translocation and insulin signalling53, implying that the status of glycolysis is sensed by proteins regulating insulin action. An enhanced glycolytic metabolic tone may alter production of reactive oxygen species, a known modulator of insulin action54, 63, 64 and insulin signalling65. Alternatively, recent approaches to map protein-metabolite interactions could identify points of allosteric crosstalk between glycolytic metabolites and insulin action proteins66, 67, hence broadening our understanding of the bidirectional communication between insulin action and metabolism.
It was striking that only several dozen insulin-regulated phosphopeptides associated significantly with glucose uptake. Since insulin triggers multiple distinct cellular outcomes, it is possible that only a subset of insulin-responsive phosphosites contribute to enhanced glucose uptake. Moreover, many of these phosphosites might only be permissive for insulin-stimulated glucose uptake and are not major regulatory nodes determining the fidelity of the process. For example, while mutation of the four primary Akt regulatory sites on Tbc1d4 blocks GLUT4 translocation68, none of these phosphosites featured strong positive correlations with glucose uptake in our analysis (Fig. S6). This implies that their phosphorylation may promote glucose uptake in a functionally permissive, switch-like fashion. We suggest that the glucose uptake-associated phosphosites we have identified will be enriched in major regulators of insulin responsiveness, necessitating future functional studies to characterise these sites and explore their involvement in insulin resistance.
Our work demonstrates that genetic and environmental variation can profoundly modulate global signalling networks and that the influences of these factors are intrinsically entwined. We show that the resulting diversity in signalling responses can be leveraged to pinpoint regulators of insulin-stimulated glucose uptake, providing a powerful methodological framework for interrogating the regulatory basis of complex biological pathways.
Limitations of this study
First, this study focused on male mice and examined only five inbred strains. This limited number of strains may mean that our association analysis was underpowered to detect some regulators of insulin responsiveness. Importantly, however, this does not imply that the regulators identified are incorrect, but only that there may be more to discover with larger cohorts. Future work should therefore extend our approach across a broader range of genetic backgrounds, as well as in female mice. Second, we only examined insulin signalling after 10 minutes, since measuring multiple timepoints would have drastically increased the number of animals and samples required. Integration of dynamic phosphoproteome data from cultured cells indicated that insulin signalling dynamics may contribute to trends in our data (Fig. 5i), suggesting the exploration of signalling at additional timepoints may be fruitful in the future. Third, mammalian tissues are a heterogeneous mixture of cell types, and differences in this mixture could result in different signalling responses measured at the whole tissue level. In our experience, the soleus can be reproducibly dissected as an intact muscle with little contamination from surrounding tissues, making it unlikely that cell type composition varied across samples due to tissue collection. However, we cannot exclude the possibility that differences in the composition of the soleus muscle across strains and diets contributed to the signalling changes we detected. Lastly, as we did not perform total proteomics in parallel to phosphoproteomics, we did not assess whether phosphosite changes were caused by differences in total protein abundance. While we cannot ignore the influence of changes in protein abundance, in our previous studies of insulin signalling in adipocytes31 or human skeletal muscle21 in which deep proteomes were measured in parallel, we found little global correlation between changes in protein phosphorylation and protein abundance in both unstimulated and insulin-stimulated conditions, suggesting the contribution of protein abundance to phosphosite changes is likely minimal.
D.E.J. is an Australian Research Council (ARC) Laureate Fellow. The content is solely the responsibility of the authors and does not necessarily represent the official views of the ARC. The authors also acknowledge the facilities, and the scientific and technical assistance of the Sydney Mass Spectrometry Facility and the Laboratory Animal Services at the Charles Perkins Centre, University of Sydney.
Conceptualization: JvG, SWCM, SJH, DEJ. Methodology: JvG, HBC, SWCM, ADV, MP, SJH. Formal analysis: JvG. Investigation: JvG, SWCM, HBC, ADV, MP, JS, SM, MEN, SJH. Resources: JS. Writing – Original Draft: JvG. Writing – Review and Editing: All authors. Visualization: JvG. Supervision: SWCM, SJH, DEJ. Project Administration: JvG, SWCM, SJH, DEJ. Funding Acquisition: DEJ
Declaration of Interests
The authors declare no competing interests.
Data Access Statement
All raw and Spectronaut processed phosphoproteomics data have been deposited in the PRIDE proteomeXchange repository and will be made publicly available upon publication. Processed data are available as supplementary tables.
All code used to analyse data and produce figures has been uploaded to https://github.com/JulianvanGerwen/GxE_muscle_phos
Most statistical analysis was performed in the R programming environment using RStudio (R version: 4.2.1, RStudio version: 2022.07.1 Build 554). Analysis of GLUT4-HA-L6 myotube Pfkfb3 expression, 2DG uptake, and ECAR was performed in GraphPad Prism (version: 9.3.1).
Male C57BL/6J (C57Bl6J), BXH9/TyJ (BXH9), BXD34/TyJ (BXD34), and CAST/EiJ (CAST) mice were purchased from Australian BioResources (Moss Vale, NSW, Australia) while NOD/ShiLtJ (NOD) mice were purchased from Animal Resources Centre (Murdoch, WA, Australia). Mice were at most 9 weeks old upon arrival. Mice were housed at 23 °C on a 12 h light/dark cycle in cages of 2-5, with free access to food and water. At 12-16 weeks of age mice were randomly allocated to a standard CHOW diet (13% calories from fat, 65% from carbohydrate, 22% from protein; “Irradiated Rat and Mouse Diet”, Specialty Feeds, Glen Forrest, WA, Australia) or a high-fat high-sucrose diet made in house (HFD; 45% calories from fat (40% calories from lard), 35% from carbohydrate (14% calories from starch), and 22% from protein) and sacrificed exactly 6 weeks later. The number of mice in each group are C57Bl6J: 8 CHOW, 10 HFD; NOD: 10 CHOW, 10 HFD; BXH9: 8 CHOW, 9 HFD; CAST: 9 CHOW, 9 HFD; BXD34: 10 CHOW, 11 HFD. Procedures were carried out with approval from the University of Sydney Animal Ethics Committee following guidelines issued by NHMRC (Australia).
Assessment of body composition
Body composition of individual mice was assessed using the EchoMRI-900 to determine lean mass 1 day before a glucose tolerance test and 5-6 days before euthanasia. Analysis was performed as per the manufacturer’s specifications.
Glucose tolerance test
On the day of a glucose tolerance test mice were fasted for 6 h (0800-1400). Mice were then orally gavaged with 20% (w/v) glucose in water at 2 g/kg lean mass, and blood glucose was measured from the tail vein using a glucometer 0, 15, 30, 45, 60, and 90 min after the gavage. At 0 and 15 min, 5 µL blood was also collected into an Insulin Mouse Ultra-Sensitive ELISA plate (Crystal Chem USA, Elk Grove Village, Illinois, USA). Blood insulin concentration was measured according to the manufacturer’s protocol, using linear extrapolation from an insulin standard curve. The area of the blood glucose curve (AOC) was calculated by:
Where i represents the ith timepoint at which blood glucose was measured, n represents the last timepoint, ti represents the time (min) of the ith timepoint, and Gi represents blood glucose concentration (mM) at the ith timepoint.
In vivo insulin stimulation
On the day of the procedure mice were fasted for 2h (1100-1300). Mice were then anaesthetised by intraperitoneal injection of sodium pentobarbital at 80 mg/kg body mass. To counter anaesthesia-associated declines in body temperature, mice were wrapped in aluminium foil and placed on a heating pad at 37°C. After 15 min anaesthesia, mice were injected retro-orbitally as previously described69 with 50 µL plasma replacement (B. Braun, Melsungen, DEU) containing 10 µCi 3H-2DG and saline or insulin (0.75 U/kg lean mass). Blood glucose was measured from the tail vein using a glucometer (AccuCheck, Roche Diabetes Care, NSW, Australia) 1 min prior to injection and 1, 5, 7.5, and 10 min after injection. Simultaneously, 5 µL blood was collected into 95 µL 0.9% NaCl on ice to measure 3H-2DG blood concentration. Ten minutes after insulin injection mice were sacrificed by cervical dislocation and the soleus muscle was rapidly excised, immediately frozen in liquid nitrogen, and stored at -80°C. To measure 3H-2DG blood concentration, diluted blood samples were first centrifuged at 10,000 xg for 10 min to pellet blood cells. Supernatant (70 µL) was collected and combined with 3 mL liquid scintillation cocktail (Perkin Elmer, Massachusetts, USA: 6013321) to allow the measurement of 3H with a Tri-Carb 2810TR Liquid Scintillation Counter (Perkin Elmer, Massachusetts, USA).
Skeletal muscle lysis
Frozen muscle tissue was powdered by grinding in a mortar and pestle chilled with liquid nitrogen and dry ice. To lyse powdered tissue, 200 µL lysis buffer (4% (w/v) sodium deoxycholate, 100 mM Tris-HCl pH 8.5) was added followed by 10 s vortexing. Samples were then sonicated at 4 °C at 90% power using pulses of 2 s on, 5 s off for a total time of 1 min. Samples were then immediately boiled at 95 °C with 1,500 rpm shaking for 5 min and sonicated for a further 2 min (4 °C, 90% power, 5 s on and 5 s off) to ensure complete lysis. Lysate was then centrifuged at 20,000 xg for 5 min and 180 µL supernatant was collected. Cysteine residues were reduced and alkylated by adding 40 mM chloroacetamide and 10 mM TCEP at pH 7. Lysate was incubated for 5 min at 45 °C with 1,500 rpm shaking and then incubated for a further 40 min at room temperature without shaking.
Next, 800 µL 100% chloroform and 1,600 µL 100% methanol were added following 30 s sonication at 90% power. LC/MS grade water (800 µL) was added following 5 min mixing at 1,000 rpm. Lysate was centrifuged for 5 min at 2,000 xg to induce a phase separation. The majority of the aqueous phase (2,400 µL) was removed and 2,000 µL was reserved for 3H-2DG quantification. Next, 2,400 µL 100% methanol was added following 30 s mixing at 800 rpm and centrifugation at 2,000 xg for 5 min. The supernatant was discarded, and the protein pellet was air-dried for 5 min. Protein was reconstituted in 200 µL lysis buffer, sonicated at 60% power for 15 s using a tip-probe sonicator and boiled for 5 min in a thermomixer at 95 °C with 1,500 rpm shaking.
Determining 3H-2DG uptake into muscle tissue
Anion exchange chromatography was used to quantify 3H-p2DG, representing 3H-2DG that has been taken up by cells. For quantification of total (phosphorylated and unphosphorylated) 3H-2DG, 375 µL lysate aqueous phase was combined with 1,125 µL water and reserved. For quantification of unphosphorylated 3H-2DG, 1,000 µL lysate aqueous phase was added to 300 µL 37.5% (w/v) AG 1-X8 anion exchange resin (Bio-Rad, Hercules, CA, USA: 1401441) and washed with 3 mL water to remove p2DG. Liquid scintillation cocktail (3 mL) was then added to 1,500 µL total and unphosphorylated 3H-2DG solutions, and 3H-2DG was measured using a Tri-Carb 2810R Liquid Scintillation Counter. Unphosphorylated and total 3H-2DG scintillation counts were subtracted to quantify 3H-p2DG.
3H-2DG blood concentration at 1, 5, 7.5, and 10 min after injection was fit to an exponential curve 𝑦 = 𝐶𝑝(0)𝑒−𝐾𝑝𝑡 where 𝐶𝑝(0) represents the predicted initial tracer concentration (DPM/µL) and 𝐾𝑝 represents the rate of tracer disappearance from the blood (1/min), to model the disappearance of 3H-2DG from the blood as it is taken up and trapped by peripheral tissues70. 𝐶𝑝(1) was removed when it was abnormally low (𝐶𝑝(1) < 𝐶𝑝(5), 𝐶𝑝(5) − 𝐶𝑝(1) > 0.5 × (𝐶𝑝(5) − 𝐶𝑝(7.5))), which likely indicates insufficient diffusion of circulating 3H-2DG into the tail vein. To account for different rates of blood 3H-2DG disappearance across mice, 3H-2DG uptake was calculated as a rate constant70:
Where 𝑡 represents the time after injection that the animal was sacrificed (min) and 𝐶𝑖 (𝑡) represents the concentration of 3H-p2DG in the tissue harvested at time 𝑡 (DPM/mg tissue).
Phosphoproteomics sample preparation
Phosphopeptides were isolated using the EasyPhos protocol27 with minor modifications. Protein (C57Bl6J and NOD: 755 µg, BXH9 and BXD34: 511 µg, CAST: 364 µg) was digested into peptides by incubation in 1% (w/w) Trypsin and LysC on a thermomixer at 37°C with 1,500 rpm shaking for 14 h. Following digestion, 400 µL 100% isopropanol and 100 µL EasyPhos enrichment buffer (48% (v/v) TFA, 8 mM KH2PO4) were sequentially added with mixing (1,500 rpm, 30 s) after each addition. Lysate was centrifuged at 20,000 xg for 15 min to pellet insoluble material and transferred to a deep well plate. The EasyPhos protocol was then followed from step 1227.
Liquid chromatography-tandem mass spectrometry (LC-MS/MS)
Enriched phosphopeptides in MS loading buffer (2% ACN, 0.3% TFA) were loaded onto in-house fabricated 55 cm columns (75 µM I.D.), packed with 1.9 µM C18 ReproSil Pur AQ particles (Dr. Maisch HPLC GmbH, Ammerbuch, DEU) with a Dionex U3000 HPLC (Thermo Fisher Scientific), interfaced with an Orbitrap Exploris 480 mass spectrometer (Thermo Fisher Scientific). Column temperature was maintained at 60°C using a column oven (Sonation lab solutions, Biberach, DEU), and peptides were separated using a binary buffer system comprising 0.1% formic acid (buffer A) and 80% ACN plus 0.1% formic acid (buffer B) at a flow rate of 400 nL/min. A gradient of 3–19% buffer B was employed over 40 min followed by 19–41% buffer B over 20 min, resulting in approximately 1 h gradients. Peptides were analysed with one full scan (350–1,400 m/z, R = 120,000) at a target of 3e6 ions, followed by 48 data-independent acquisition MS/MS scans (350–1,022 m/z) with higher-energy collisional dissociation (target 3e6 ions, max injection time 22 ms, isolation window 14 m/z, 1 m/z window overlap, normalised collision energy 25%), and fragments were detected in the Orbitrap (R = 15,000).
MS raw data processing
Raw spectral data were analysed using Spectronaut (v16.0.220606.53000). Data were searched using directDIA against the Mouse UniProt Reference Proteome database (January 2022 release), with default settings (precursor and protein Qvalue cutoffs 0.01, Qvalue filtering, MS2 quantification), with “PTM localization” filtering turned on (threshold 0.5), and the inbuilt peptide collapse function.
Phosphoproteomics data processing
Phosphopeptide intensities were log2 transformed and median normalised. Non-class I phosphopeptides (localisation score ≤ 0.75) were then removed. Finally, for each phosphopeptide, outlier values were removed that had a log2 intensity > 5 and were < 6 log2 intensity units lower than the phosphopeptide median. Log2 fold changes between conditions were computed using condition medians.
Identifying insulin-regulated phosphopeptides
To allow comparison across conditions, phosphopeptides were filtered for those highly quantified in most strain-diet combinations. For a given phosphopeptide, this filtering was performed on two levels. Firstly, each of the 10 strain-diet combinations were retained if there were ≥ 3 insulin-stimulated values and ≥ 3 unstimulated values. Then, the phosphopeptide itself was retained if ≥ 8 strain-diet combinations had passed the previously filtering. Phosphopeptides were then fit to a three-way ANOVA with all interaction terms (“aov” in the R package “stats”) and an F-test was performed assessing the main effect of insulin stimulation. To correct for multiple hypothesis testing p-values were converted into q-values (R package “qvalue”71). The log2(insulin/unstimulated) fold change with the greatest magnitude across strain-diet combinations was then calculated (max log2(insulin/unstimulated)). Phosphopeptides were considered “insulin-regulated” when q < 0.05 and if insulin increased or decreased phosphorylation by > 1.5-fold in at least one strain-diet combination (max log2(insulin/unstimulated) > 0.58 or < -0.58).
Calculation of insulin response values
For all phosphopeptides the distribution of “insulin responses” in each strain-diet combination was calculated. Specifically, within each strain-diet combination all insulin-stimulated values were normalised by subtracting the unstimulated median.
Identifying strain and diet effects
For each insulin-regulated phosphopeptide a one-way ANOVA was performed modelling the insulin response as a function of mouse strain within the CHOW diet. p-values were converted to q-values. For significant phosphopeptides (q < 0.05), t-tests were performed comparing the insulin response of C57Bl6J to each of the other four strains. t-test p-values were converted to q-values and considered significant when q < 0.05. To ensure that significant differences between a strain and C57Bl6J were of a meaningful magnitude, the strain’s log2(insulin/unstimulated) was compared to the C57Bl6J log2(insulin/unstimulated). In general, if the absolute difference between the two was greater than 0.58 this was accepted. However, this threshold was decreased for phosphopeptides with weaker insulin regulation. Specifically, the difference was retained if it passed the following filtering:
An insulin-regulated phosphopeptide was considered to have a “Strain effect” if the insulin response in at least one strain was different to C57Bl6J using the q-value and log2 fold-change criteria described above.
Uniform diet and StrainxDiet effects
For each insulin-regulated phosphopeptide a two-way ANOVA was performed modelling the insulin response as a function of strain, diet, and their interaction. The p-values for the Diet and StrainxDiet terms were converted to q-values. Whenever the StrainxDiet term was significant (q < 0.05), additional tests were performed to identify specific strains in which the insulin response differed between CHOW and HFD. If the StrainxDiet term was not significant but the Diet term was significant, a separate filtering procedure was performed.
When the StrainxDiet term was significant, t-tests were performed to compare the CHOW insulin response to the HFD insulin response within each strain. When a t-test was significant (q < 0.05), the log2(insulin/unstimulated) filtering procedure described for “Strain effects” was applied comparing CHOW and HFD fold changes. Insulin-regulated phosphopeptides were considered to have a “StrainxDiet effect” if there was a difference between CHOW and HFD in at least one strain.
When only the Diet term was significant, the log2(insulin/unstimulated) filtering procedure described for “Strain effects” was applied, comparing the mean log2(insulin/unstimulated) across strains within CHOW, to the mean across HFD. Insulin-regulated phosphopeptides that passed this filter were considered to have a “Uniform diet effect”.
Curated insulin signalling subnetwork
A subnetwork of insulin-regulated phosphosites was curated by compiling all sites on proteins from a previously published knowledge pathway-derived list of canonical insulin signalling proteins21. Several phosphosites and proteins that were not detected as insulin regulated were included due to their importance in the insulin signalling pathway. Additionally, all in vivo substrates of canonical insulin-regulated kinases (Akt, mTOR, AMPK, Raf, Mek1/2, Erk1/2, p90RSK/Rsk2, p70S6K, Pdk1, INSR) annotated in PhosphositePlus were included32. Annotations from orthologous phosphosites were pooled across species using PhosphositePlus Site Group IDs. Phosphosite regulatory roles from PhosphositePlus were indicated after manual validation by literature search. Proteins were assigned to functional groups (e.g. mRNA processing, lipid metabolism) based on their Uniprot descriptions.
Kinase substrate enrichment analysis
Kinase-substrate annotations were collated from PhosphositePlus and mapped into phosphoproteomics data using Site Group IDs. Only annotations supported by in vivo evidence were used. Annotations for kinase isoforms (e.g. Akt1, Akt2, Akt3) were merged. Substrate annotations for GSK3 were supplemented with a recent list of 274 putative GSK3 substrates determined by phosphoproteomics and motif analysis31. Autophosphorylation sites and promiscuous phosphosites targeted by ≥ 4 kinases were removed. Kinase substrate enrichment analysis (KSEA) was then performed with the “ksea” function from the R package “ksea”37 (version: 0.1.2) using insulin response data and 1,000 permutations to determine empirical p-values. Only phosphopeptides quantified in ≥ 50% of samples and with ≥ 1 insulin response value in all strain-diet combinations were used. In each sample kinases with < 5 quantified substrates were excluded, and only kinases with significant enrichment (p < 0.05) in ≥ 5 samples were used in subsequent analysis. To identify Strain effects on kinase activity, one-way ANOVAs were performed on CHOW KSEA enrichment scores. To identify Uniform diet or StrainxDiet effects, two-way ANOVAs were performed on KSEA enrichment scores testing the effects of strain, diet, and their interaction. p-values were adjusted by the Benjamini-Hochberg procedure.
Insulin signalling subnetwork analysis
Weighted gene correlation network analysis (WGCNA39, 40) was performed with the “blockwiseModules” function from the R package “WGCNA” (version 1.71) using the insulin response values of all insulin-regulated phosphopeptides. Default parameters were used except for power = 3 (determined as recommended in 39), deepSplit = 3, minModuleSize = 15, reassignThreshold = 0, and mergeCutHeight = 0.25. Subnetwork eigengenes were extracted and termed “eigenpeptides”. One-sided fisher’s exact tests were performed to assess the enrichment of Gene Ontology (GO)
Biological Processes, GO Cellular Compartments (R package “org.Mm.eg.db” version 3.15.072), and kinase substrates in each subnetwork relative to the entire phosphoproteome. Only pathways containing three or more subnetwork phosphoproteins were tested. Kinase substrate enrichment was performed using the same annotations as KSEA. P-values were adjusted within each analysis by the Benjamini-Hochberg procedure. Subnetwork phosphopeptides were mapped into insulin signalling temporal clusters defined in our previous study of insulin signalling dynamics10, using PhosphositePlus Site Group IDs. The timepoint at which each cluster appeared to reach its maximum insulin-stimulated value was used as a measure of insulin response speed.
Glucose uptake correlations
For each insulin-regulated phosphopeptide, Pearson’s correlation tests were performed to assess the linear association between 3H-2DG uptake in insulin-stimulated mice and phosphopeptide insulin response values or unnormalized insulin-stimulated log2 intensity. Phosphopeptides were considered correlated with 3H-2DG uptake when q < 0.1 and their Pearson’s correlation coefficient was of substantial magnitude (r > 0.35 or r < -0.35). Pearson’s correlation tests were also performed comparing insulin-stimulated 3H-2DG uptake to KSEA enrichment scores in individual mice or using the median in each strain-diet combination.
GLUT4-HA-L6 myoblasts73 were grown in α-MEM supplemented with 10% fetal bovine serum in a humidified chamber at 37 °C, 10% CO2. Differentiation was induced by changing media to α-MEM supplemented with 2% horse serum for 5 days.
Platinum-E (Plat-E) retroviral packaging cells were grown to confluency and transfected with 10 μg total DNA: either pBabe puromycin empty vector, pBabe puromycin PFKFB3 or pWZL neomycin HA-GLUT4. After 48 h retroviral media was collected and passed through a 0.45 μm filter. L6 myotubes were grown to confluence and retrovirally transfected with 2 mL of HA-GLUT4 neomycin viral media in the presence of 10 μg/ml polybrene. The following morning, cells were split into growth media containing neomycin (800 μg/ml) and passaged until only successfully transfected cells remained. These cells were then grown to confluence and retrovirally transfected again with 2 mL of either empty vector puromycin viral media or PFKFB3 puromycin viral media in the presence of 10 μg/ml of polybrene. The following morning, cells were split into growth media containing both neomycin (800 μg/ml) and puromycin (2 μg/ml) and passaged until only successfully transfected cells remained in culture.
GLUT4-HA-L6 myotubes were incubated overnight (16 h) in either BSA-conjugated 125 µM palmitate or BSA vehicle control. Cells were then washed in ice-cold PBS and lysed by scraping directly into 55 °C Laemmli sample buffer with 10 % (tris 2-carboxyethyl phosphine; TCEP). Samples were sonicated for 24 s (3s on/3s off) and heated at 65 °C for 5 minutes. Samples were then resolved by SDS-PAGE as previously described17, transferred onto PVDF membranes and blocked in TBS-T (0.1% Tween in Tris-buffered saline) containing 5% skim milk for 1 h. Membranes were then washed 3 x 10 min in TBS-T and incubated overnight in primary antibodies against Pfkfb3 (Proteintech Group; 13763-1-AP) and α-tubulin (Cell Signalling Technologies #2125; diluted 1:1000). The following day membranes were washed 3 x 10 min in TBS-T and incubated for 1 h in species-appropriate fluorescent antibodies. Imaging and densitometry were performed using LI-COR Image Studio.
Extracellular acidification rate
The extracellular acidification rate (ECAR) of GLUT4-HA-L6 cells myotubes was measured using Seahorse XFp miniplates and a Seahorse XF HS Mini Analyzer (Seahorse Bioscience, Copenhagen, Denmark) as previously described74. Cells incubated in palmitate or BSA control were washed twice with Krebs-Ringer Bicarbonate Buffer (Sigma, K4002) and once with standard cell culture media without bicarbonate (XF-DMEM, pH 7.4). Cells were then incubated in XF-DMEM without glucose at 37°C for 1 h in a non-CO2 incubator, followed by assaying in the XFp Analyzer. ECAR was measured after a 12-minute equilibration period followed by mix/wait/read cycles of 3/0/3 min. After stabilizing the baseline rates, compounds were injected to reach a final concentration of: 10 mM glucose, 5 μg/mL oligomycin, and 50 mM 2-deoxyglucose (2-DG), allowing estimation of glucose-driven glycolysis (glucose ECAR - basal ECAR), glycolytic capacity (oligomycin ECAR - 2DG ECAR), and non-glycolytic acidification (equal to 2DG ECAR). Data were normalized to protein concentration and presented as a percentage of total ECAR.
2DG uptake in GLUT4-HA-L6 myotubes
2-deoxyglucose (2DG) uptake into GLUT4-HA-L6 myotubes was performed as previously described with modifications73, 75. Cells were incubated overnight in αMEM supplemented with either BSA-coupled 125 µM palmitic acid or BSA vehicle control before being washed 3x with 37 °C HEPES-buffered saline (HBS). Cells were then incubated in HBS supplemented with 10 µM unlabelled 2-deoxyglucose and either 0 or 100 nM insulin at 37°C for 15 min. Cells were then incubated for a further 5 min following the addition of 0.5 µCi/ml [3H]-2-deoxyglucose in HBS. Cells were then washed on ice 5x with ice-cold PBS and lysed in 1 M NaOH. For non-specific background uptake, 1 well per condition was pre-treated with cytochalasin B. Counts were determined by Perkin Elmer Quantulus GCT Liquid Scintillation Counter (Perkin Elmer, Waltham, MA, USA). 2DG uptake was expressed relatively to protein concentration as determined by bicinchoninic acid (BCA) assay after neutralisation with 1 M HCl and subtraction of non-specific uptake.
Table S1: Muscle phosphoproteomics
(Page 1 “01_quantification”) Normalized LFQ intensities of class I phosphopeptides. (Page 2 “02_analysis”) Statistical analysis of phosphoproteome data.
Table S2: Insulin signalling subnetworks
WGCNA-derived subnetworks of insulin-regulated phosphopeptides.
- 1.Protein Phosphorylation: A Major Switch Mechanism for Metabolic RegulationTrends Endocrinol. Metab 26:676–687
- 2.Illuminating the dark phosphoproteomeSci. Signal 12https://doi.org/10.1126/scisignal.aau8645
- 3.Global, in vivo, and site-specific phosphorylation dynamics in signaling networksCell 127:635–648
- 4.Phosphoproteomic analysis reveals interconnected system-wide responses to perturbations of kinases and phosphatases in yeastSci. Signal 3
- 5.Phosphoproteomics of Acute Cell Stressors Targeting Exercise Signaling Networks Reveal Drug Interactions Regulating Protein SecretionCell Rep 29:1524–1538
- 6.The regulatory landscape of the yeast phosphoproteome
- 7.Biochemical and cellular properties of insulin receptor signallingNat. Rev. Mol. Cell Biol 19:31–44
- 8.MAPK signalling in cellular metabolism: stress or wellness?EMBO Rep 11:834–840
- 9.Rac1 signaling is required for insulin-stimulated glucose uptake and is dysregulated in insulin-resistant murine and human skeletal muscleDiabetes 62:1865–1875
- 10.Dynamic adipocyte phosphoproteome reveals that Akt directly regulates mTORC2Cell Metab 17:1009–1020
- 11.High-throughput phosphoproteomics reveals in vivo insulin signaling dynamicsNat. Biotechnol 33:990–995
- 12.Dissection of the insulin signaling pathway via quantitative phosphoproteomicsProc. Natl. Acad. Sci. U. S. A 105:2451–2456
- 13.The aetiology and molecular landscape of insulin resistanceNat. Rev. Mol. Cell Biol https://doi.org/10.1038/s41580-021-00390-6
- 14.Muscle and adipose tissue insulin resistance: malady without mechanism?J. Lipid Res 60:1720–1732
- 15.Gene–environment interactions in human diseasesNat. Rev. Genet 6:287–298
- 16.Systems genetics approaches to understand complex traitsNat. Rev. Genet 15:34–48
- 17.Systems-level analysis of insulin action in mouse strains provides insight into tissue- and pathway-specific interactions that drive insulin resistanceCell Metab 34:227–239
- 18.Mouse strain-dependent variation in obesity and glucose homeostasis in response to high-fat feedingDiabetologia 56:1129–1139
- 19.The impact of genomic variation on protein phosphorylation states and regulatory networksMol. Syst. Biol 18
- 20.Multi-omics analysis identifies drivers of protein phosphorylationGenome Biol 24
- 21.Personalized phosphoproteomics identifies functional signalingNat. Biotechnol https://doi.org/10.1038/s41587-021-01099-9
- 22.Genome-wide association studies for complex traits: consensus, uncertainty and challengesNat. Rev. Genet 9:356–369
- 23.The Hybrid Mouse Diversity Panel: a resource for systems genetics analyses of metabolic and cardiovascular traitsJ. Lipid Res 57:925–942
- 24.High dietary fat and sucrose result in an extensive and time-dependent deterioration in health of multiple physiological systems in miceJ. Biol. Chem 293:5731–5745
- 25.The Triumvirate: β-Cell, Muscle, Liver: A Collusion Responsible for NIDDMDiabetes 37:667–687
- 26.Fiber types in mammalian skeletal musclesPhysiol. Rev 91:1447–1531
- 27.High-throughput and high-sensitivity phosphoproteomics with the EasyPhos platformNat. Protoc 13:1897–1916
- 28.Rapid and site-specific deep phosphoproteome profiling by data-independent acquisition without the need for spectral librariesNat. Commun 11
- 29.A Cell-Autonomous Signature of Dysregulated Protein Phosphorylation Underlies Muscle Insulin Resistance in Type 2 DiabetesCell Metab 32:844–859
- 30.Signaling defects associated with insulin resistance in non-diabetic and diabetic individuals and modification by sexJ. Clin. Invest https://doi.org/10.1172/JCI151818
- 31.Phosphoproteomics reveals rewiring of the insulin signaling network and multi-nodal defects in insulin resistanceNat. Commun 14:1–20
- 32.PhosphoSitePlus, 2014: mutations, PTMs and recalibrationsNucleic Acids Res 43:D512–20
- 33.The functional landscape of the human phosphoproteomeNat. Biotechnol 38:365–373
- 34.Overexpression of vesicle-associated membrane protein (VAMP) 3, but not VAMP2, protects glucose transporter (GLUT) 4 protein translocation in an in vitro model of cardiac insulin resistanceJ. Biol. Chem 287:37530–37539
- 35.ZRF1 is a novel S6 kinase substrate that drives the senescence programmeEMBO J 36:736–750
- 36.Oxysterols inhibit phosphatidylcholine synthesis via ERK docking and phosphorylation of CTP:phosphocholine cytidylyltransferaseJ. Biol. Chem 280:21577–21587
- 37.Benchmarking substrate-based kinase activity inference using phosphoproteomic dataBioinformatics 33:1845–1851
- 38.Impaired Akt phosphorylation in insulin-resistant human muscle is accompanied by selective and heterogeneous downstream defectsDiabetologia 56:875–885
- 39.A General Framework for Weighted Gene Co-Expression Network AnalysisStatistical Applications in Genetics and Molecular Biology 4https://doi.org/10.2202/1544-6115.1128
- 40.WGCNA: an R package for weighted correlation network analysisBMC Bioinformatics 9
- 41.Insulin stimulation of glucose uptake in skeletal muscles and adipose tissues in vivo is NO dependentAm. J. Physiol 274:E692–9
- 42.Role of Nitric Oxide in Insulin Secretion and Glucose MetabolismTrends Endocrinol. Metab 31:118–130
- 43.Activation of nitric oxide synthase in endothelial cells by Akt-dependent phosphorylationNature 399:601–605
- 44.Afadin is a scaffold protein repressing insulin action via HDAC6 in adipose tissueEMBO Rep 20
- 45.AMPK and Exercise: Glucose Uptake and Insulin SensitivityDiabetes Metab. J 37:1–21
- 46.Phosphorylation and activation of heart 6-phosphofructo-2-kinase by protein kinase B and other protein kinases of the insulin signaling cascadesJ. Biol. Chem 272:17269–17275
- 47.Phosphorylation and activation of heart PFK-2 by AMPK has a role in the stimulation of glycolysis during ischaemiaCurr. Biol 10:1247–1255
- 48.Balancing glycolytic flux: the role of 6-phosphofructo-2-kinase/fructose 2,6-bisphosphatases in cancer metabolismCancer Metab 1
- 49.GLUT4 exocytosisJ. Cell Sci 124:4147–4159
- 50.Acute selective glycogen synthase kinase-3 inhibition enhances insulin signaling in prediabetic insulin-resistant rat skeletal muscleAm. J. Physiol. Endocrinol. Metab 288:E1188–94
- 51.Short-term in vitro inhibition of glycogen synthase kinase 3 potentiates insulin signaling in type I skeletal muscle of Zucker Diabetic Fatty ratsMetabolism 56:931–938
- 52.Selective glycogen synthase kinase 3 inhibitors potentiate insulin activation of glucose transport and utilization in vitro and in vivoDiabetes 52:588–595
- 53.Kinome Screen Identifies PFKFB3 and Glucose Metabolism as Important Regulators of the Insulin/Insulin-like Growth Factor (IGF)-1 Signaling PathwayJ. Biol. Chem 290:25834–25846
- 54.Insulin resistance is a cellular antioxidant defense mechanismProc. Natl. Acad. Sci. U. S. A 106:17787–17792
- 55.Multiple levels of cyclin specificity in cell-cycle controlNat. Rev. Mol. Cell Biol 8:149–160
- 56.New insights into activation and function of the AMPKNat. Rev. Mol. Cell Biol 24:255–272
- 57.ERK signalling: a master regulator of cell behaviour, life and fateNat. Rev. Mol. Cell Biol 21:607–632
- 58.Ribosomal Protein S6 Phosphorylation: Four Decades of ResearchInt. Rev. Cell Mol. Biol 320:41–73
- 59.An atlas of substrate specificities for the human serine/threonine kinomeNature 613:759–766
- 60.Trends in kinase drug discovery: targets, indications and inhibitor designNat. Rev. Drug Discov 20:839–861
- 61.Patient-specific logic models of signaling pathways from screenings on cancer biopsies to prioritize personalized combination therapiesMol. Syst. Biol 16
- 62.Patient-specific Boolean models of signalling networks guide personalised treatmentsElife 11https://doi.org/10.7554/eLife.72626
- 63.Reactive oxygen species have a causal role in multiple forms of insulin resistanceNature 440:944–948
- 64.Mitochondrial oxidative stress causes insulin resistance without disrupting oxidative phosphorylationJ. Biol. Chem 293:7315–7328
- 65.Global redox proteome and phosphoproteome analysis reveals redox switch in AktNat. Commun 10
- 66.A Map of Protein-Metabolite Interactions Reveals Principles of Chemical CommunicationCell 172:358–372
- 67.Protein-metabolite interactomics of carbohydrate metabolism reveal regulation of lactate dehydrogenaseScience 379:996–1003
- 68.Insulin-stimulated phosphorylation of a Rab GTPase-activating protein regulates GLUT4 translocationJ. Biol. Chem 278:14599–14602
- 69.Retro-orbital injections in miceLab Anim 40:155–160
- 70.Investigation of the effect of insulin upon regional brain glucose metabolism in the rat in vivoEndocrinology 107:1827–1832
- 71.A direct approach to false discovery ratesJ. R. Stat. Soc. Series B Stat. Methodol 64:479–498
- 72.org. Mm. eg. db: Genome wide annotation for Mouse. R package version 3.2. 3. Bioconductor
- 73.Interleukin-6 increases insulin-stimulated glucose disposal in humans and glucose uptake and fatty acid oxidation in vitro via AMP-activated protein kinaseDiabetes 55:2688–2697
- 74.Mitochondrial oxidants, but not respiration, are sensitive to glucose in adipocytesJ. Biol. Chem 295:99–110
- 75.β-Catenin is required for optimal exercise- and contraction-stimulated skeletal muscle glucose uptakeJ. Physiol 599:3897–3912
- 1.Personalized phosphoproteomics identifies functional signalingNat. Biotechnol https://doi.org/10.1038/s41587-021-01099-9
- 2.PhosphoSitePlus, 2014: mutations, PTMs and recalibrationsNucleic Acids Res 43:D512–20