Skip to main content
Download PDF

Skeletal muscle methylome and transcriptome integration reveals profound sex differences related to muscle function and substrate metabolism

Clinical Epigenetics volume 13, Article number: 202 (2021) Cite this article

Abstract

Nearly all human complex traits and diseases exhibit some degree of sex differences, with epigenetics being one of the main contributing factors. Various tissues display sex differences in DNA methylation; however, this has not yet been explored in skeletal muscle, despite skeletal muscle being among the tissues with the most transcriptomic sex differences. For the first time, we investigated the effect of sex on autosomal DNA methylation in human skeletal muscle across three independent cohorts (Gene SMART, FUSION, and GSE38291) using a meta-analysis approach, totalling 369 human muscle samples (222 males and 147 females), and integrated this with known sex-biased transcriptomics. We found 10,240 differentially methylated regions (DMRs) at FDR < 0.005, 94% of which were hypomethylated in males, and gene set enrichment analysis revealed that differentially methylated genes were involved in muscle contraction and substrate metabolism. We then investigated biological factors underlying DNA methylation sex differences and found that circulating hormones were not associated with differential methylation at sex-biased DNA methylation loci; however, these sex-specific loci were enriched for binding sites of hormone-related transcription factors (with top TFs including androgen (AR), estrogen (ESR1), and glucocorticoid (NR3C1) receptors). Fibre type proportions were associated with differential methylation across the genome, as well as across 16% of sex-biased DNA methylation loci (FDR < 0.005). Integration of DNA methylomic results with transcriptomic data from the GTEx database and the FUSION cohort revealed 326 autosomal genes that display sex differences at both the epigenome and transcriptome levels. Importantly, transcriptional sex-biased genes were overrepresented among epigenetic sex-biased genes (p value = 4.6e−13), suggesting differential DNA methylation and gene expression between male and female muscle are functionally linked. Finally, we validated expression of three genes with large effect sizes (FOXO3A, ALDH1A1, and GGT7) in the Gene SMART cohort with qPCR. GGT7, involved in antioxidant metabolism, displays male-biased expression as well as lower methylation in males across the three cohorts. In conclusion, we uncovered 8420 genes that exhibit DNA methylation differences between males and females in human skeletal muscle that may modulate mechanisms controlling muscle metabolism and health.

Introduction

Various diseases, including but not limited to muscular dystrophy, cardiomyopathies, and cardiovascular disease [1], display sex differences in prevalence, onset, progression, or severity. To improve treatment for such diseases, it is crucial to uncover the molecular basis for the sex differences and their consequences on organ function [2]. Sexually differentiated traits and phenotypes stem from a combination of factors, including genetics (gene variants-by-sex interactions [3], XY chromosome complements [4,5,6,7], genomic imprinting [8]), the hormonal milieu [9, 10], and gene regulation [11], with the latter likely contributing the most [1].

Recently, a large-scale study from the Genotype-Tissue Expression (GTEx) consortium unravelled mRNA expression differences between the sexes that are not driven by sex chromosomes, across all tissues. Skeletal muscle was particularly divergent between the sexes, as gene expression profiles in this tissue could predict sex with high specificity and sensitivity [4]. These transcriptomic differences underpin the numerous physiological differences in skeletal muscle between males and females, such as differences in substrate metabolism [12,13,14]. For example, females oxidize more lipids and less carbohydrates and amino acids during endurance exercise, and albeit depending on training status, tend to have a higher proportion of type I (slow-twitch) muscle fibres [15], all of which inherently contribute to enhanced fatigue-resistance in female skeletal muscle [16]. As such, females exhibit higher mRNA and protein levels of lipid oxidation-related genes than males [13]. Interestingly, the top gene set corresponding to sex-biased genes in the GTEx study corresponded to targets of the epigenetic writer polycomb repressive complex 2 (PRC2) and its associated epigenetic mark (H3K27me3) [4]. This suggests that the sex-specific deposition of epigenetic marks may be one of the sources of sex differences in gene expression. Epigenetics is a regulatory system that influences gene expression and is modulated by the genetic sequence and environmental stimuli [17]. DNA methylation is currently the best-characterized epigenetic modification and has been shown to differ between males and females in various tissues, such as pancreatic islets [18], blood [19, 20], and more recently cultured myoblasts and myotubes [21]. While there is ample evidence for transcriptomic sex differences in skeletal muscle [4, 11, 12, 22,23,24], it is unclear whether sex differences exist in the DNA methylome of skeletal muscle tissue, and what factors contribute to these differences. Epigenome-wide association studies (EWAS) are ideal for investigating the impact of sex on genome-wide DNA methylation profiles.

Differences in skeletal muscle fibre-type proportions and circulating levels of sex hormones may contribute to epigenomic and transcriptomic differences between the sexes. Males and females are exposed to differing levels of sex hormones across the lifespan [25], which are primarily ascribed to the reproductive function, but their importance to non-reproductive functions, such as skeletal muscle [26], is becoming more apparent [27, 28]. Sex steroid hormones engage through their specific ligand-receptors [4, 11, 29] and influence transcription and phenotype in a tissue- and sex-specific manner [9, 30,31,32,33]. Not only can these receptors be differentially expressed between sexes [25], they also show sex-biased gene targeting patterns due to intrinsic differences in sex hormone levels [11, 34]. Human muscle consists of three distinct muscle fibre types: type I (slow-twitch oxidative), type IIa (fast-twitch oxidative), and type IIx (fast-twitch) fibres. Males tend to have a higher proportion of fast-twitch type IIa muscle fibres in various muscle groups compared to females [35, 36]. The different fibre types exhibit different methylation patterns [37], as well as different contractile and metabolic properties [38]. Thus, differences in fibre type proportions may underlie sex differences in skeletal muscle DNA methylation.

We performed a large-scale EWAS meta-analysis to explore sex differences in the DNA methylome of human skeletal muscle tissue, using three datasets from our own laboratory and open-access databases (n = 369 individuals; 217 males and 152 females). We established a list of sites (CpG) and regions showing DNA methylation differences between males and females, and explored their genomic context. We then investigated whether transcription factor binding, muscle fibre type distribution, and circulating sex hormone levels explained the observed sex differences in DNA methylation. Next, we integrated the sex-biased DNA methylation with known sex-biased mRNA expression from the GTEx consortium, and inferred the potential downstream effects on skeletal muscle function. Lastly, we validated our findings with transcriptomic data from one cohort used in the meta-analysis and targeted qPCR from another cohort.

Results

Males show profound genome-wide autosomal hypomethylation compared with females in human skeletal muscle

The DNA methylation meta-analysis was conducted on 369 individuals from three datasets (217 males and 152 females). We focused exclusively on the 22 autosomes as sex chromosomes represent a particular case whereby epigenetically driven X-chromosome inactivation takes place exclusively in females. All of the Gene SMART participants were apparently healthy, while the FUSION participants were either healthy or diagnosed with type 2 diabetes mellitus (T2D), and the GSE38291 consisted of monozygotic twins discordant for T2D (Table 1).

Table 1 Characteristics of participants in each data set included in the DNA methylation meta-analysis
Full size table

We found 56,813 differentially methylated positions (DMPs, single CpG sites) between males and females, spread across the 22 autosomes, at a meta-analysis false discovery rate (FDR) < 0.005 (Fig. 1, Additional file 3: Table S2). Ninety-four percent of DMPs were hypomethylated in males compared with females (Fig. 1a). Henceforth throughout the text, “hypo-DMPs” and “hyper-DMPs” refer to positions which are differentially methylated in males compared to females. On average, DNA methylation levels differed by + 2.8% (hyper-DMPs) and − 3.5% (hypo-DMPs) between males and females, with the largest effect sizes reaching + 15.2% and − 35.7%. To determine the degree to which sex explained the variability in DNA methylation at the sex-specific loci, we performed a principal component analysis (PCA) on these loci in each cohort independently, and ran a linear model between sex and the principal component (PC, also referred to as “dimension”) which visually distinguished sex (Fig. 1b; see “Methods”). Participants clustered according to sex when only focusing on the 56,813 DMPs, suggesting that sex explained a substantial amount of variance at the DMPs. In the FUSION (PC1), Gene SMART (PC2), and GSE38291 (PC2) datasets, sex explained 30, 53, and 66% of variability in the relative PC at the sex-DMPs, respectively. Participants did not cluster according to sex when looking at the whole methylome (Additional file 1: Figure S4), which is expected given that ~ 10% of the autosomal genome displayed sex-specific DNA methylation.

Fig. 1
figure 1

Differentially methylated positions (DMPs) with sex in skeletal muscle. a Volcano plot of DNA methylation differences between males compared to females. Each point represents a tested CpG (633,645 in total) and those that appear in colour are DMPs at a meta-analysis false discovery rate < 0.005; red DMPs are hypermethylated in males compared with females; blue DMPs are hypomethylated in males compared with females. The x-axis represents the amount of DNA methylation difference between the sexes, and the y-axis represents statistical significance (higher = more significant). Two DMPs that were present in all three studies and showed the largest effect size are labelled with the respective CpG and boxplots of β-values from each study appear to the right (hyper-DMP) and left (hypo-DMP). b Principal component analysis plots of the methylation values at the DMPs; each point on the graph represents an individual; males denoted in green, females denoted in orange

Full size image

Each data set had a unique study design that required considering different independent variables known to affect DNA methylation, such as age [39] and type 2 diabetes (T2D) [40] in the analysis. We included each of these factors in the relevant linear models, but noted that sex was associated with T2D in the FUSION dataset (Table 1), meaning that male participants from the FUSION cohort more commonly had T2D than females. Therefore, it is possible that the sex-related signal captured in this dataset was partially confounded by T2D. We repeated the meta-analysis excluding T2D participants from the FUSION cohort, but results remained unchanged (Additional file 1: Figure S7).

Fibre type proportions, but not circulating hormones, were associated with differential methylation at loci with sex-specific DNA methylation

Males typically show a greater proportion of type II muscle fibres compared with females [15], and type II fibres exhibit hypomethylation compared to type I fibres [37]. Therefore, we hypothesized that the observed DNA methylation sex differences, specifically the hypomethylation in males, may be a result of differing fibre type distributions between males and females. We first estimated type I fibre proportions in the Gene SMART cohort via immunohistochemistry (Additional file 1: Figure S2D) and the FUSION cohort via RNA-seq (data available for each dataset in Additional file 14: Table S13). In both the Gene SMART and FUSION cohorts, females had higher proportions of type I fibres than males (Additional file 1: Figure S2B/C). We could not directly add fibre type proportions to the linear model as a covariate, since fibre type proportions are not a confounder (i.e. a factor that influences both sex and DNA methylation independently), but may be a direct downstream effect of sex, in turn affecting DNA methylation. Adding fibre type proportions in the model would therefore distort the association between sex and DNA methylation. To overcome this issue, we stratified the cohorts by sex, added fibre type proportions to the model as a covariate and identified DNA methylation patterns associated with fibre type proportions. We then meta-analysed the results to find CpGs robustly associated with fibre type proportions across both cohorts and all sexes (see “Methods”). We identified 16,275 CpGs associated with fibre type proportions (Additional file 1: Figure S2A). When restricting the analysis to the loci exhibiting sex-biased DNA methylation, 8805 (15.5%) of those were associated with fibre type proportions (FDR < 0.005). Effect sizes ranged from − 0.28% to + 0.30% DNA methylation difference per % increase in type I fibre content (Fig. 2a).

