- Open Access
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)
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.
Various diseases, including but not limited to muscular dystrophy, cardiomyopathies, and cardiovascular disease , 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 . Sexually differentiated traits and phenotypes stem from a combination of factors, including genetics (gene variants-by-sex interactions , XY chromosome complements [4,5,6,7], genomic imprinting ), the hormonal milieu [9, 10], and gene regulation , with the latter likely contributing the most .
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 . 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 , all of which inherently contribute to enhanced fatigue-resistance in female skeletal muscle . As such, females exhibit higher mRNA and protein levels of lipid oxidation-related genes than males . 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) . 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 . 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 , blood [19, 20], and more recently cultured myoblasts and myotubes . 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 , which are primarily ascribed to the reproductive function, but their importance to non-reproductive functions, such as skeletal muscle , 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 , 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 , as well as different contractile and metabolic properties . 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.
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).
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.
Each data set had a unique study design that required considering different independent variables known to affect DNA methylation, such as age  and type 2 diabetes (T2D)  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 , and type II fibres exhibit hypomethylation compared to type I fibres . 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).
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 . 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 . 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) . 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”).
Given the role of transcription factors (TFs) in regulating chromatin accessibility and thus effecting downstream gene expression , 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.
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 ), 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.  (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).
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 . 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).
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 . 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).
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 , including blood [54, 55] and pancreatic islets . 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 , and as females tend to have a higher proportion of type I fibres than males . 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 . 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 , 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 .
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 , and SMAD3 , 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  and gene expression . Differences in sex hormone levels between developing males and females are already evident in utero , 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  and transgender people  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 . 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 .
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 . 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 , 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 .
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 . 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 .
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 . 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.  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 ; 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 ; and GGT7 for its role in antioxidant activity . 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 . 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  and that elderly females have higher skeletal muscle FOXO3 expression than younger females . 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 . 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 .
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.
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 , the Finland–United States Investigation of NIDDM Genetics (FUSION) study from the dbGAP repository (phs000867.v1.p1) , and the GSE38291 dataset from the Gene Expression Omnibus (GEO) platform . 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
The pre-processing of DNA methylation data was performed according to the bioinformatics pipeline developed for the Bioconductor project . 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  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 . 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 .
We adjusted each EWAS for bias and inflation using the empirical null distribution as implemented in BACON . Inflation and bias in EWAS are caused by unmeasured technical and biological confounding, such as population substructure, batch effects, and cellular heterogeneity . 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 . The results were consistent with the inflation factors and biases reported in an EWAS of age in blood . 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 . 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 , 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 , 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:
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:
For the GSE38291 study, the linear model was of the form:
All results were adjusted for multiple testing using the Benjamini and Hochberg correction , and all CpGs showing an FDR < 0.005 were considered significant . DMRs were identified using the DMRcate package with default settings . 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  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 .
To annotate the CpGs to the genome (see “Availability of Data and Materials”), we integrated a comprehensive annotation of Illumina HumanMethylation arrays  with chromatin states from the Roadmap Epigenomics Project  and the latest GeneHancer information . 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  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 . 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 , 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 . 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 . 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 . 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 ; 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 . GTEx effect sizes are represented as mash posterior effect sizes , 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.
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).
Khramtsova EA, Davis LK, Stranger BE. The role of sex in the genomics of human complex traits. Nat Rev Genet. 2018;1:173–90.
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.
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.
Oliva M, et al. The impact of sex on gene expression across human tissues. Science. 2020;369(6509):eaba3066.
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.
Arnold AP. Y chromosome’s roles in sex differences in disease. Proc Natl Acad Sci. 2017;114(15):3787–9.
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.
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.
Pihlajamaa P, et al. Tissue-specific pioneer factors associate with androgen receptor cistromes and transcription programs. EMBO J. 2014;33(4):312–26.
Varlamov O, Bethea CL, Roberts CT Jr. Sex-specific differences in lipid and glucose metabolism. Front Endocrinol. 2015;5:241.
Lopes-Ramos CM, et al. Sex differences in gene expression and regulatory networks across 29 human tissues. Cell Rep. 2020;31(12):107795.
Maher AC, et al. Sex differences in global mRNA content of human skeletal muscle. PLoS ONE. 2009;4(7):e6335.
Tarnopolsky M. Sex differences in exercise metabolism and the role of 17-beta estradiol. Med Sci Sports Exerc. 2008;40(4):648–54.
Landen S, et al. Genetic and epigenetic sex-specific adaptations to endurance exercise. Epigenetics. 2019;14(6):523–35.
Haizlip KM, Harrison BC, Leinwand LA. Sex-based differences in skeletal muscle kinetics and fiber-type composition. Physiology (Bethesda). 2015;30(1):30–9.
Hunter SK. Sex differences in human fatigability: mechanisms and insight to physiological responses. Acta Physiol (Oxf). 2014;210(4):768–89.
Moore LD, Le T, Fan G. DNA methylation and its basic function. Neuropsychopharmacology. 2013;38(1):23–38.
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.
Singmann P, et al. Characterization of whole-genome autosomal differences of DNA methylation between men and women. Epigenet Chromatin. 2015;8(1):1–13.
Liu J, et al. A study of the influence of sex on genome wide methylation. PLoS ONE. 2010;5(4):e10028.
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.
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.
Welle S, Tawil R, Thornton CA. Sex-related differences in gene expression in human skeletal muscle. PLoS ONE. 2008;3(1):e1385.
Gershoni M, Pietrokovski S. The landscape of sex-differential transcriptome and its consequent selection in human adults. BMC Biol. 2017;15(1):7.
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.
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.
Pihlajamaa P, Sahu B, Jänne OA. Determinants of receptor- and tissue-specific actions in androgen signaling. Endocr Rev. 2015;36(4):357–84.
Varlamov O, Bethea CL, Roberts CT Jr. Sex-specific differences in lipid and glucose metabolism. Front Endocrinol (Lausanne). 2014;5:241.
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.
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.
Wu Y, et al. Identification of androgen response elements in the insulin-like growth factor I upstream promoter. Endocrinology. 2007;148(6):2984–93.
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.
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.
Mayne BT, et al. Large scale gene expression meta-analysis reveals tissue-specific, sex-biased gene expression in humans. Front Genet. 2016;7:183.
Maughan R, Watson JS, Weir J. Strength and cross-sectional area of human skeletal muscle. J Physiol. 1983;338(1):37–49.
Carter S, et al. Changes in skeletal muscle in males and females following endurance training. Can J Physiol Pharmacol. 2001;79(5):386–92.
Begue G, et al. DNA methylation assessment from human slow-and fast-twitch skeletal muscle fibers. J Appl Physiol. 2017;122(4):952–67.
Zierath JR, Hawley JA. Skeletal muscle fiber type: influence on contractile and metabolic properties. PLoS Biol. 2004;2(10):e348.
Voisin S, et al. Meta-analysis of genome-wide DNA methylation and integrative OMICs in human skeletal muscle. bioRxiv, 2020.
Davegårdh C, et al. DNA methylation in the pathogenesis of type 2 diabetes in humans. Mol Metab. 2018;14:12–25.
Jones PA, Takai D. The role of DNA methylation in mammalian epigenetics. Science. 2001;293(5532):1068–70.
Kundaje A, et al. Integrative analysis of 111 reference human epigenomes. Nature. 2015;518(7539):317–30.
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.
Voss TC, Hager GL. Dynamic regulation of transcriptional states by chromatin and transcription factors. Nat Rev Genet. 2014;15(2):69–81.
Gheorghe M, et al. A map of direct TF–DNA interactions in the human genome. Nucl Acids Res. 2019;47(4):e21–e21.
Puig RR, et al. UniBind: maps of high-confidence direct TF-DNA interactions across nine species. bioRxiv, 2020.
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.
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.
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.
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.
Miller AEJ, et al. Gender differences in strength and muscle fiber characteristics. Eur J Appl Physiol. 1993;66(3):254–62.
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.
McCarthy NS, et al. Meta-analysis of human methylation data for evidence of sex-specific autosomal patterns. BMC Genom. 2014;15(1):981.
Yousefi P, et al. Sex differences in DNA methylation assessed by 450 K BeadChip in newborns. BMC Genom. 2015;16(1):911.
Inoshita M, et al. Sex differences of leukocytes DNA methylation adjusted for estimated cellular proportions. Biol Sex Differ. 2015;6(1):1–7.
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.
Suzuki M, Bird A, Suzuki MM, Bird A. DNA methylation landscapes: provocative insights from epigenomics. Nat Rev Genet. 2008;9:465–76.
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.
Rubenstein AB, et al. Single-cell transcriptional profiles in human skeletal muscle. Sci Rep. 2020;10(1):1–15.
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.
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.
Suzuki MM, Bird A. DNA methylation landscapes: provocative insights from epigenomics. Nat Rev Genet. 2008;9(6):465–76.
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.
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.
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.
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.
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.
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.
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.
Amar D, et al. Time trajectories in the transcriptomic response to exercise-a meta-analysis. Nat Commun. 2021;12(1):1–12.
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.
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.
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.
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.
Lundsgaard A-M, Kiens B. Gender differences in skeletal muscle substrate metabolism–molecular mechanisms and insulin sensitivity. Front Endocrinol. 2014;5:195.
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.
De Vries GJ. Minireview: sex differences in adult and developing brains: compensation, compensation. Compens Endocrinol. 2004;145(3):1063–8.
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.
Penny GD, et al. Requirement for Xist in X chromosome inactivation. Nature. 1996;379(6561):131–7.
Carrel L, Willard HF. X-inactivation profile reveals extensive variability in X-linked gene expression in females. Nature. 2005;434(7031):400–4.
Sidorenko J, et al. The effect of X-linked dosage compensation on complex trait variation. Nat Commun. 2019;10(1):3009.
Tukiainen T, et al. Landscape of X chromosome inactivation across human tissues. Nature. 2017;550(7675):244–8.
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.
Trolle C, et al. Widespread DNA hypomethylation and differential gene expression in turner syndrome. Sci Rep. 2016;6(1):1–14.
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.
Villicaña S, Bell JT. Genetic impacts on DNA methylation: research findings and future perspectives. Genome Biol. 2021;22(1):1–35.
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.
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.
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.
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.
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.
Smith J, et al. Promoter DNA hypermethylation and paradoxical gene activation. Trends Cancer. 2020;6(5):392–406.
Stefanetti RJ, et al. Recent advances in understanding the role of FOXO3. F1000Research. 2018. https://doi.org/10.12688/f1000research.15258.1.
Bui TT, et al. γ-Glutamyl transferase 7 is a novel regulator of glioblastoma growth. BMC Cancer. 2015;15:225.
Skelly LE, et al. Effect of sex on the acute skeletal muscle response to sprint interval exercise. Exp Physiol. 2017;102(3):354–65.
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.
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.
Baldelli S, et al. Glutathione and nitric oxide: key team players in use and disuse of skeletal muscle. Nutrients. 2019;11(10):2318.
Yan X, et al. The gene SMART study: method, study design, and preliminary findings. BMC Genomics. 2017;18(Suppl 8):821.
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.
Tian Y, et al. ChAMP: updated methylation analysis pipeline for Illumina BeadChips. Bioinformatics. 2017;33(24):3982–4.
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.
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.
Leek JT, et al. Package ‘sva’. 2014.
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.
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.
Willer CJ, Li Y, Abecasis GR. METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics. 2010;26(17):2190–1.
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.
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.
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.
Benjamin DJ, et al. Redefine statistical significance. Nature Human. Behaviour. 2018;2(1):6.
Peters TJ, et al. De novo identification of differentially methylated regions in the human genome. Epigenetics Chromatin. 2015;8(1):6.
Ferrari S, Cribari-Neto F. Beta regression for modelling rates and proportions. J Appl Stat. 2004;31(7):799–815.
Van Buuren S, Groothuis-Oudshoorn K. Multivariate imputation by chained equations in R. J Stat Softw. 2011;45:1–67.
Fishilevich S, et al. GeneHancer: genome-wide integration of enhancers and target genes in GeneCards. Database. 2017. https://doi.org/10.1093/database/bax028.
Phipson B, Maksimovic J, Oshlack A. missMethyl: an R package for analyzing data from Illumina’s HumanMethylation450 platform. Bioinformatics. 2016;32(2):286–8.
Maksimovic J, Oshlack A, Phipson B. Gene set enrichment analysis for genome-wide DNA methylation data. bioRxiv, 2020.
Schiaffino S, Reggiani C. Fiber types in mammalian skeletal muscles. Physiol Rev. 2011;91(4):1447–531.
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.
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.
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.
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.
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).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
. 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.
. Study descriptions. Description of participants, study design, muscle collection, and data pre-processing/analysis.
. 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.
. 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.
. 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.
. 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.
. 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.
. 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.
. 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.
. 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.
. 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.
. 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.
. 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).
. 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)
. 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.
About this article
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
- DNA methylation
- Sex differences
- Skeletal muscle
- Gene expression
- Epigenome-wide study (EWAS)