Inverse directions of association of higher physical activity and higher insulin resistance with human skeletal muscle cell type abundance and fiber-type-level gene expression

Abstract

To investigate the interplay between physical activity and cardiometabolic traits in human skeletal muscle, we characterized gene expression and chromatin accessibility across skeletal muscle cell types in 263 Finnish individuals from the FUSION Tissue Biopsy Study. We analyzed skeletal muscle single-nucleus RNA-seq data (168,309 nuclei, 23,849 genes), ATAC-seq data (242,069 nuclei, 927,588 peaks), and bulk RNA-seq data (22,309 genes). Lower insulin resistance (HOMA-IR) and higher total physical activity were both associated with higher proportions of Type 1 nuclei and lower proportions of Type 2x nuclei. We identified cell-type-level and tissue-level gene expression-trait and gene set-trait associations for cardiometabolic and physical activity traits, and a smaller proportion of cell-type-level chromatin accessibility-trait associations. Traits typically associated with better health-lower trait values of cardiometabolic traits (BMI, HOMA-IR, normal glucose tolerance vs. type 2 diabetes, 2-hour plasma glucose) and higher physical activity levels (total and vigorous)-were associated with higher expression of energy metabolism genes and lower expression of signaling pathway genes across muscle fiber types, total pseudobulk, and to some extent in bulk tissue. For HOMA-IR and physical activity, these directions of association remained when adjusting for both traits in the same model, indicating apparently independent associations in the same pathways.

Keywords: Type 2 diabetes; chromatin accessibility; gene expression; insulin resistance; physical activity; single nucleus; skeletal muscle; snATAC-seq; snRNA-seq.

Figure 1.. Correlation of trait values and differences in cell type abundance in human skeletal muscle.

(A) Heatmap of Spearman correlations of trait values in 263 individuals. Sample characteristics for each trait can be found in Table S1. Male vs. female and T2D vs. NGT variables were coded as binary 0/1 variables. Physical activity is abbreviated as PHYS ACT. Numbers in bold represent correlations significantly different than 0 tested at a Bonferroni-adjusted significance level of 0.05/28. (B) UMAP projection of 410,378 RNA and ATAC nuclei across 263 individuals with both snATAC and snRNA nuclei. Each point is a nucleus. (C) Cell type proportions for 263 individuals with both snATAC and snRNA nuclei separated by sex and insulin resistance level (HOMA-IR≤2.9 or >2.9). (D) Fold change and 95% confidence intervals comparing the number of snRNA and snATAC nuclei by trait, adjusting for the number of total nuclei (FDR<5% across all tests). Sample sizes for each trait can be found in the methods.

Figure 2.. Gene expression–trait associations across cell types.

(A) Number of significant (FDR<5%) genes for each trait across muscle fiber types, total pseudobulk, and bulk tissue. (B) Heatmap of Spearman correlations of signed (by beta coefficient) −log10 p-value of gene expression–trait associations between cell types for (top left) HOMA-IR and (bottom right) vigorous physical activity. Numbers in bold represent correlations significantly different than 0 tested at a Bonferroni-adjusted significance level of 0.05/132. (C, D, F, G, H) Scatter plots of signed −log10 p-values of gene expression–HOMA-IR associations in (C) Type 2a vs. 2x, (D) Type 1 vs. 2x, (F) Type 1 vs. total pseudobulk, (G) Type 1 vs. bulk tissue, and (H) total pseudobulk vs. bulk tissue. Each point is a gene. (E) Upset plot of the number of significant gene expression–HOMA-IR associations in the muscle fiber types. Number of genes tested and sample sizes for each trait and cell type combination can be found in Tables S15 and S16.

Figure 3.. Gene expression–trait associations across traits.

(A, B) Scatter plots of signed −log10 p-values of gene expression–trait associations for HOMA-IR vs. vigorous physical activity (VIG PHYS ACT) when only adjusting for base covariates in (A) Type 1 and (B) Type 2x. Each point is a gene. (C, D) Scatter plots of signed −log10 p-values of gene expression–trait associations for HOMA-IR vs. vigorous physical activity (VIG PHYS ACT) when jointly adjusting (adj.) for both traits in addition to base covariates in (C) Type 1 and (D) Type 2x. Each point is a gene. (E) Heatmap of Spearman correlations of signed (by beta coefficient) −log10 p-values of gene expression–trait associations between traits in (top left) Type 1 and (bottom right) Type 2x from models only adjusting for base covariates. Numbers in bold represent correlations significantly different than 0 tested at a Bonferroni-adjusted significance level of 0.05/224. Number of genes tested and sample sizes for each trait and cell type combination can be found in Tables S15 and S16.

Figure 4.. Gene set–trait associations across traits.

(A) Number of significant (FDR<5%) gene sets for gene set–trait associations across muscle fiber types, total pseudobulk, and bulk. Gene sets tested are Gene Ontology biological processes. (B, C) Scatter plots of signed −log10 p-values of gene set–trait associations for HOMA-IR vs. vigorous physical activity when only adjusting for base covariates in gene expression models in (B) Type 1 and (C) Type 2x. Each point is a gene set. (D, E) Scatter plots of signed −log10 p-values of gene set–trait associations for HOMA-IR vs. vigorous physical activity when jointly adjusting (adj.) for both traits in gene expression models in (D) Type 1 and (E) Type 2x. Each point is a gene set. (F) Heatmap of Spearman correlations of signed (by beta coefficient) −log10 p-values of gene set–trait associations between traits in (top left) Type 1 and (bottom right) Type 2x from gene expression models only adjusting for base covariates. Numbers in bold represent correlations significantly different than 0 tested at a Bonferroni-adjusted significance level of 0.05/224. Sample sizes for each trait and cell type combination can be found in Table S16.

Figure 5.. Selected significant Gene Ontology (GO) terms across data types and traits.

(A) Gene set results for seven GO terms across all traits and muscle fibers, total pseudobulk, and bulk tissue. A filled-in star indicates that the gene set–trait association is significant (FDR<5%). (B) Histograms of the signed −log10 p-values of gene expression–trait associations of aerobic respiration genes with HOMA-IR and vigorous physical activity in the muscle fiber types. A star indicates that the aerobic respiration–trait association was significant for that trait in that cell type. The number indicates the percentage of genes with positive (red square) or negative (blue circle) gene expression–trait associations for the gene set. (C, D) Scatter plots of the signed −log10 p-values of gene expression–trait associations for HOMA-IR vs. vigorous physical activity for (C) aerobic respiration genes in Type 2x and (D) MAPK cascade genes in Type 1. Arrows indicate the direction of effect of the gene set–trait association from panel A for each trait. Sample sizes for each trait and cell type combination can be found in Table S16.