Fig. 2
figure 2

Fibre type-related DNA methylation loci across sex-biased DNA methylation loci. a Meta-analysis effect size (x-axis) and meta-analysis significance (y-axis) for the 56,813 tested sex-biased CpGs. Hypomethylated (blue) and hypermethylated (red) point represent differentially methylated positions (DMPs) at false discovery rate (FDR) < 0.005. One hyper- and one hypo-DMPs which showed the largest effect sizes are labelled with the respective CpG, with boxplots of β-values per sex and scatter plots of β-values relative to type I fibre proportion from the Gene SMART (b, c) and FUSION (e, f) cohorts. Females are represented in orange and males in green. d, g Forest plots for the given CpG, showing effect size and confidence intervals for each sex in each study

Full size image

We then aimed to determine whether circulating sex hormone levels underlie the observed DNA methylation sex differences. We analysed estrogen (as estradiol, E2), testosterone (T), free testosterone (Free T), and sex-hormone binding globulin (SHBG) levels using mass spectrometry (Free T derived from calculation) of blood serum in males and females in the Gene SMART cohort. Males and females significantly differed in all four hormone levels (Additional file 13: Table S12). To avoid collinearity with sex, we separated males and females for the association between sex-differential DNA methylation and circulating hormones. We assessed whether each of the four hormone levels was associated with DNA methylation across all of the CpGs and across the sex-DMPs in each sex by adjusting the linear model for a given hormone. In both males and females, circulating free testosterone, testosterone, estrogen, and SHBG levels were not highly associated with DNA methylation (less than five DMPs; FDR < 0.005) of neither all of the CpGs tested nor the sex-DMPs previously identified (Additional file 1: Figure S3). To limit the potentially confounding effect of fluctuating ovarian hormone levels on DNA methylation, female muscle biopsies were collected in the early follicular phase of the menstrual cycle (days 1–7) and blood serum was tested for follicle stimulating hormone (FSH), luteinizing hormone (LH), and progesterone (as well as E2 as previously mentioned). None of the four hormones (when adjusting for each hormone separately or when adjusting for the first two principal components; see ‘methods’), were associated with differential methylation of all the CpGs or the sex-DMPs. This suggests that variations in ovarian hormone levels in the early follicular phase did not confound our results.

Males show profound genome-wide autosomal hypomethylation compared with females in human skeletal muscle

Since the association between DNA methylation and gene expression depends on the genomic context, we inspected the genomic distribution of the DMPs to gain insights into their potential function [41]. We compared the distribution of hyper-, hypo-, and non-DMPs (CpGs included in the analysis which did not present sex-biased DNA methylation) among the various chromatin states in human skeletal muscle using the Roadmap Epigenomics Project [42]. DMPs were not randomly distributed in the chromatin states (χ2 p value < 2.2 × 10−16, Fig. 3a); specifically, hypo-DMPs were enriched in enhancers and depleted in transcription start sites (Additional file 1: Figure S1A), while hyper-DMPs were not enriched in any chromatin states given their scarcity. It should be noted that the Roadmap Epigenomics Project characterizes both male and female skeletal muscle chromatin states regions, and there are 536 regions across 369 unique genes where male and female chromatin states differ (across many tissues including skeletal muscle) [43]. Therefore, we performed the chromatin state enrichment analysis on both the male and female chromatin state annotation in skeletal muscle, which yielded equivalent findings. We noted that DMPs were overrepresented in loci whose chromatin states differ between males and females (38.7% of DMPs vs. 32.4% of non-DMPs are in chromatin states that differ between males and females), which means that the odds of a DMP being located in a sex-differing chromatin state increased by a factor of 1.3 compared with a non-DMP (OR = 0.76, 95% confidence interval = 0.75–0.77, Fisher test p value < 2.2e−16) (Fig. 3b). DMPs were also enriched in CpG island shores and depleted in CpG islands (χ2 p value < 2.2e−16) (Fig. 3c, Additional file 1: Figure S1B). DMPs were over-represented among CpGs previously reported to be differentially methylated between male and female myoblasts (766 CpGs; p value < 0.0001) as well as myotubes (644 CpGs; p value < 0.0001) (see “Methods”).

Fig. 3
figure 3

Genomic context of sex-differentially methylated positions. a Distribution of hyper-/hypo-DMPs and non-DMPs with respect to chromatin states (male skeletal muscle annotation). Blue is hypomethylated in males and red is hypermethylated in males. Red and blue add up to all of the sex-DMPs. Black denotes the rest of the CpG sites from the analysis which are not DMPs. Asterisks represent a greater contribution to the significant relationship between DMP status and chromatin state (Additional file 1: Figure S1A). b Distribution of sex-DMPs and non-DMPs at loci whose chromatin states differ between male and female skeletal muscle. Purple denotes all DMPs (hypo and hyper combined) and black denotes non-DMPs. c Distribution of sex-DMPs and non-DMPs in relation to CpG islands. Asterisks represent a greater contribution to the significant relationship between DMP status and CpG island location (Additional file 1: Figure S1B). d Bee swarm plot showing enrichment of transcription factor binding sites (TFBSs)(−(log10(p value) using Fisher’s exact tests) on the y-axis for differentially methylated positions (DMPs) according to UniBind [45]. The names of the top 10 enriched TFs are denoted by the colour key; brown denotes non-significant TFs. The various data sets for the same TFs are graphed with the corresponding colour

Full size image

Given the role of transcription factors (TFs) in regulating chromatin accessibility and thus effecting downstream gene expression [44], as well as the recent studies identifying sex differences in TF targeting patterns [4, 11], we next tested whether DMPs were enriched for the experimentally validated binding sites (TFBSs) of 268 TFs from 518 different cell and tissue types [45, 46]. The DMPs were overrepresented in the binding sites of 41 TFs (p value < 0.005, Fig. 3d, Additional file 15: Table S14), including hormone-related TFs such as androgen (AR), estrogen (ESR1), and glucocorticoid (NR3C1) receptors.

Gene set enrichment analysis of differentially methylated genes

We next performed Gene set enrichment analysis (GSEA) on the DMGs, as GSEA using epigenomic features may reveal distinct enriched pathways that may not display gene expression differences [11, 43]. We performed GSEA on both the DMRs and DMPs (Fig. 4). GSEA on the DMRs revealed enrichment of several Gene Ontology (GO) terms, one Reactome pathway (“muscle contraction”), but no Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Additional file 11: Table S10) (FDR < 0.005). However, GSEA on the DMPs revealed enrichment across all three databases (Additional file 6: Table S5, Additional file 8: Table S7 and Additional file 10: Table S9). Most of the enriched GO terms are biological process (BP) terms, many of which relate to anatomical structure development as well as many muscle-related processes. Furthermore, DMPs were enriched for GO terms related to substrate metabolism such as lipid, protein, and carbohydrate derivative metabolic processes. Nine-hundred and twenty-five genes of the 1407 genes involved in KEGG metabolic pathways were differentially methylated, representing many aspects of substrate metabolism, although the pathway was only significant when analysing the DMPs.

Fig. 4
figure 4

Gene set enrichment analysis of the differentially methylated genes. a Top enriched Gene Ontology (GO) terms, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and Reactome pathways from GSEA of DMRs and DMPs. Columns represent (1) name (2) number of significant/total number of genes (3) false discovery rate (FDR) of over-representation of term or pathway. Full list of terms and pathways in Additional file 6: Table S5, Additional file 7: Table S6, Additional file 8: Table S7, Additional file 9: Table S8, Additional file 10: Table S9, Additional file 11: Table S10. b Sankey diagram of muscle contraction-related pathways across the three GSEA databases tested and genes within those pathways that were both differentially methylated and expressed (in GTEx and FUSION) between males and females. Numbers next to pathways denote the number of enriched genes in the pathway; numbers next to genes denote the number of pathways (from the ones displayed) that the gene belongs to. Diagram created using SankeyMATIC

Full size image

Genes with sex-biased methylation exhibit sex-biased DNA expression in human skeletal muscle

Differentially methylated genes (DMGs) were determined by identifying the genes annotated to differentially methylated regions (DMRs), as DMRs remove spatial redundancy (CpG sites ~ 500 base pairs apart are typically highly correlated [47]), and may provide more robust and functionally important information than DMPs [48, 49]. We identified 10,240 DMRs (Stouffer, harmonic mean of the individual component FDRs (HMFDR), and Fisher p value < 0.005). These DMRs were annotated to 8420 unique autosomal genes (including non-coding genes) as annotated by Zhou et al. [50] (see “Methods”; Additional file 4: Table S3).

To gain insights into the potential downstream effects of sex-biased DNA methylation on gene expression, we integrated results from the EWAS meta-analysis of sex with genes whose mRNA expression levels are known to differ between males and females. We used version 8 of the Genotype-Tissue Expression (GTEx) database which contains 803 RNA-sequencing profiles in human skeletal muscle (n = 543 males and n = 260 females). There were 2689 sex-differentially expressed genes (DEGs) on the autosomes in skeletal muscle (accessed from GTEx portal on 08/26/2020). Of the 2689 DEGs, 973 (~ 36%) were in common with DMGs from our cohorts (Fig. 5, Additional file 3: Table S2), including the gene Gamma-Glutamyltransferase 7 (GGT7) (Fig. 6). We confirmed an enrichment of DMRs across sex-biased genes (hypergeometric test p value = 4.6 × 10−13), suggesting that the overlap between sex-differentially methylated genes and sex-differentially expressed genes is larger than what would be expected by chance alone. To gain insight on the relationship between DNA methylation and gene expression of sex-biased genes, we assessed the direction of correlation between DMRs that are annotated to either promoter (TssA and TssAFlnk) or enhancer (Enh and EnhG) regions and their given gene expression (Fig. 5c, d). Sixty-two and 59% of DMRs in promoter and enhancer regions, respectively, were inversely correlated with gene expression. (From GTEx transcriptome data, similar results were yielded with the FUSION transcriptome data.) The inverse correlation between DNA methylation at both promoter and enhancer regions with gene expression was more than would be expected to occur by random chance (10,000 random permutations; p value < 0.0001 and p value = 0.0009, respectively; Additional file 1: Figure S5).

