Participant characteristics

Table 1 shows characteristics of the 119 participants of which 89 were female and 30 male. The mean (SD) age was 78.24 (2.59), body mass index (BMI) (kg/m2) 27.25 (4.06), HOMA2-IR 1.03 (0.58), fasting insulin (mU/L) 7.73 (4.42), fasting glucose (mmol/l) 5.43 (0.65) and HbA1c 37.79 (5.33). Of the 119 participants with DNA methylation data, 115 had phenotypic data for HOMA2-IR and 114 for HbA1c. Of the HOMA2-IR characterised participants, 115 had phenotypic data for grip strength, 111 for appendicular lean mass (ALMi) and 115 for gait speed. For HbA1c-characterised participants, 114 had phenotypic data for grip strength, 110 for ALMi and 114 for gait speed.

Identification of differentially methylated CpGs in myoblasts associated with HOMA2-IR

There were 38 differentially methylated CpGs (dmCpGs) significantly (FDR < 0.05) associated with HOMA2-IR (Fig. 1 A, B, C, D, E, Table 2, Supplementary Table 1). The three dmCpGs with the strongest associations were cg08944484, located within the 3′ untranslated region (UTR) of the saccharine dehydrogenase (SCCPDH) gene on chromosome 1 (FDR = 3.27 × 10−03), cg01351409 located in an intergenic region of chromosome 7 (FDR = 3.37 × 10−03) and cg20625334 found within the gene body of signal transducing adapter molecule 2 (STAM2) (FDR = 3.37 × 10−03).

Fig. 1

A Volcano plot showing the dmCpGs with respect to HOMA2-IR (FDR < 0.05 = red, top 10 FDR < 0.05 labelled). B Manhattan plot showing the dmCpGs with respect to HOMA2-IR (FDR < 0.05 = red, top 10 FDR < 0.05 labelled). C Locations of dmCpGs with respect to HOMA2-IR (FDR < 0.05). D Overlap of dmCpGs (FDR < 0.05) with respect to HOMA2-IR and HbA1c. E Overlap of top 100 dmCpG-associated genes with respect to HOMA2-IRand HbA1c

Table 2 HOMA2-IR-associated dmCpGs (FDR < 0.05) in skeletal muscle myoblast cells from elderly individuals

To gain an understanding of the functional significance of the methylation changes associated with HOMA2-IR, the top 100 dmCpG-associated genes were inputted into STRING to generate a PPI network. The PPI enrichment P value for the overall network for HOMA2-IR was 4.3 × 10−02 indicating biological connection between the proteins. Selection of the largest gene cluster identified within the network was found to be enriched for tissue homeostasis (FDR 6.60 × 10−04) anatomical structure homeostasis (FDR 7.50 × 10−04) and regulation of JNK cascade (FDR 1.90 × 10−03) within gene ontology (GO) pathways for biological process, and AMPK signalling (FDR 1.41 × 10−02) and Insulin signalling (FDR 1.41 × 10−02) within the KEGG ontology (Table 3, Supplementary Table 2). To identify potential causal networks and upstream regulators that may mediate changes in IR, the HOMA2-IR dmCpG-associated genes were analysed using IPA. The top networks were cellular movement, connective tissue development and function, and skeletal and muscular system development and function (Fig. 2, Supplementary Table 3), with the top upstream regulators identified as T cell factor (TCF) (p 3.80E − 04), Clock Circadian Regulator (CLOCK) (p 1.58 × 10−03) and MAPK (p 3.76 × 10−03) (Supplementary Table 4).

Table 3 HOMA2-IR-enriched pathways in myoblasts from aged skeletal muscle; top 10 Gene Ontology (GO) biological process and top 10 KEGG pathways shownFig. 2

Network diagram showing the top IPA gene network of cellular movement, connective tissue development and function and skeletal and muscular system development and function and associated molecules. Grey nodes indicate molecules within the input dataset, while white nodes are generated by IPA core analysis and represent molecules not in the input dataset added to the network from knowledge base to generate the interaction network

Influence of BMI on the HOMA2-IR associated DNA methylation signature

To assess the influence of BMI on the HOMA2-IR associated methylation signature in myoblasts, BMI was added as a covariate to the regression model; after adjustment for BMI, there were 6 dmCpGs associated (FDR ≤ 0.05) with HOMA2-IR. The strongest associations were with cg13451048 located within the gene body of exoribonuclease family member 3 (ERI3) on chromosome 1 (FDR 3.60 × 10−03), and cg01351409 and cg22007860 located in intergenic regions of chromosome 7 (FDR 4.20 × 10−03) and chromosome 3 (FDR 1.40 × 10−02), respectively (Supplementary Table 5). Of the 6 dmCpGs associated with HOMA2-IR after adjustment for BMI, 3 overlapped with the dmCpGs identified without adjustment for BMI, these being cg13451048 within ERI3, and the intergenic CpGs cg01351409 and cg10416784.

