Integrative analysis of gene expression, DNA methylation, physiological traits, and genetic variation in human skeletal muscle

Abstract

We integrate comeasured gene expression and DNA methylation (DNAme) in 265 human skeletal muscle biopsies from the FUSION study with >7 million genetic variants and eight physiological traits: height, waist, weight, waist-hip ratio, body mass index, fasting serum insulin, fasting plasma glucose, and type 2 diabetes. We find hundreds of genes and DNAme sites associated with fasting insulin, waist, and body mass index, as well as thousands of DNAme sites associated with gene expression (eQTM). We find that controlling for heterogeneity in tissue/muscle fiber type reduces the number of physiological trait associations, and that long-range eQTMs (>1 Mb) are reduced when controlling for tissue/muscle fiber type or latent factors. We map genetic regulators (quantitative trait loci; QTLs) of expression (eQTLs) and DNAme (mQTLs). Using Mendelian randomization (MR) and mediation techniques, we leverage these genetic maps to predict 213 causal relationships between expression and DNAme, approximately two-thirds of which predict methylation to causally influence expression. We use MR to integrate FUSION mQTLs, FUSION eQTLs, and GTEx eQTLs for 48 tissues with genetic associations for 534 diseases and quantitative traits. We identify hundreds of genes and thousands of DNAme sites that may drive the reported disease/quantitative trait genetic associations. We identify 300 gene expression MR associations that are present in both FUSION and GTEx skeletal muscle and that show stronger evidence of MR association in skeletal muscle than other tissues, which may partially reflect differences in power across tissues. As one example, we find that increased RXRA muscle expression may decrease lean tissue mass.

Keywords: DNA methylation; eQTL; gene expression; mQTL; skeletal muscle.

Fig. 1.

Association of FUSION physiological traits with skeletal muscle gene expression and DNAme, controlling for estimated fiber type proportions. Analysis performed with base covariates (sex, age, sample collection site, smoking status, and molecular-trait specific technical covariates), plus tissue type (base+tissue, green bars and points), or base+tissue plus fiber type covariates (base+tissue+fiber, orange bars and points). (A) Percentage of DNAme sites or genes (x axis) associated with each physiological trait (y axis; FDR ≤ 1%). (B) Scatter plot of the number of DNAme sites (y axis) and number of genes (x-axis) associated with each physiological trait adjusting for base+tissue or for base+tissue+fiber covariates (results for a given trait connected with black line).

Fig. 2.

FAM179A-cg09001591 causal analysis. (A, Top) SNV association with cg09001591 (orange) or FAM179A (green). (A, Bottom facet) cg09001591 DNAme site (orange lollipop) and FAM179A gene body (green line). (B) Scatter plot of cg09001591percentage DNAme (x axis) and FAM179A transcripts per million (TPM; y axis). (C) Scatter plot of mQTL effect sizes with SEs (x axis) by eQTL effect sizes with SEs (y axis) for SNVs used in the HEIDI test. The black dashed line is the estimated MR effect based on the top QTL SNV (black triangle). (D) Percent DNAme and gene expression variance explained (y axis) by rs1867944 genotypes. (E) Scatter plot of residual cg09001591 DNAme (adjusted for PEER factors used in eQTM mapping; x axis) and residual FAM179A gene expression (adjusted for PEER factors; y axis). Linear regression line for eQTM association overall (black) and colored by the rs1867944 genotype (TT, green; TC orange; CC, purple; Left). Box plots and linear regression line (additive model) of residual cg09001591 DNAme by rs1867944 genotype (facet M). Box plot and regression line as for M, except with adjustment of residual cg09001591 DNAme by residual FAM179A gene expression (facet M|E). Box plots and linear regression line (additive model) of residual FAM179A gene expression by rs1867944 genotype (facet E). Box plot and regression line as for E except with adjustment of residual FAM179A gene expression by residual cg09001591 DNAme (facet E|M).

Fig. 3.

MR association of UK Biobank GWAS of trunk predicted mass and RXRA-eQTL results in FUSION skeletal muscle and GTEx tissues. (A) MR association (x axis) for UK Biobank trunk predicted mass with the top RXRA-eQTL SNV from FUSION (rs6583658; square) or the top GTEx tissue-specific RXRA-eQTL SNV (triangle) across FUSION skeletal muscle and GTEx tissues (y axis). Power to detect an RXRA MR association (color of tissue name; Methods). ACC, anterior cingulate cortex; GE, gastroesophageal, NAc, nucleus accumbens basal. (B, Top) UK Biobank trunk predicted mass–SNV association (gray points); MR association P values for RXRA (dark blue diamond) and nearby, protein-coding genes (light blue diamond; diamonds drawn at the TSS of the gene). (Middle) FUSION SNV-gene expression association results for RXRA and other nearby genes. (Bottom) Genes in the region. (C) Scatter plot of FUSION RXRA-eQTL effect sizes and SEs (x axis) and trunk predicted mass GWAS effect sizes and SEs (y axis) for SNVs used in the HEIDI test. The black dashed line is the estimated MR effect based on the top QTL SNV (black triangle).