Fig. 5
figure 5

Integration of differentially methylated genes and differentially expressed genes. a Venn diagram of the overlap between differentially methylated genes (DMGs; derived from DMRs), differentially expressed genes derived from GTEx (DEGs GTEx), and differentially expressed genes derived from FUSION (DEGs FUSION) between males and females. b Subset of 12 genes with consistently large effect sizes or of biological relevance to skeletal muscle. c Correlation between the effect sizes of DMRs in enhancer regions and the effect sizes of gene expression of the relative annotated gene (for GTEx sex-biased genes). Quadrant percentages indicate the percentage DMRs/DEGs that fall into each quadrant. d Correlation between the effect sizes of DMRs in promoter regions and the effect sizes of gene expression of the relative annotated gene (for GTEx sex-biased genes). Quadrant percentages indicate the percentage DMRs/DEGs that fall into each quadrant

Full size image
Fig. 6
figure 6

adapted from GTEx portal, n = 803). e GGT7 RNAseq expression in the FUSION males and females (FPKM—fragments per kilobase of transcript per million) (n = 274). f GGT7 qPCR expression in a subset of Gene SMART males and females (AU—arbitrary units; 2−∆Ct) (n = 30)

Differential DNA methylation and expression of GGT7 between males and females. a UCSC gene track of GGT7. From top to bottom: base pair scale in black, GENCODE gene tracks transcript variants in blue, GeneHancer regulatory element annotations in light blue, hyper-DMRs tracks in red, hypo-DMRs tracks in blue. b Heatmap of the Gene SMART study (beta values adjusted for all confounders except sex) across the 3 CpGs included in the GGT7 hypo-DMR selected in blue lines and labelled with mean DMR effect size (n = 65). Each row represents an individual; green denotes males and orange denotes females; ordered by similarity to other individuals. Each column corresponds to a CpG in the DMR, ordered by genomic location and corresponding to 5C. Blue denotes hypomethylation; red denotes hypermethylation. c Distribution of DNA methylation (beta values) in males and females, for the three CpGs in the DMR, matching 5B (n = 65). d GGT7 RNAseq expression (TPM—transcripts per million) in males and females of the GTEx (

Full size image

Validation of GTEx sex-biased genes in the cohorts used for methylation analysis

We sought to validate the sex-biased gene expression obtained from GTEx in a subset of the samples used for methylation analysis since the DMGs and DEGs analyses were obtained from different muscle groups. (The DMGs of the current study are from the vastus lateralis, while the GTEx DEGs are from the gastrocnemius.) Although both are skeletal muscle tissue from the leg, there may be differences in muscle phenotypes in differing muscle groups [51]. Analysis of RNA sequencing data from the FUSION cohort revealed 3751 autosomal genes with sex-biased expression (FDR < 0.005). The FDR threshold we chose for the FUSION gene expression data was more stringent than the GTEx local false sign rate threshold (lfsr < 0.05), yet, ~ 34% of the genes which were both DEGs in GTEx and DMGs were also DEGs in the FUSION cohort, totalling 326 genes (hereinto referred to as ‘overlapping genes’) (Fig. 5a). Given that both the GTEx and FUSION cohorts include participants of relatively older ages, we sought to confirm the mRNA levels in the younger cohort in the analysis (the Gene SMART) for three genes that displayed sex differences at both the mRNA and DNA methylation levels (GGT7, FOXO3, and ALDH1A1) (Table 2, Additional file 12: Table S11, Additional file 1: Figure S6).

Table 2 Gene expression and DNA methylation differences between males and females for three genes across the cohorts used in the analysis
Full size table

DNA methylation and gene expression of GGT7, FOXO3 and ALDH1A1 consistently differ between males and females in human skeletal muscle

Three-hundred twenty-six genes exhibited differential methylation in the meta-analysis and differential expression among the GTEx and FUSION cohorts, termed ‘overlapping genes’. Of those genes, we tested three for gene expression levels, GGT7, Forkhead Box O3 (FOXO3), and Aldehyde Dehydrogenase 1 Family Member A1 (ALDH1A1), in the younger cohort included in the DNA methylation analysis (Gene SMART) given the effect that age has on skeletal muscle gene expression [52]. These three genes showed a large effect size in gene expression and DNA methylation, displayed moderate gene expression levels in skeletal muscle relative to other tissues, and/or contained numerous DMPs and DMRs (Table 2). The direction of sex-biased expression was consistent for GGT7 and ALDH1A1 across GTEx, FUSION, and Gene SMART cohorts (GGT7: GTEx lfsr < 2.2e−16, FUSION FDR = 1.3e−45, Gene SMART p value = 0.0003; ALDH1A1: GTEx lfsr < 2.2e−16; FUSION FDR = 2.3e−8, Gene SMART p value = 0.001), while the direction was opposite for FOXO3 (FUSION and GTEx FOXO3 expression lower in males, Gene SMART FOXO3 expression higher in males (GTEx lfsr = 0.01, FUSION FDR = 0.001, Gene SMART p value = 0.0008)). As a specific example of the extent of sex differences across the different layers of analysis, GGT7 displays male-biased expression in skeletal muscle (GTEx lfsr < 2.2e−16; FUSION FDR = 1.3e−45, Gene SMART p value = 0.0003) as well as lower methylation in males at DMPs and DMRs annotated to GGT7 (max DMR: Fisher p value < 0.00−15, max beta value effect size = −28.5%, mean beta value effect size = −20.4%) (Fig. 6).

Discussion

We conducted a large-scale meta-analysis of DNA methylation differences between males and females in skeletal muscle, and integrated them with transcriptomic data. We revealed that males display profound genome-wide hypomethylation compared with females, and that hormone-related TFBSs and muscle fibre type proportions underlie the observed DNA methylation sex differences. We also showed that many sex-biased genes found in GTEx exhibit sex-biased DNA methylation, which was partially confirmed in the FUSION cohort. We then validated the gene expression (qPCR) levels of three genes with large DNA methylation and expression differences between the sexes across cohorts, and confirmed that GGT7 and ALDH1A1 had higher gene expression in males. Finally, we showed that the DMGs are overwhelmingly involved in muscle contraction, as well as substrate metabolism and anatomical structure-related pathways.

In the present study, the overwhelming majority (94%) of the DMPs were hypomethylated in males. Interestingly, global autosomal hypomethylation in males has been observed in various other tissues [53], including blood [54, 55] and pancreatic islets [18]. We aimed to uncover some of the biological mechanisms at the root of these epigenetic sex-differences by assessing the contributions of fibre type proportions and sex hormone levels to the observed DNA methylation sex differences. We hypothesized that differences in fibre type proportions between sexes may partly explain our findings [55,56,57], as studies report that type I fibres are hypermethylated compared with type II fibres [37], and as females tend to have a higher proportion of type I fibres than males [15]. Consistent with this, we observed that females had higher proportions of type I muscle fibres than males and that type I fibre content was mostly associated with DNA hypermethylation. Importantly, 16% of the loci exhibiting sex-biased DNA methylation were also associated with fibre type proportions. This suggests that at those CpGs, differences in DNA methylation between the sexes are due to the inherent sex differences in fibre type proportions. Nonetheless, the vast majority of the loci that exhibit sex-biased DNA methylation (84%; 48,008 CpGs) differ regardless of the sex differences in fibre type proportions. A recent study on the FUSION cohort, adjusted for fibre type proportions, found that it explains a substantial portion of the variability in DNA methylation for many metabolic phenotypes of interest [58]. Skeletal muscle DNA methylation analyses are performed on whole muscle due to the cost and technical limitations of isolating muscle cell types. Differing non-muscle cell types may be present in a muscle biopsy sample, and it is currently unknown how much of the muscle DNA methylation profile may actually be representing other cell types [59,60,61]. Bioinformatics deconvolution methods have not yet been developed for bulk skeletal muscle DNA methylation. Considering that the DNA methylation differences between cell types are large [62], future studies should aim at determining DNA methylation patterns of the different muscle fibre and cell types so that bulk muscle DNA methylation data can be adjusted for the appropriate cell and fibre proportions.

None of the circulating sex hormones were associated with differential methylation across all CpGs, nor across the sex-DMPs in males or females. However, the range of each hormone within each sex may not be large enough to draw out the effect of varying levels of each hormone on the methylome. In the current study, hormone levels were measured from blood, while DNA methylation was measured from skeletal muscle. DNA methylation patterns are highly tissue-specific [63, 64], and sex hormone levels in the circulation are not necessarily correlated with those intramuscularly. Moreover, intramuscular, and not circulating, sex hormone levels may be correlated with muscular function [65, 66]. A recent review emphasizes the importance of measuring intramuscular sex hormone levels when assessing muscle-related properties in females [67].

The enrichment of hormone-related TFBS among the sex-DMPs suggests that lifelong exposure to differing hormone levels significantly contributes to the observed sex differences in skeletal muscle DNA methylation. Although, the alternative is possible, in which DNA methylation differences influence sex hormone levels. In Unibind, ChIP-seq data in skeletal muscle were limited to one TF (CTCF), so the enrichment of TFBSs among sex DMPs may have limited functional significance in skeletal muscle. Nonetheless, many of the TFs that showed strong enrichment in the present study, such as AR [68, 69], ESR1 [69], and SMAD3 [70], are expressed in skeletal muscle and have important roles in muscle phenotype. Of note, multiple testing assumptions cannot be met with Unibind; therefore, no correction can be made, and therefore, false positives may emerge. Two recent studies leveraging the GTEx database identified sex differences in TF targeting patterns across several human tissues, including skeletal muscle, which contribute to sex-biased gene regulatory networks [11] and gene expression [4]. Differences in sex hormone levels between developing males and females are already evident in utero [71], making it challenging to design an experiment in humans that disentangles the effect of long-term hormonal exposure from biological sex, and other related factors, on cell function. Studies have utilized menopausal females [72] and transgender people [73] receiving hormone replacement therapy (HRT) to investigate the influence of long-term exposure to sex hormones on various phenotypes and risk of diseases. For example, HRT for one postmenopausal monozygotic twin and not the other has positive effects on regulation of muscle contraction and myonuclei organization, suggesting that estrogen has direct effects on muscle function [74]. Nevertheless, uncovering the genomic regions that display sex-differential methylation as well as contain hormone-responsive TFBSs, provides insight on which genomic regions, hormones, and TFs are discerning male and female skeletal muscle.

The enrichment for muscle contraction-related pathways among the DMGs across GO, KEGG, and Reactome suggests that the sex-differential DNA methylation has functional relevance in skeletal muscle function. Furthermore, enrichment of substrate metabolism pathways among the DMGs suggests that the observed sex differences in carbohydrate, lipid, and protein metabolism [13, 75] may have a molecular, more specifically an epigenetic, basis. This is corroborated with results from transcriptomic studies, which report that skeletal muscle female-biased genes are enriched for pathways involved in fatty acid metabolism, while male-biased genes are enriched for pathways involved in protein catabolism [22].

Although not well understood, the sex chromosome complement may also influence autosomal DNA methylation patterns. In cultured fibroblasts, the presence of Sex-determining Region Y (SRY) is associated with lower autosomal methylation levels [76,77,78]. Additionally, a higher number of the X chromosomes, in the absence of SRY, leads to increased methylation levels at a specific sex-differentially methylated autosomal region [78]. This could be attributed to allele dosage compensation, a female-specific process that silences one of the X chromosomes in a cell [79, 80]. Approximately one-third of genes ‘escape’ inactivation, remain transcriptionally active in XX cells, [80,81,82], and have been suggested to affect autosomal DNA methylation via their histone marks [78, 83]. Moreover, females with Turner syndrome (partially/fully missing one X) and monosomy X have lower global methylation than XX females, but higher than XY males [84, 85]. The effect of sex chromosome complement on autosomal DNA methylation in skeletal muscle has yet to be explored. Finally, genetic variants (copy number variants (CNVs) and single nucleotide polymorphisms (SNPs)) may affect DNA methylation status at specific loci [86], termed methylation-quantitative trait loci (meQTL). However, with the exception of a female bias for large, rare CNVs, there are no large sex differences in autosomal SNP minor allele frequencies [1].

The relationship between DNA methylation and gene expression is complex; DNA methylation at promoters, enhancers, and 1st exons is generally believed to enhance gene silencing, while DNA methylation at gene bodies can sometimes be associated with increased gene expression [87,88,89,90,91]. Using a permutation test, we showed that DNA methylation differences between the sexes at promoters and enhancers were more often associated with lower gene expression than would be expected by chance alone. DNA methylation differences between the sexes were also particularly prominent in chromatin states that are known differ between males and females. This suggests that DNA methylation differences between males and females reflect alterations in chromatin activity, and differential epigenetic states and expression are likely functionally connected. In line with this, chromatin states that differ between the sexes have been shown to be enriched for sex-biased genes across various tissues, including skeletal muscle [43]. However, it is not yet possible to assess whether the relationship reflects correlation or direct causality. There is still debate around whether epigenomic features drive regulatory processes or are merely a consequence of transcription factor binding [43].

We identified 326 genes with consistent differential skeletal muscle DNA methylation and expression across 1172 individuals altogether (369 individuals from three cohorts for DNA methylation and 1077 individuals from two cohorts for gene expression). Although we found profound global DNA hypomethylation in males, of the overlapping genes there were equivalent numbers of genes over- and under-expressed in males compared with females for both GTEx and FUSION. Indeed, hypermethylation is not always associated with decreased gene expression [92]. The substantial overlap between differentially methylated genes and differentially expressed genes highlights many genes that may be of interest for their roles in muscle-related processes. Furthermore, the significant overlap between the present study and Davegårdh et al. [21] implies that there are specific loci and genes that show DNA methylation differences between the sexes both in muscle cells and in muscle tissue. We focused on three of these genes that displayed a large DNA methylation difference between males and females, are highly expressed in skeletal muscle, or play a role in skeletal muscle function: FOXO3 for its role in ageing, longevity, and regulating the cell cycle [93]; ALDH1A1 for its role in aldehyde oxidation and because sex differences in skeletal muscle mRNA levels have been reported, suggesting that males might be able to metabolize aldehydes (i.e. alcohol) more efficiently than females [12]; and GGT7 for its role in antioxidant activity [94]. Of these three genes (which were validated across GTEx, FUSION, and Gene SMART), FOXO3 and GGT7 have also been reported to exhibit differential methylation between male and female myoblasts as well as myotubes [21]. GGT7 and ALDH1A1 showed consistently higher expression levels in males, while FOXO3 showed opposite sex-biased expression in the young versus the old cohorts. FOXO3 expression was lower in males in the older cohorts (GTEx and FUSION), and higher in males in the younger cohort (Gene SMART). Other studies have shown that males have higher FOXO3 expression in young skeletal muscle [95] and that elderly females have higher skeletal muscle FOXO3 expression than younger females [96]. While FOXO3 skeletal muscle gene expression differs between males and females, it seems that the direction is opposite in young and old individuals, which emphasizes the caution that should be used when interpreting sex differences across a large age range of individuals. Interestingly, FOXO3 was hypomethylated in skeletal muscle with age in a recent study from our group [97]. Regions with differential DNA methylation in the three validated genes were located in enhancers (FOXO3 and ALDH1A1) and regions of strong transcription (GGT7), suggesting that subsequent transcription may be altered. As an example, the promoter, 1st exon, and gene body of GGT7 were hypomethylated in males and males had higher GGT7 expression. GGT7 is highly expressed in skeletal muscle and metabolises glutathione, which is a ubiquitous “master antioxidant” that contributes to cellular homeostasis [98].

In conclusion, we showed that the DNA methylation of hundreds of genes differs between male and female human skeletal muscle. We uncovered important biological factors underlying sex-specific skeletal muscle DNA methylation. Integration of the DNA methylome and transcriptome, as well as gene expression validation, identifies sex-specific genes associated with muscle metabolism and function. Uncovering the molecular basis of sex differences across different tissues will aid in the characterization of muscle phenotypes in health and disease. The effects of other upstream drivers on sex differences in the muscle methylome, such as non-muscle cell type, the XY chromosomes, and genetic variants still need to be explored. Molecular mechanisms that display sex differences in skeletal muscle may help uncover novel targets for therapeutic interventions.

Methods

Datasets

We conducted a meta-analysis of three independent epigenome-wide association studies (EWAS) of sex including the Gene Skeletal Muscle Adaptive Response to Training (SMART) study from our laboratory [99], the Finland–United States Investigation of NIDDM Genetics (FUSION) study from the dbGAP repository (phs000867.v1.p1) [58], and the GSE38291 dataset from the Gene Expression Omnibus (GEO) platform [100]. Detailed participant characteristics, study design, muscle collection, data pre-processing, and data analysis specifications for each study are in Additional file 2: Table S1. Briefly, all studies performed biopsies on the vastus lateralis muscle, all participants were of Caucasian descent (except 1 individual of mixed Caucasian/aboriginal decent), and included either healthy or healthy and T2D individuals aged 18–80 years. The Gene SMART study was approved by the Victoria University Human Ethics Committee (HRE13-223), and written informed consent was obtained from each participant. NIH has approved our request [#96795-2] for the dataset general research use in the FUSION tissue biopsy study.

DNA extraction and methylation method- gene SMART study samples

Genomic DNA was extracted from the samples using the AllPrep DNA/RNA Mini Kit (Qiagen, 80204) following the user manual guidelines. Global DNA methylation profiling was generated with the Infinium MethylationEPIC BeadChip Kit (Queensland University of Technology and Diagenode, Austria). The first set was analysed in 2017 and contained samples (n = 50) from 25 males which were scrambled on the chips to ensure that timepoint and age were not confounded with batch and chip positioning in the array. The second set was analysed in 2019 and contained samples (n = 80) from 20 males and 20 females which were scrambled on the chips to ensure that timepoint, sex, and age were not confounded with batch and chip positioning in the array. The genome-wide DNA methylation pattern was analysed with the Infinium MethylationEPIC BeadChip array.

Bioinformatics and statistical analysis of DNA methylation

Pre-processing

The pre-processing of DNA methylation data was performed according to the bioinformatics pipeline developed for the Bioconductor project [101]. Raw methylation data were pre-processed, filtered and normalized across samples. Probes that had a detection p value of > 0.01, located on X and Y chromosomes or cross-hybridizing, or related to a SNP frequent in European populations, were removed. It is important to note that the list of cross-hybridizing probes was supplied manually [102] as the list supplied to the ChAMP package was outdated. Specifically, there are thousands of probes in the Illumina microarrays that cross-hybridize with the X-chromosome and may lead to false discovery of autosomal sex-associated DNA methylation [103]. The BMIQ algorithm was used to correct for the Infinium type I and type II probe bias. β-values were corrected for both batch and position in the batch using ComBat [104].

Statistical analysis

We adjusted each EWAS for bias and inflation using the empirical null distribution as implemented in BACON [105]. Inflation and bias in EWAS are caused by unmeasured technical and biological confounding, such as population substructure, batch effects, and cellular heterogeneity [106]. The inflation factor is higher when the expected number of true associations is high (as it is for age); it is also greater for studies with higher statistical power [105]. The results were consistent with the inflation factors and biases reported in an EWAS of age in blood [105]. Smaller cohorts were weighed less in the meta-analysis than bigger cohorts as the pooled estimates are based on standard error (which takes sample size into account). Results from the independent EWAS were combined using an inverse variance weighted meta-analysis with METAL [107]. We used METAL since it does not require all DNA methylation datasets to include every CpG site on the HumanMethylation arrays. Given that the FUSION and Gene SMART datasets were analysed with the EPIC array, while the GSE38291 dataset was analysed with the 27 K array, there are a substantial number of probes with data available only from the FUSION and Gene SMART studies. For robustness, we only included CpGs present in at least two of the three cohorts after pre-processing (633,645 CpGs). We used a fixed effects (as opposed to random effects) meta-analysis, assuming one true effect size of sex on DNA methylation, which is shared by all the included studies. Nevertheless, Cochran's Q-test for heterogeneity was performed to test whether effect sizes were homogeneous between studies (a heterogeneity index (I2) > 50% reflects heterogeneity between studies). We noted one sample from the FUSION cohort had a discordant sex prediction, using the getSex function from minfi [108], and reran analyses without that sample but obtained identical pathways and differentially methylated genes.

To identify DMPs, we used linear models as implemented in the limma package in R [109], using the participants’ ID as a blocking variable to account for the repeated measures design (for twin (GSE38291) and duplicate samples (Gene SMART), using DuplicateCorrelation). The main sources of variability in methylation varied depending on the cohort and were included as independent variables in the linear model accordingly. For the Gene SMART study, the linear model was of the form:

$${\text{DNAm}}\sim\,{\text{sex}}*{\text{timepoint}} + {\text{age}} + {\text{analysis}}\,{\text{set}}$$

Timepoint: before and after four weeks of high-intensity interval training

Analysis set: one set analysed in 2017 and one in 2019

For the FUSION study, the linear model was of the form:

$${\text{DNAm}}\sim\,{\text{sex}} + {\text{BMI}} + {\text{smoking}}\,{\text{status}} + {\text{oral}}\,{\text{glucose}}\,{\text{tolerance}}\,{\text{test}}\,{\text{status}}$$

For the GSE38291 study, the linear model was of the form:

$${\text{DNAm}}\sim\,{\text{sex}} + {\text{age}} + {\text{diabetes}}\,{\text{status}}$$

All results were adjusted for multiple testing using the Benjamini and Hochberg correction [110], and all CpGs showing an FDR < 0.005 were considered significant [111]. DMRs were identified using the DMRcate package with default settings [112]. DMRs with Stouffer, Fisher, and harmonic mean of the individual component FDRs (HMFDR) statistics < 0.005 were deemed significant. Effect sizes are reported as mean differences in DNA methylation (%) between the sexes. PCAs were analysed using the package factoMiner, and were performed separately per cohort given that technical variability substantially contributes to variability in DNA methylation. The percentage reported in which sex explains the variability in DNA methylation at DMPs was determined by extracting the adjusted R2 from the linear model: PC ~ sex; the PC which visually distinguished males and females in the PCA in each cohort was used in the linear model. To assess whether fibre type proportions differ between males and females in the Gene SMART and FUSION cohorts, we used a beta regression model [113] using the betareg package in R. We then included type I fibre ratio as a covariate in the linear models. To investigate whether circulating hormones is associated with DNA methylation of the sex-DMPs, we included each sex hormone as a covariate in a linear model, separately. One male from the Gene SMART cohort had missing circulating hormone values; and one two females from the Gene SMART cohort had missing type I fibre proportions. For these three individuals, missing values were imputed with the mice package in R [114].

To annotate the CpGs to the genome (see “Availability of Data and Materials”), we integrated a comprehensive annotation of Illumina HumanMethylation arrays [50] with chromatin states from the Roadmap Epigenomics Project [42] and the latest GeneHancer information [115]. DMGs were determined by the annotations of the DMRs. For chromatin state enrichment, DMPs that were annotated to two differing chromatin states were removed for simplicity and because there were very few such DMPs. GSEA on KEGG and GO databases was performed on DMRs and DMPs using the goregion and gometh (gsameth for Reactome) functions in the missMethyl R package [116] 117. All analyses were performed using the R software version 4.0.3.

Enrichment of TFBSs among the identified DMPs was performed using the enrichment analysis tool in http://unibind.uio.no/ which utilizes the runLOLA function of the R package LOLA [46]. The analysis type selected was enrichment with background; and background provided included all tested CpGs. Analysis using CpGs in DMRs versus DMPs yielded similar results. The outputted p value is reported; without applying an FDR correction since many TFBS datasets are not independent so the assumptions behind the FDR correction cannot be met. TFs from all tissues were included as only one TF, to date, has been validated among skeletal muscle tissues or cells.

To assess whether sex-DMPs were over-represented among CpGs which have been reported to be differentially methylated in male and female myoblasts and myotubes [21], we first overlapped the sex-DMPs with CpGs with sex-differential DNA methylation in myoblasts. We then randomly sampled the same number of CpGs which were differentially methylated in myoblasts from the total CpGs in the 450 K array, as this was the Illumina array used for the study, and overlapped those randomly sampled CpGs with the sex-DMPs. We performed this random sampling 10,000 times to obtain a p value signifying whether the amount of overlap between the two studies could have occurred by chance. We then repeated this same method for the myotubes.

Fibre types: meta-analysis and derivation from immunohistochemistry and RNA-seq

Myosin heavy chain is currently the best available marker for fibre typing [118]. Gene SMART muscle sections were frozen in optimum-cutting temperature (OCT) medium by holding the sample with OCT in liquid-nitrogen cooled isopentane until frozen. Samples were stored in − 80 °C until they were sectioned at 8 µM with a cryostat. The IHC protocol was performed as is described elsewhere [119]. Briefly, sections were blocked in 4% goat serum (Life Technologies). Primary antibodies BA-F8 (MHCI), BF-35 (MHCIIA), and 6H1 (MHCIIX) were purchased from DSHB, Iowa. Secondary antibodies goat anti-mouse IgG2b 350, goat anti-mouse IgG1 488, and goat anti-mouse IgM 555 were purchased from Invitrogen. Some samples were fixed in paraformaldehyde for other analyses, and for those samples, an antigen-retrieval protocol consisting of a 10-min incubation at 50 °C of Proteinase K diluted in MilliQ (1:1000) and subsequent 1-min washes was performed before the IHC. Imaging was performed on the Olympus BX51.

To determine type I fibre proportions in the FUSION cohort, we followed the validated method as performed by the original study on the FUSION cohort [58]. Briefly, we derived type I fibre proportions from the RNA-seq expression data (TPMs) for type I (MYH7), type IIA (MYH2), and type IIX (MYH1). We calculated the ratio of MYH7 out of the total. We then included this ratio in the + FT linear model.

To determine whether the inherent sex differences in fibre type proportions underlie the sex differences in DNA methylation, we separated the males and females of the Gene SMART and FUSION cohorts and performed a meta-analysis on the four groups (FUSION females, FUSION males, Gene SMART females, Gene SMART males). Given that females displayed significantly higher type I fibre proportions than males in both cohorts, we could not simply include type I fibre content in a linear model performed on a mixed sex cohort as two issues would arise: (1) collinearity of fibre type with sex, (2) differences in fibre type proportions may be a downstream effect of sex. Dividing the cohorts by sex, conducting a meta-analysis, and selecting the sex-biased DMPs and performing an FDR adjustment among those cites allowed us to address whether fibre type proportion is associated with DNA methylation at sex-biased DNA methylation loci. The fibre type meta-analysis was performed with the same methodology of the sex meta-analysis as described above.

Controlling for the female menstrual cycle

Various contraceptives have different dosage, administration patterns, and different hormone combinations causing variability in metabolism and gene expression [120]; therefore, only females not taking any form of hormonal contraceptives were recruited for the Gene SMART study. Furthermore, to minimise the effect of fluctuating hormone levels, females were required to have a regular menstrual cycle (27–35 days), and all samples were aimed to be collected during the early follicular phase (day 1-day 8 of cycle), with few exceptions due to logistics. Estrogen, progesterone, follicle stimulating hormone (FSH), and luteinizing hormone (LH) were measured in blood serum. Given the intricate fluctuations of the ovarian hormones, these four hormones were combined into a principal component analysis and the first two PCs, which explained the majority of the variability, were each added into the linear model.

Blood serum hormones

The hormone assays were completed in the accredited pathology laboratory at Monash Health, Australia. Estradiol (E2) and Progesterone assays are competitive binding immunoenzymatic assays performed on the Unicel DXI 800 system (Beckman Coulter). FSH assay is based on Microparticle Enzyme Immunoassay (MEIA) and is carried out on the Unicel DXI 800 system (Beckman Coulter). The LH and sex-hormone binding globulin (SHBG) assays were performed using a sequential two-step immunoenzymatic (“sandwich”) assay carried out on a Unicel DXI 800 (Beckman Coulter). Testosterone was measured using the HPLC–tandem mass spectrometry method using a liquid sample extraction (AB Sciex Triple Quad 5500 liquid chromatography–tandem mass spectrometry). Free androgen index was calculated as (total testosterone × 100)/SHBG. Free testosterone was calculated by the Södergard free testosterone calculation (36).

Integration of DNA Methylation and Gene Expression

The Genotype-Tissue Expression (GTEx) Project sex-biased data were downloaded from the GTEx Portal on 08/26/2020 and filtered for skeletal muscle samples. The enrichment of DMG for GTEx DEGs was done by supplying the list of sex-biased genes to the gsameth function in the missMethyl R package [116, 117], which performs a hypergeometric test and outputs a p value and which genes are significant within the supplied list, taking into account biases due to the number of CpG sites per gene and the number of genes per probe on the EPIC array. Therefore, caution should be taken when interpreting the number of DMPs reported per DMG. The analysis for direction of correlation between DNA methylation and gene expression was performed by randomly shuffling DNA methylation effect sizes and performing 10,000 permutations to assess how often a negative correlation occurs. This analysis was performed for both GTEx and FUSION transcriptome data and yielded similar results; data presented reflect results from the integration of differential methylation with differential GTEx expression. Significance reported for GTEx sex-biased genes is represented as the local false sign rate (lfsr) which is analogous to FDR [121]. GTEx effect sizes are represented as mash posterior effect sizes [121], in which positive values indicate male-biased genes and negative values indicate female-biased genes. FUSION and Gene SMART gene expression significance statistics are represented as FDR and p value, respectively, and effect sizes as fold changes for both cohorts.

Validation of top genes with qPCR

Skeletal muscle previously stored at − 80 °C was lysed with the RLT buffer Plus buffer (Qiagen) and beta-mercaptoethanol using the TissueLyser II (Qiagen, Australia). DNA was extracted using the AllPrep DNA/RNA Mini Kit following the manufacturer guidelines (Qiagen, Australia). RNA yield and purity were assessed using the spectrophotometer (NanoDrop One, Thermo Fisher). RNA was reverse transcribed to cDNA using a commercially available iScript Reverse Transcriptase supermix (cat #1708841) and C1000 Touch Thermal Cycler (Bio-Rad, Hercules, CA, USA). Complementary DNA samples were stored at − 20 °C until further analysis. Quantitative real-time PCR was performed using SYBR Green Supermix (Bio-Rad, Hercules, CA) and gene-specific primers (listed in Additional file 12: Table S11). Primers were either adapted from existing literature or designed using Primer-BLAST (http://www.ncbi.nlm.nih.gov/tools/primer-blast/) to include all splice variants, and were purchased from Integrated DNA Technologies. Ten microliter reactions comprised of SYBR, and optimized concentrations of forward and reverse primers (Additional file 12: Table S11 for primer conditions), nuclease-free water and 8 ng of cDNA were run in triplicate using an automated pipetting system (epMotion M5073, Eppendorf, Hamburg, Germany), with no-template negative controls on a 384-well plate in a thermo-cycler (QuantStudio 12 K Flex System, Thermo Fisher Scientific, Australia). Gene expression was normalized to the geometric mean expression of the two most stable housekeeping genes, as determined by Ref finder, TATAA-box binding protein (TPB), and 18 s rRNA, which did not differ between sexes (Additional file 12: Table S11). Data are presented as the fold change in males compared to females, using 2−∆∆CT.

Significance

The importance of uncovering biological sex differences and their translation to physiology has become increasingly evident. Using a large-scale meta-analysis of three cohorts, we perform the first comparison of genome-wide skeletal muscle DNA methylation between males and females, and identify thousands of genes that display sex-differential methylation. We then explore intrinsic biological factors that may be underlying the DNA methylation sex differences, such as fibre type proportions and sex hormones. Leveraging the GTEx database, we identify hundreds of genes with both sex-differential expression and DNA methylation in skeletal muscle. We further confirm the sex-biased genes with gene expression data from two cohorts included in the methylation meta-analysis. Our study integrates genome-wide sex-biased DNA methylation and expression in skeletal muscle, shedding light on distinct sex differences in skeletal muscle.

Availability of data and materials

The dataset generated and analysed during the current study is available in the GEO repository, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE171140. The R codes can be found on Github for constructing the annotation table (https://github.com/sarah-voisin/Annotation-of-Illumina-HumanMethylation-arrays) as well as the main analysis (https://github.com/shanie-landen/sex-specific-skeletal-muscle-DNA-methylation.git).

References

  1. Khramtsova EA, Davis LK, Stranger BE. The role of sex in the genomics of human complex traits. Nat Rev Genet. 2018;1:173–90.

    Google Scholar 

  2. Mamlouk GM, et al. Sex bias and omission in neuroscience research is influenced by research model and journal, but not reported NIH funding. Front Neuroendocrinol. 2020;57:100835.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Rawlik K, Canela-Xandri O, Tenesa A. Evidence for sex-specific genetic architectures across a spectrum of human complex traits. Genome Biol. 2016;17(1):166.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  4. Oliva M, et al. The impact of sex on gene expression across human tissues. Science. 2020;369(6509):eaba3066.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Zore T, Palafox M, Reue K. Sex differences in obesity, lipid metabolism, and inflammation—A role for the sex chromosomes? Mol Metab. 2018;15:35–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Arnold AP. Y chromosome’s roles in sex differences in disease. Proc Natl Acad Sci. 2017;114(15):3787–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Arnold AP, Chen X, Itoh Y. What a difference an X or Y makes: sex chromosomes, gene dose, and epigenetics in sexual differentiation. In: Regitz-Zagrosek V, editor. Sex and gender differences in pharmacology. Berlin: Springer; 2013. p. 67–88.

    Chapter  Google Scholar 

  8. Golden LC, et al. Parent-of-origin differences in DNA methylation of X chromosome genes in T lymphocytes. Proc Natl Acad Sci. 2019;116(52):26779–87.

    Article  CAS  PubMed Central  Google Scholar 

  9. Pihlajamaa P, et al. Tissue-specific pioneer factors associate with androgen receptor cistromes and transcription programs. EMBO J. 2014;33(4):312–26.

    CAS  PubMed  PubMed Central  Google Scholar 

  10. Varlamov O, Bethea CL, Roberts CT Jr. Sex-specific differences in lipid and glucose metabolism. Front Endocrinol. 2015;5:241.

    Article  Google Scholar 

  11. Lopes-Ramos CM, et al. Sex differences in gene expression and regulatory networks across 29 human tissues. Cell Rep. 2020;31(12):107795.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Maher AC, et al. Sex differences in global mRNA content of human skeletal muscle. PLoS ONE. 2009;4(7):e6335.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  13. Tarnopolsky M. Sex differences in exercise metabolism and the role of 17-beta estradiol. Med Sci Sports Exerc. 2008;40(4):648–54.

    Article  CAS  PubMed  Google Scholar 

  14. Landen S, et al. Genetic and epigenetic sex-specific adaptations to endurance exercise. Epigenetics. 2019;14(6):523–35.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Haizlip KM, Harrison BC, Leinwand LA. Sex-based differences in skeletal muscle kinetics and fiber-type composition. Physiology (Bethesda). 2015;30(1):30–9.

    CAS  PubMed Central  Google Scholar 

  16. Hunter SK. Sex differences in human fatigability: mechanisms and insight to physiological responses. Acta Physiol (Oxf). 2014;210(4):768–89.

    Article  CAS  PubMed Central  Google Scholar 

  17. Moore LD, Le T, Fan G. DNA methylation and its basic function. Neuropsychopharmacology. 2013;38(1):23–38.

    Article  CAS  PubMed  Google Scholar 

  18. Hall E, et al. Sex differences in the genome-wide DNA methylation pattern and impact on gene expression, microRNA levels and insulin secretion in human pancreatic islets. Genome Biol. 2014;15(12):522.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  19. Singmann P, et al. Characterization of whole-genome autosomal differences of DNA methylation between men and women. Epigenet Chromatin. 2015;8(1):1–13.

    Article  CAS  Google Scholar 

  20. Liu J, et al. A study of the influence of sex on genome wide methylation. PLoS ONE. 2010;5(4):e10028.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  21. Davegårdh C, et al. Sex influences DNA methylation and gene expression in human skeletal muscle myoblasts and myotubes. Stem Cell Res Ther. 2019;10(1):26.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  22. Lindholm ME, et al. The human skeletal muscle transcriptome: sex differences, alternative splicing, and tissue homogeneity assessed with RNA sequencing. Faseb J. 2014;28(10):4571–81.

    Article  CAS  PubMed  Google Scholar 

  23. Welle S, Tawil R, Thornton CA. Sex-related differences in gene expression in human skeletal muscle. PLoS ONE. 2008;3(1):e1385.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Gershoni M, Pietrokovski S. The landscape of sex-differential transcriptome and its consequent selection in human adults. BMC Biol. 2017;15(1):7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. Vingren JL, et al. Effect of resistance exercise on muscle steroid receptor protein content in strength-trained men and women. Steroids. 2009;74(13–14):1033–9.

    Article  CAS  PubMed  Google Scholar 

  26. Alexander SE, Pollock AC, Lamon S. The effect of sex hormones on skeletal muscle adaptation in females. Eur J Sport Sci. 2021. https://doi.org/10.1080/17461391.2021.1921854.

    Article  PubMed  Google Scholar 

  27. Pihlajamaa P, Sahu B, Jänne OA. Determinants of receptor- and tissue-specific actions in androgen signaling. Endocr Rev. 2015;36(4):357–84.

    Article  CAS  PubMed  Google Scholar 

  28. Varlamov O, Bethea CL, Roberts CT Jr. Sex-specific differences in lipid and glucose metabolism. Front Endocrinol (Lausanne). 2014;5:241.

    Google Scholar 

  29. Fornes O, et al. JASPAR 2020: update of the open-access database of transcription factor binding profiles. Nucleic Acids Res. 2020;48(D1):D87-d92.

    CAS  PubMed  Google Scholar 

  30. Rajan P, et al. Identification of novel androgen-regulated pathways and mRNA isoforms through genome-wide exon-specific profiling of the LNCaP transcriptome. PLoS ONE. 2011;6(12):e29088.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Wu Y, et al. Identification of androgen response elements in the insulin-like growth factor I upstream promoter. Endocrinology. 2007;148(6):2984–93.

    Article  CAS  PubMed  Google Scholar 

  32. Hevener AL, et al. The impact of ERα action on muscle metabolism and insulin sensitivity - Strong enough for a man, made for a woman. Mol Metab. 2018;15:20–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Wiik A, et al. Expression of both oestrogen receptor alpha and beta in human skeletal muscle tissue. Histochem Cell Biol. 2009;131(2):181–9.

    Article  CAS  PubMed  Google Scholar 

  34. Mayne BT, et al. Large scale gene expression meta-analysis reveals tissue-specific, sex-biased gene expression in humans. Front Genet. 2016;7:183.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  35. Maughan R, Watson JS, Weir J. Strength and cross-sectional area of human skeletal muscle. J Physiol. 1983;338(1):37–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Carter S, et al. Changes in skeletal muscle in males and females following endurance training. Can J Physiol Pharmacol. 2001;79(5):386–92.

    Article  CAS  PubMed  Google Scholar 

  37. Begue G, et al. DNA methylation assessment from human slow-and fast-twitch skeletal muscle fibers. J Appl Physiol. 2017;122(4):952–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Zierath JR, Hawley JA. Skeletal muscle fiber type: influence on contractile and metabolic properties. PLoS Biol. 2004;2(10):e348.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  39. Voisin S, et al. Meta-analysis of genome-wide DNA methylation and integrative OMICs in human skeletal muscle. bioRxiv, 2020.

  40. Davegårdh C, et al. DNA methylation in the pathogenesis of type 2 diabetes in humans. Mol Metab. 2018;14:12–25.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  41. Jones PA, Takai D. The role of DNA methylation in mammalian epigenetics. Science. 2001;293(5532):1068–70.

    Article  CAS  PubMed  Google Scholar 

  42. Kundaje A, et al. Integrative analysis of 111 reference human epigenomes. Nature. 2015;518(7539):317–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Yen A, Kellis M. Systematic chromatin state comparison of epigenomes associated with diverse properties including sex and tissue type. Nat Commun. 2015;6(1):1–13.

    Article  CAS  Google Scholar 

  44. Voss TC, Hager GL. Dynamic regulation of transcriptional states by chromatin and transcription factors. Nat Rev Genet. 2014;15(2):69–81.

    Article  CAS  PubMed  Google Scholar 

  45. Gheorghe M, et al. A map of direct TF–DNA interactions in the human genome. Nucl Acids Res. 2019;47(4):e21–e21.

    Article  PubMed  CAS  Google Scholar 

  46. Puig RR, et al. UniBind: maps of high-confidence direct TF-DNA interactions across nine species. bioRxiv, 2020.

  47. Guo S, et al. Identification of methylation haplotype blocks aids in deconvolution of heterogeneous tissue samples and tumor tissue-of-origin mapping from plasma DNA. Nat Genet. 2017;49(4):635–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. VanderKraats ND, et al. Discovering high-resolution patterns of differential DNA methylation that correlate with gene expression changes. Nucleic Acids Res. 2013;41(14):6816–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Schlosberg CE, VanderKraats ND, Edwards JR. Modeling complex patterns of differential DNA methylation that associate with gene expression changes. Nucleic Acids Res. 2017;45(9):5100–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Zhou W, Laird PW, Shen H. Comprehensive characterization, annotation and innovative use of Infinium DNA methylation BeadChip probes. Nucleic Acids Res. 2017;45(4):e22–e22.

    PubMed  Google Scholar 

  51. Miller AEJ, et al. Gender differences in strength and muscle fiber characteristics. Eur J Appl Physiol. 1993;66(3):254–62.

    Article  CAS  Google Scholar 

  52. Su J, et al. A novel atlas of gene expression in human skeletal muscle reveals molecular changes associated with aging. Skeletal Muscle. 2015;5(1):1–12.

    Article  CAS  Google Scholar 

  53. McCarthy NS, et al. Meta-analysis of human methylation data for evidence of sex-specific autosomal patterns. BMC Genom. 2014;15(1):981.

    Article  CAS  Google Scholar 

  54. Yousefi P, et al. Sex differences in DNA methylation assessed by 450 K BeadChip in newborns. BMC Genom. 2015;16(1):911.

    Article  CAS  Google Scholar 

  55. Inoshita M, et al. Sex differences of leukocytes DNA methylation adjusted for estimated cellular proportions. Biol Sex Differ. 2015;6(1):1–7.

    Article  CAS  Google Scholar 

  56. Bauer M. Cell-type-specific disturbance of DNA methylation pattern: a chance to get more benefit from and to minimize cohorts for epigenome-wide association studies. Int J Epidemiol. 2018;47(3):917–27.

    Article  PubMed  Google Scholar 

  57. Suzuki M, Bird A, Suzuki MM, Bird A. DNA methylation landscapes: provocative insights from epigenomics. Nat Rev Genet. 2008;9:465–76.

    Article  CAS  PubMed  Google Scholar 

  58. Taylor DL, et al. Integrative analysis of gene expression, DNA methylation, physiological traits, and genetic variation in human skeletal muscle. Proc Natl Acad Sci. 2019;116(22):10883–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Rubenstein AB, et al. Single-cell transcriptional profiles in human skeletal muscle. Sci Rep. 2020;10(1):1–15.

    Article  CAS  Google Scholar 

  60. De Micheli AJ, et al. A reference single-cell transcriptomic atlas of human skeletal muscle tissue reveals bifurcated muscle stem cell populations. Skeletal muscle. 2020;10(1):1–13.

    Article  CAS  Google Scholar 

  61. Domenig SA, Palmer AS, Bar-Nur O. Stem cell-based and tissue engineering approaches for skeletal muscle repair. In: Eberli D, Lee SJ, Traweger A, editors. Organ tissue engineering. Berlin: Springer; 2020. p. 1–62.

    Google Scholar 

  62. Suzuki MM, Bird A. DNA methylation landscapes: provocative insights from epigenomics. Nat Rev Genet. 2008;9(6):465–76.

    Article  CAS  PubMed  Google Scholar 

  63. Pai AA, et al. A genome-wide study of DNA methylation patterns and gene expression levels in multiple human and chimpanzee tissues. PLoS Genet. 2011;7(2):e1001316.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Zhou J, et al. Tissue-specific DNA methylation is conserved across human, mouse, and rat, and driven by primary sequence conservation. BMC Genom. 2017;18(1):1–17.

    Article  Google Scholar 

  65. Pöllänen E, et al. Differential influence of peripheral and systemic sex steroids on skeletal muscle quality in pre-and postmenopausal women. Aging Cell. 2011;10(4):650–60.

    Article  PubMed  CAS  Google Scholar 

  66. Pöllänen E, et al. Intramuscular sex steroid hormones are associated with skeletal muscle strength and power in women with different hormonal status. Aging Cell. 2015;14(2):236–48.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  67. Alexander SE, Pollock AC, Lamon S. The effect of sex hormones on skeletal muscle adaptation in females: influence of sex hormones on female muscle physiology. Eur J Sport Sci. 2021. doi: https://doi.org/10.1080/17461391.2021.1921854.

  68. Sinha-Hikim I, et al. Androgen receptor in human skeletal muscle and cultured muscle satellite cells: up-regulation by androgen treatment. J Clin Endocrinol Metab. 2004;89(10):5245–55.

    Article  CAS  PubMed  Google Scholar 

  69. Ghanim H, et al. Diminished androgen and estrogen receptors and aromatase levels in hypogonadal diabetic men: reversal with testosterone. Eur J Endocrinol. 2018;178(3):277–83.

    Article  CAS  PubMed  Google Scholar 

  70. Amar D, et al. Time trajectories in the transcriptomic response to exercise-a meta-analysis. Nat Commun. 2021;12(1):1–12.

    Article  CAS  Google Scholar 

  71. Meakin AS, et al. Let’s talk about placental sex, baby: understanding mechanisms that drive female-and male-specific fetal growth and developmental outcomes. Int J Mol Sci. 2021;22(12):6386.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Skinner BD, et al. A systematic review and meta-analysis examining whether changing ovarian sex steroid hormone levels influence cerebrovascular function. Front Physiol. 2021;12:687591.

    Article  PubMed  PubMed Central  Google Scholar 

  73. Kotamarti VS, et al. Risk for venous thromboembolism in transgender patients undergoing cross-sex hormone treatment: a systematic review. J Sex Med. 2021;18:1280–91.

    Article  CAS  PubMed  Google Scholar 

  74. Qaisar R, et al. Hormone replacement therapy improves contractile function and myonuclear organization of single muscle fibres from postmenopausal monozygotic female twin pairs. J Physiol. 2013;591(9):2333–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Lundsgaard A-M, Kiens B. Gender differences in skeletal muscle substrate metabolism–molecular mechanisms and insulin sensitivity. Front Endocrinol. 2014;5:195.

    Article  Google Scholar 

  76. Wijchers PJ, et al. Sexual dimorphism in mammalian autosomal gene regulation is determined not only by sry but by sex chromosome complement as Well. Dev Cell. 2010;19(3):477–84.

    Article  CAS  PubMed  Google Scholar 

  77. De Vries GJ. Minireview: sex differences in adult and developing brains: compensation, compensation. Compens Endocrinol. 2004;145(3):1063–8.

    Article  CAS  Google Scholar 

  78. Ho B, et al. X chromosome dosage and presence of SRY shape sex-specific differences in DNA methylation at an autosomal region in human cells. Biol Sex Differ. 2018;9(1):10.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  79. Penny GD, et al. Requirement for Xist in X chromosome inactivation. Nature. 1996;379(6561):131–7.

    Article  CAS  PubMed  Google Scholar 

  80. Carrel L, Willard HF. X-inactivation profile reveals extensive variability in X-linked gene expression in females. Nature. 2005;434(7031):400–4.

    Article  CAS  PubMed  Google Scholar 

  81. Sidorenko J, et al. The effect of X-linked dosage compensation on complex trait variation. Nat Commun. 2019;10(1):3009.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  82. Tukiainen T, et al. Landscape of X chromosome inactivation across human tissues. Nature. 2017;550(7675):244–8.

    Article  PubMed  PubMed Central  Google Scholar 

  83. Grafodatskaya D, et al. Multilocus loss of DNA methylation in individuals with mutations in the histone H3 lysine 4 demethylase KDM5C. BMC Med Genomics. 2013;6:1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Trolle C, et al. Widespread DNA hypomethylation and differential gene expression in turner syndrome. Sci Rep. 2016;6(1):1–14.

    Article  CAS  Google Scholar 

  85. Sharma A, et al. DNA methylation signature in peripheral blood reveals distinct characteristics of human X chromosome numerical aberrations. Clin Epigenetics. 2015;7(1):1–15.

    Article  CAS  Google Scholar 

  86. Villicaña S, Bell JT. Genetic impacts on DNA methylation: research findings and future perspectives. Genome Biol. 2021;22(1):1–35.

    Article  CAS  Google Scholar 

  87. Rauch TA, Wu X, Zhong X, Riggs AD, Pfeifer GP. A human B cell methylome at 100-base pair resolution. Proc Natl Acad Sci USA. 2009;106(3):671-8. https://doi.org/10.1073/pnas.0812399106.

  88. Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, Nery JR, Lee L, Ye Z, Ngo QM, Edsall L, Antosiewicz-Bourget J, Stewart R, Ruotti V, Millar AH, Thomson JA, Ren B, Ecker JR. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature. 2009;462(7271):315–22. https://doi.org/10.1038/nature08514.

  89. Ball MP, Li JB, Gao Y, Lee JH, LeProust EM, Park IH et al. Targeted and genome-scale strategies reveal gene-body methylation signatures in human cells. Nat Biotechnol. 2009;27:361–68.

  90. Aran D, Toperoff G, Rosenberg M, Hellman A. Replication timing-related and gene body-specific methylation of active human genes. Hum Mol Genet. 2011;20:670–80.

  91. Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012;13(7):484–92. https://doi.org/10.1038/nrg3230.

  92. Smith J, et al. Promoter DNA hypermethylation and paradoxical gene activation. Trends Cancer. 2020;6(5):392–406.

    Article  CAS  PubMed  Google Scholar 

  93. Stefanetti RJ, et al. Recent advances in understanding the role of FOXO3. F1000Research. 2018. https://doi.org/10.12688/f1000research.15258.1.

    Article  PubMed  PubMed Central  Google Scholar 

  94. Bui TT, et al. γ-Glutamyl transferase 7 is a novel regulator of glioblastoma growth. BMC Cancer. 2015;15:225.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  95. Skelly LE, et al. Effect of sex on the acute skeletal muscle response to sprint interval exercise. Exp Physiol. 2017;102(3):354–65.

    Article  CAS  PubMed  Google Scholar 

  96. Welle S, et al. Skeletal muscle gene expression profiles in 20–29 year old and 65–71 year old women. Exp Gerontol. 2004;39(3):369–77.

    Article  CAS  PubMed  Google Scholar 

  97. Voisin S, et al. Meta-analysis of genome-wide DNA methylation and integrative omics of age in human skeletal muscle. J Cachexia Sarcopenia Muscle. 2021;12(4):1064–78.

    Article  PubMed  PubMed Central  Google Scholar 

  98. Baldelli S, et al. Glutathione and nitric oxide: key team players in use and disuse of skeletal muscle. Nutrients. 2019;11(10):2318.

    Article  CAS  PubMed Central  Google Scholar 

  99. Yan X, et al. The gene SMART study: method, study design, and preliminary findings. BMC Genomics. 2017;18(Suppl 8):821.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  100. Ribel-Madsen R, et al. Genome-wide analysis of DNA methylation differences in muscle and fat from monozygotic twins discordant for type 2 diabetes. PLoS ONE. 2012;7(12):e51302.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  101. Tian Y, et al. ChAMP: updated methylation analysis pipeline for Illumina BeadChips. Bioinformatics. 2017;33(24):3982–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  102. Pidsley R, et al. Critical evaluation of the Illumina MethylationEPIC BeadChip microarray for whole-genome DNA methylation profiling. Genome Biol. 2016;17(1):1–17.

    Article  CAS  Google Scholar 

  103. Chen Y-A, et al. Cross-reactive DNA microarray probes lead to false discovery of autosomal sex-associated DNA methylation. The American Journal of Human Genetics. 2012;91(4):762–4.

    Article  CAS  PubMed  Google Scholar 

  104. Leek JT, et al. Package ‘sva’. 2014.

  105. van Iterson M, van Zwet EW, Heijmans BT. Controlling bias and inflation in epigenome-and transcriptome-wide association studies using the empirical null distribution. Genome Biol. 2017;18(1):1–13.

    CAS  Google Scholar 

  106. Leek JT, et al. Tackling the widespread and critical impact of batch effects in high-throughput data. Nat Rev Genet. 2010;11(10):733–9.

    Article  CAS  PubMed  Google Scholar 

  107. Willer CJ, Li Y, Abecasis GR. METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics. 2010;26(17):2190–1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  108. Aryee MJ, et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics. 2014;30(10):1363–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  109. Smyth GK, et al. Limma: linear models for microarray data. In: Gentleman R, et al., editors. Bioinformatics and computational biology solutions using R and bioconductor. New York: Springer; 2005. p. 397–420.

    Chapter  Google Scholar 

  110. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B (Methodological). 1995;57:289–300.

    Google Scholar 

  111. Benjamin DJ, et al. Redefine statistical significance. Nature Human. Behaviour. 2018;2(1):6.

    Google Scholar 

  112. Peters TJ, et al. De novo identification of differentially methylated regions in the human genome. Epigenetics Chromatin. 2015;8(1):6.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  113. Ferrari S, Cribari-Neto F. Beta regression for modelling rates and proportions. J Appl Stat. 2004;31(7):799–815.

    Article  Google Scholar 

  114. Van Buuren S, Groothuis-Oudshoorn K. Multivariate imputation by chained equations in R. J Stat Softw. 2011;45:1–67.

    Article  Google Scholar 

  115. Fishilevich S, et al. GeneHancer: genome-wide integration of enhancers and target genes in GeneCards. Database. 2017. https://doi.org/10.1093/database/bax028.

    Article  PubMed  PubMed Central  Google Scholar 

  116. Phipson B, Maksimovic J, Oshlack A. missMethyl: an R package for analyzing data from Illumina’s HumanMethylation450 platform. Bioinformatics. 2016;32(2):286–8.

    Article  CAS  PubMed  Google Scholar 

  117. Maksimovic J, Oshlack A, Phipson B. Gene set enrichment analysis for genome-wide DNA methylation data. bioRxiv, 2020.

  118. Schiaffino S, Reggiani C. Fiber types in mammalian skeletal muscles. Physiol Rev. 2011;91(4):1447–531.

    Article  CAS  PubMed  Google Scholar 

  119. Bloemberg D, Quadrilatero J. Rapid determination of myosin heavy chain expression in rat, mouse, and human skeletal muscle using multicolor immunofluorescence analysis. PLoS ONE. 2012;7(4):e35273–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  120. Godsland IF, et al. The effects of different formulations of oral contraceptive agents on lipid and carbohydrate metabolism. N Engl J Med. 1990;323(20):1375–81.

    Article  CAS  PubMed  Google Scholar 

  121. Urbut SM, et al. Flexible statistical methods for estimating and testing effects in genomic studies with multiple conditions. Nat Genet. 2019;51(1):187–95.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We would like to thank Dr Alba Moreno-Asso from Victoria University for her assistance with the qPCR analysis and Dr Andrew Garnham for his expertise performing the muscle biopsies. The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. The data used for the analyses described in this manuscript were obtained from the GTEx Portal on 08/26/2020.

Funding

We are grateful for the support of the Australian National Health and Medical Research Council (NHMRC) via Sarah Voisin’s Early Career Research Fellowship (APP11577321) and Nir Eynon’s NHMRC Investigator Grant (APP1194159), as well as the Australian Research Council (ARC) (DP190103081). The Gene SMART study and Nick Harvey's stipend was supported by the Collaborative Research Network for Advancing Exercise and Sports Science (201202) from the Department of Education and Training, Australia. This research was also supported by infrastructure purchased with Australian Government EIF Super Science Funds as part of the Therapeutic Innovation Australia—Queensland Node project (LRG).

Author information

Author notes

  1. Sarah Voisin and Nir Eynon are co-senior authors

Authors and Affiliations

  1. Institute for Health and Sport (iHeS), Victoria University, PO Box 14428, Melbourne, VIC, 8001, Australia

    Shanie Landen, Macsue Jacques, Danielle Hiam, Javier Alvarez-Romero, Sarah Voisin & Nir Eynon

  2. Faculty of Health Sciences and Medicine, Bond University, Gold Coast, QLD, 4226, Australia

    Nicholas R. Harvey & Kevin J. Ashton

  3. Centre for Genomics and Personalised Health, Genomics Research Centre, School of Biomedical Sciences, Queensland University of Technology (QUT), 60 Musk Ave., Kelvin Grove, QLD, 4059, Australia

    Nicholas R. Harvey, Larisa M. Haupt & Lyn R. Griffiths

  4. Institute for Physical Activity and Nutrition, School of Exercise and Nutrition Sciences, Deakin University, Geelong, Australia

    Danielle Hiam & Séverine Lamon

Authors
  1. Shanie Landen
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Macsue Jacques
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Danielle Hiam
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Javier Alvarez-Romero
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Nicholas R. Harvey
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Larisa M. Haupt
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Lyn R. Griffiths
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Kevin J. Ashton
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Séverine Lamon
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Sarah Voisin
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Nir Eynon
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

S.L., S.V., and N.E. contributed to conceptualization. S.L., and S.V. helped in methodology. S.L., S.V. contributed to investigation. S.L. contributed to formal analysis. S.L., M.J., D.H., J.A.R., N.R.H., L.M.H., L.R.G., K.J.A., and N.E. contributed to resources. S.L., S.V., D.H., and N.E. helped in writing—original draft. S.L., S.V., D.H., S.Lamon, N.R.H., L.M.H., K.J.A., and N.E. contributed to writing—review and editing. S.V., L.R.G., and N.E. helped in funding acquisition. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Nir Eynon.

Ethics declarations

Consent for publication

Not applicable.

Competing interest

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

. Supplementary figures. S1-Correlation plots of chi2 tests of genomic locations. S2-Fibre type proportion analysis. S3-Quantile-quantile plot of −log10 transformed p-values for hormone meta-analysis in the Gene SMART study and sex meta-analysis in the three datasets included in the epigenome-wide association study. S4-Principal component analysis (PCA) of beta values of all tested CpGs across the participants in the FUSION, Gene SMART, and GSE38291 cohorts. S5-Distribution of the 10,000 random permutations for a negative correlation between DNA methylation and gene expression. S6-Gene expression for FOXO3 and ALDH1A1 validated in 3 cohorts. S7-Comparison of results from the full meta-analysis and from a meta-analysis excluding T2D participants in FUSION.

Additional file 2

. Study descriptions. Description of participants, study design, muscle collection, and data pre-processing/analysis.

Additional file 3

. DMPs. Differentially methylated positions between males and females in the meta-analysis FDR<0.005. Corresponding chromosome, genomic location, annotated genes, male and female chromatin states from the Roadmaps Epigenomics Project, and genes annotated by GeneHancer. Positive effect size indicates higher DNA methylation in males compared to females.

Additional file 4

. DMRs. Differentially methylated regions between males and females in the meta-analysis Stouffer, HMFDR, and Fisher <0.005. Corresponding chromosome, genomic location, width of DMR, number of CpGs in DMR, statistics (Stouffer, Harmonic mean false discovery rate (HMFDR), and Fisher statistic), maximum and mean DMR effect sizes, and annotated genes. Positive effect size indicates higher DNA methylation in males compared to females.

Additional file 5

. Overlapping genes. Genes which displayed sex-biased gene expression in GTEx and FUSION as well as sex-biased DNA methylation (according to DMRs) in the meta-analysis. Corresponding chromosome, ensemble gene ID, gene name, GTEx mash posterior effect size, GTEx local false sign rate threshold, FUSION mRNA fold change, FUSION mRNA FDR, and number of DMPs per gene. Positive effect sizes indicate higher DNA methylation or expression in males compared to females.

Additional file 6

. GO (DMPs). Gene Ontology terms identified with GSEA using the differentially methylated positions. Type of GO term- biological process (BP), molecular function (MF), and cellular component (CC), name of GO term, N represents the total number of genes in the GO term, DE represents the number of differentially methylated genes in the GO term, and SigGenesInSet are the differentially methylated genes in the GO term.

Additional file 7

. GO (DMRs). Gene Ontology terms identified with GSEA using the differentially methylated regions. Type of GO term- biological process (BP), molecular function (MF), and cellular component (CC), name of GO term, N represents the total number of genes in the GO term, DE represents the number of differentially methylated genes in the GO term, and SigGenesInSet are the differentially methylated genes in the GO term.

Additional file 8

. KEGG (DMPs). Kyoto Encyclopedia of Genes and Genomes pathways identified with GSEA using the differentially methylated positions. Description of KEGG pathway, N represents the total number of genes in the KEGG pathway, DE represents the number of differentially methylated genes in the KEGG pathway, and SigGenesInSet are the differentially methylated genes in the KEGG pathway.

Additional file 9

. KEGG (DMRs). Kyoto Encyclopedia of Genes and Genomes pathways identified with GSEA using the differentially methylated regions. Description of KEGG pathway, N represents the total number of genes in the KEGG pathway, DE represents the number of differentially methylated genes in the KEGG pathway, and SigGenesInSet are the differentially methylated genes in the KEGG pathway.

Additional file 10

. Reactome (DMPs). Reactome pathways identified with GSEA using the differentially methylated positions. Description of Reactome pathway, N represents the total number of genes in the Reactome pathway, DE represents the number of differentially methylated genes in the Reactome pathway, and SigGenesInSet are the differentially methylated genes in the Reactome pathway.

Additional file 11

. Reactome (DMRs). Reactome pathways identified with GSEA using the differentially methylated regions. Description of Reactome pathway, N represents the total number of genes in the Reactome pathway, DE represents the number of differentially methylated genes in the Reactome pathway, and SigGenesInSet are the differentially methylated genes in the Reactome pathway.

Additional file 12

. PCRs. Results from qPCR of GGT7, FOXO3, and ALDH1A1 in Gene SMART cohort. Blank cells indicate sample that had a standard deviation of greater than 1 Ct between triplicates. Housekeeping genes tested are Cyclophilin, 18s rRNA, and TBP. Primer sequences and PCR conditions used.

Additional file 13

. Concentrations of circulating testosterone (nmol/L), free testosterone (pmol/L), sex hormone-binding globulin (nmol/L), and estrogen (pmol/L) in males and females from the Gene SMART cohort. Includes 20 females (20 at rest before four weeks of exercise training, same 20 at rest after four weeks of exercise training, and six of those at rest before a control period) and 44 males (we did not have data for one male included in the DNA methylation analysis).

Additional file 14

. Data available for each of the datasets included in the DNA methylation meta-analysis. Immunohistochemistry (IHC); sex hormone-binding globulin (SHBG), free testosterone (Free T), testosterone (T), estradiol (E2)

Additional file 15

. List of transcription factors (TFs) included in analysis for enrichment of transcription factor binding sites (TFBSs) among differentially methylated positions (DMPs). The current UniBind database tests a total of 268 unique TFs from 518 different cell types.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Landen, S., Jacques, M., Hiam, D. et al. Skeletal muscle methylome and transcriptome integration reveals profound sex differences related to muscle function and substrate metabolism. Clin Epigenet 13, 202 (2021). https://doi.org/10.1186/s13148-021-01188-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13148-021-01188-1

Keywords

Clinical Epigenetics

ISSN: 1868-7083

Contact us