Pathway analysis of the top 100 dmCpG-associated genes after adjustment for BMI showed enrichment of regulation of epidermal growth factor receptor (FDR 1.01 × 10−02) and regulation of protein kinase B signalling (FDR 1.89 × 10−02) and anatomical structure morphogenesis (FDR 4.06 × 10−02) within the GO pathways for biological process, and relaxin signalling (FDR 5.98 × 10−05) and endocrine resistance (FDR 2.00 × 10−04) amongst the KEGG ontology (FDR 4.00 × 10−04) (Supplementary Table 6).

Influence of sarcopenia on HOMA2-IR-associated DNA methylation signature

As sarcopenia has been implicated in the development of insulin resistance, we adjusted the regression model for each of the 3 definitional components of sarcopenia, namely ALMi, grip strength and gait speed to determine how these measures of muscle mass, strength and function may influence the HOMA2-IR-associated DNA methylation signal in myoblasts (Supplementary Tables 7, 8 and 9). After addition of ALMi as a covariate in the model, the number of dmCpGs associated with HOMA2-IR was reduced from 38 to 7. The top dmCpGs were cg01351409, an intergenic CpG located on chromosome 7, cg18292904 located within the macrophage-stimulating protein receptor gene (MSTR1, also known as RON kinase) and cg14878421 an intergenic CpG located on chromosome 4 (Supplementary Table 7). Only 5 of the 38 dmCpGs identified in the unadjusted analysis remained significantly associated with HOMA2-IR; these were cg20625334 (STAM2), cg13451048 (ERI3) and the intergenic CpGs: cg01351409, cg10703001 and cg00623421, which all showed the same direction of association, effect size and significance. Adjustment for grip strength in the model also attenuated the methylation signal associated with HOMA2-IR with 15 dmCpGs associated with HOMA2-IR after adjustment; 12 of which overlapped with the dmCpGs from the unadjusted analysis (Supplementary Table 8), while the addition of gait speed to the model resulted in the association of 23 dmCpGs with HOMA2-IR, of which 18 of the dmCpGs overlapped with those identified in the unadjusted analyses (Supplementary Table 9).

PPI analysis revealed that the top pathways enriched amongst biological process GO category were regulation of protein kinase B signalling (FDR 2.20 × 10−05) and regulation of protein phosphorylation (FDR 4.78 × 10−05) after adjustment for ALMi (Supplementary Table 10), regulation of protein kinase B signalling (FDR 2.00 × 10−03) and positive regulation of protein kinase B signalling (FDR 4.00 × 10−03) after adjustment for grip strength (FDR 2.00 × 10−03) (Supplementary Table 11), and generation of neurons (FDR 1.60 × 10−03) and positive regulation of kinase activity (FDR 1.90 × 10−03) after adjustment for gait speed (Supplementary Table 12). The top KEGG pathways were MAPK (FDR 8.30 × 10−05) and insulin signalling (FDR 8.30 × 10−05) after adjustment for ALMi, Prolactin Signalling (FDR 2.90 × 10−03) and longevity regulation (FDR 1.34 × 10−02) after adjustment for grip strength, and insulin signalling (FDR4.70 × 10−03) and bacterial invasion of epithelial cells (FDR 4.70 × 10−03) after adjustment for gait speed.

HbA1c was associated with differential DNA methylation of skeletal muscle myoblasts

To identify methylation changes associated with longer term glycemia, we examined associations between DNA methylation and HbA1c levels. Forty-nine dmCpGs were associated with HbA1c (FDR < 0.05) (Supplementary Table 13), with the top 3 CpGs being cg19477361 located in a CpG island and 5′UTR region of guanine nucleotide-binding protein subunit gamma-7 (GNG7) (FDR 2.42 × 10−03), cg13451048 located in the gene body of exoribonuclease family member 3 (ERI3) on chromosome 1 (FDR 3.06 × 10−03), and cg22337620 located in the gene body of the DLGAP1 gene which encodes for disk large-associated protein 1 (DAP-1), also known as guanylate kinase-associated protein (GKAP) (FDR 3.81 × 10−03). Only one dmCpG associated with both HOMA2-IR and HbA1c, which was cg13451048 within ERI3 (Fig. 1D, Supplementary Table 14). After addition of BMI in the regression model, there were 48 dmCpGs associated with HbA1c (Supplementary Table 15), with 36 of the dmCpGs overlapping with the dmCpGs in the unadjusted model. The addition of ALMi as a covariate resulted in 44 dmCpGs associated with Hb1Ac, with 18 dmCpGs in common with those from the unadjusted analysis (Supplementary Table 16), while the addition of grip strength resulted in 70 dmCpGs associated with Hb1Ac levels with 37 in common with the unadjusted analyses (Supplementary Table 17); adjustment for gait speed led to 64 dmCpGs associated with HbA1c, with 38 overlapping with those identified in the unadjusted analysis (Supplementary Table 18).

Network analysis of HbA1c dmCpG-associated genes showed a PPI enrichment P value for the overall network of 1.38 × 10−02. Within the largest cluster the top GO biological process pathways were nervous system development (FDR 1.27 × 10−05), multicellular organism development (FDR 4.23 × 10−05) and cell differentiation (FDR 4.23 × 10−05), whilst in the KEGG ontology top pathways were Hippo signalling (FDR 1.60 × 10−03) and pathways in cancer (FDR 1.60 × 10−03) (Supplementary Table 19). The top pathways enriched within the biological process category, after the addition of BMI as a covariate, were cell communication (FDR 3.79 × 10−06), Signalling (FDR 3.79 × 10−06) and developmental process (FDR 3.79 × 10−06) (Supplementary Table 20). Top pathways enriched were neuron migration (FDR 8.17 × 10−05), cell communication (FDR 8.17 × 10−05) and signal transduction FDR 8.17 × 10−05) after addition of ALMi in the model (Supplementary Table 21); enzyme-linked receptor protein signalling pathways (FDR 4.90 × 10−03), cell differentiation (FDR 4.90 × 10−03) and signaling (FDR = 5.50 × 10−03) after adjustment for grip strength (Supplementary Table 22), and cell differentiation (FDR 8.98 × 10−06), generation of neurons (FDR 8.98 × 10−06) and cell communication (FDR 1.10 × 10−04) after adjustment for gait speed (Supplementary Table 23). Similar KEGG pathways were enriched in both the unadjusted and adjusted analyses.

HOMA2-IR and HbA1c are associated with multiple differentially methylated regions

Regional analysis identified DMRs associated with HOMA2-IR and HbA1c (Table 4, Supplementary Tables 24 and 25). 8 DMRs were associated with HOMA2-IR (Stouffer < 0.05), with the top 3 DMRs located within the Homobox A3 (HOXA3) gene, consisting of 8 CpGs (Stouffer 9.70 × 10−03), the T-box transcription factor (TBX1) gene consisting of 4 CpGs (Stouffer 1.0 × 10−03) and the nuclear receptor subfamily 2 group F member 2 (NR2F2) consisting of 4 CpGs (Stouffer 1.3 × 10−02), respectively (Table 4). After adjustment for BMI, 12 DMRs were associated with HOMA2-IR, with the top 3 located in CTD-2562J15.6, TBX1 and N2RF2, which were all associated with HOMA2-IR in the unadjusted analysis (Supplementary Table 26) Adjustment for each of the definitional components of sarcopenia resulted in 4 DMRs associated with HOMA2-IR; with DMRs being located within TBX1, RNF126P1, HOXA3 and NR2F2 after adjustment for either grip strength or gait speed (Supplementary Tables 27 and 28), and DMRs within NR2F2, RNF126P and TBX1 together with an intergenic region of Chr 20 after adjustment for ALMi (Supplementary Table 29). HbA1c was associated with 21 DMRs, with the top DMR located in the HLA Complex Group 9 (HGC9) gene (Stouffer 1.9 × 10−04) (Supplementary Table 25). There were no DMRs that overlapped between HOMA2-IR and HbA1c. There were 32 DMRs associated with HbA1c after adjustment for BMI (Supplementary Table 30), and 9 after adjustment for either grip strength or gait speed (Supplementary Tables 31 and 32), the top DMR in these analyses being located within HCG9 which is also the top DMR in the unadjusted analysis. Twenty-four DMRs were associated with HbA1c after adjustment for ALMi, the top DMR being located within FERM domain containing 4A (FRMD4A) (Supplementary Table 33).

Table 4 HOMA2-IR-associated DMRs and overlapping genes (Stouffer < 0.05) in skeletal muscle myoblast cells from elderly individualsDmCpGs associated with HOMA2-IR or HbA1c were not related to accelerated epigenetic ageing in skeletal muscle myoblasts

To determine whether the methylation changes associated with HOMA2-IR or HbA1c represented accelerated epigenetic ageing within skeletal muscle myoblasts, epigenetic age was calculated [39]. Epigenetic age as determined by the muscle epigenetic age estimator MEAT was strongly correlated with chronological age (P = 0.002), but there were no associations between accelerated epigenetic age and HOMA2-IR (P = 0.556) or HbA1c (P = 0.552) unadjusted, or after adjustment for BMI, ALMi, grip strength or gait speed (Supplementary Table 34).