key: cord-0304674-d7pkp5e8
authors: Landen, Shanie; Jacques, Macsue; Hiam, Danielle; Alvarez-Romero, Javier; Harvey, Nicholas R; Haupt, Larisa M.; Griffiths, Lyn R; Ashton, Kevin J; Lamon, Séverine; Voisin, Sarah; Eynon, Nir
title: Skeletal muscle methylome and transcriptome integration reveals profound sex differences related to muscle function and substrate metabolism
date: 2021-08-25
journal: bioRxiv
DOI: 10.1101/2021.03.16.435733
sha: f3369f94ac372eb7bf1734123eab3b993d1f2d91
doc_id: 304674
cord_uid: d7pkp5e8

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, 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 8,420 genes that exhibit DNA methylation differences between males and females in human skeletal muscle that may modulate mechanisms controlling muscle metabolism and health. 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 genomewide sex-biased DNA methylation and expression in skeletal muscle, shedding light on distinct sex differences in skeletal muscle.

underlying DNA methylation sex differences and found that circulating hormones were not 48 associated with differential methylation at sex-biased DNA methylation loci, however these 49 sex-specific loci were enriched for binding sites of hormone-related transcription factors (with 50 top TFs including androgen (AR), estrogen (ESR1), and glucocorticoid (NR3C1) receptors). 51

Fibre type proportions were associated with differential methylation across the genome, as well 52 as across 16 % of sex-biased DNA methylation loci (FDR < 0.005). Integration of DNA 53 methylomic results with transcriptomic data from the GTEx database and the FUSION cohort 54 revealed 326 autosomal genes that display sex differences at both the epigenome and 55 transcriptome levels. Importantly, transcriptional sex-biased genes were overrepresented 56 Introduction 65 66

Sex differences are evident in nearly all complex traits. Various diseases, including but 67 not limited to cancer, muscular dystrophy, and 4] , display sex differences in 68 prevalence, onset, progression, or severity. To improve treatment for such diseases, it is crucial 69 to uncover the molecular basis for the sex differences and their consequences on organ function 70 [5] . Sexually differentiated traits and phenotypes stem from a combination of factors, including 71 genetics (gene variants-by-sex interactions [6] , XY chromosome complements [7] [8] [9] [10] , genomic 72 imprinting [11] ), the hormonal milieu [12, 13] , and gene regulation [14] , with the latter likely 73 contributing the most [3] . 74

Recently, a large-scale study from the Genotype-Tissue Expression (GTEx) consortium 75 unravelled mRNA expression differences between the sexes that are not driven by sex 76 chromosomes, across all tissues. Skeletal muscle was particularly divergent between the sexes, 77 as gene expression profiles in this tissue could predict sex with high specificity and sensitivity 78 [7] . These transcriptomic differences underpin the numerous physiological differences in 79 skeletal muscle between males and females, such as differences in substrate metabolism [15-80 17] . For example, females oxidize more lipids and less carbohydrates and amino acids during 81 CpG (633, 645 in total) and those that appear in color 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 labeled 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.

Males typically show a greater proportion of type II muscle fibres compared with 157 females [18] , and type II fibres exhibit hypomethylation compared to type I fibres [40] . 158 Therefore, we hypothesized that the observed DNA methylation sex differences, specifically 159 the hypomethylation in males, may be a result of differing fibre type distributions between 160 males and females. We first estimated type I fibre proportions in the Gene SMART cohort via 161 immunohistochemistry (Supplementary Figure 2D ) and the FUSION cohort via RNA-seq (data 162 available for each dataset in Supplementary Table 13 ). In both the Gene SMART and FUSION 163 cohorts, females had higher proportions of type I fibres than males (Supplementary Figure  164 2B/C). We could not directly add fibre type proportions to the linear model as a covariate, since 165 fibre type proportions are not a confounder (i.e. a factor that influences both sex and DNA 166 methylation independently), but may be a direct downstream effect of sex, in turn affecting 167 DNA methylation. Adding fibre type proportions in the model would therefore distort the 168 association between sex and DNA methylation. To overcome this issue, we stratified the 169 cohorts by sex, added fibre type proportions to the model as a covariate and identified DNA 170 methylation patterns associated with fibre type proportions. We then meta-analysed the results 171 to find CpGs robustly associated with fibre type proportions across both cohorts and all sexes 172 (see "Methods"). We identified 16,275 CpGs associated with fibre type proportions 173 (Supplementary Figure 2) . When restricting the analysis to the loci exhibiting sex-biased DNA 174 methylation, 8,805 (15.5%) of those were associated with fibre type proportions (FDR < 175 0.005). Effect sizes ranged from -0.28% to +0.30% DNA methylation difference per % increase 176 in type I fibre content (Figure 2A) . 177

We then aimed to determine whether circulating sex hormone levels underlie the 178 observed DNA methylation sex differences. We analysed oestrogen (as oestradiol, E2), 179 testosterone (T), free testosterone (Free T), and sex-hormone binding globulin (SHBG) levels 180 using mass-spectrometry (Free T derived from calculation) of blood serum in males and 181 females in the Gene SMART cohort. Males and females significantly differed in all four 182 hormone levels (Supplementary Table 12 ). To avoid collinearity with sex, we separated males 183 and females for the association between sex-differential DNA methylation and circulating 184 hormones. We assessed whether each of the four hormone levels was associated with DNA 185 methylation across all of the CpGs and across the sex-DMPs in each sex by adjusting the linear 186 model for a given hormone. In both males and females, circulating free testosterone, 187 testosterone, oestrogen, and SHBG levels were not highly associated with DNA methylation 188 (less than five DMPs; FDR < 0.005) of neither all of the CpGs tested nor the sex-DMPs 189 previously identified (Supplementary Figure 3) . To limit the potentially confounding effect of 190 fluctuating ovarian hormone levels on DNA methylation, female muscle biopsies were 191 collected in the early follicular phase of the menstrual cycle (days 1-7) and blood serum were 192 tested for follicle stimulating hormone (FSH), luteinizing hormone (LH), and progesterone (as 193 well as E2 as previously mentioned). None of the four hormones (when adjusting for each 194 hormone separately or when adjusting for the first two principal components; see 'methods'), 195 were associated with differential methylation of all the CpGs or the sex-DMPs. This suggests 196 that variations in ovarian hormone levels in the early follicular phase did not confound our 197 results. 198 199

Since the association between DNA methylation and gene expression depends on the 202 genomic context, we inspected the genomic distribution of the DMPs to gain insights into their 203 potential function [44] . We compared the distribution of hyper-, hypo-, and non-DMPs among 204 the various chromatin states in human skeletal muscle using the Roadmap Epigenomics project 205

[45]. DMPs were not randomly distributed in the chromatin states (χ 2 p-value < 2.2 x 10 -16 , 206 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-DMP which showed the largest effect sizes are labeled 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. which yielded equivalent findings. We noted that DMPs were overrepresented in loci whose 214 chromatin states differ between males and females (38.7 % of DMPs vs. 32.4% of non-DMPs 215 are in chromatin states that differ between males and females), which means that the odds of a 216 DMP being located in a sex-differing chromatin state increased by a factor of 1.3 compared 217 with a non-DMP (OR = 0.76, 95% confidence interval = 0.75-0.77, Fisher test p-value < 2.2e-218 16) ( Figure 2B ). DMPs were also enriched in CpG island shores and depleted in CpG islands 219 (χ 2 p-value < 2.2e-16) ( Figure 2C, Supplementary figure 1B) . To gain insights into the potential downstream effects of sex-biased DNA methylation 241 on gene expression, we integrated results from the EWAS meta-analysis of sex with genes 242 whose mRNA expression levels are known to differ between males and females. We 

We sought to validate the sex-biased gene expression obtained from GTEx in a subset 279 of the samples used for methylation analysis since the DMGs and DEGs analyses were obtained 280 from different muscle groups (the DMGs of the current study are from the vastus lateralis while 281 the GTEx DEGs are from the gastrocnemius). Although both are skeletal muscle tissue from 282 the leg, there may be differences in muscle phenotypes in differing muscle groups [52] . 283

Analysis of RNA sequencing data from the FUSION cohort revealed 3,751 autosomal genes 284 with sex-biased expression (FDR < 0.005). The FDR threshold we chose for the FUSION gene 285 expression data was more stringent than the GTEx local false sign rate threshold (lfsr < 0.05), 286 yet, ~34% of the genes which were both DEGs in GTEx and DMGs were also DEGs in the 287 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. FUSION cohort, totalling 326 genes (hereinto referred to as `overlapping genes`) ( Figure 3A) . 288

Given that both the GTEx and FUSION cohorts include participants of relatively older ages, 289 we sought to confirm the mRNA levels in the younger cohort in the analysis (the Gene 290 SMART) for three genes that displayed sex differences at both the mRNA and DNA 291 methylation levels (GGT7, FOXO3, and ALDH1A1) ( Table 2, Supplementary Table 11 

Three-hundred twenty-six genes exhibited differential methylation in the meta-analysis 300 and differential expression among the GTEx and FUSION cohorts, termed `overlapping genes`. 301

Of those genes, we tested three for gene expression levels, GGT7, Forkhead Box O3 (FOXO3), 302

and Aldehyde Dehydrogenase 1 Family Member A1 (ALDH1A1), in the younger cohort 303 included in the DNA methylation analysis (Gene SMART) given the effect that age has on 304 skeletal muscle gene expression [53] . These three genes showed a large effect size in gene 305 expression and DNA methylation, displayed moderate gene expression levels in skeletal 306 muscle relative to other tissues, and/or contained numerous DMPs and DMRs ( 

We conducted a large-scale meta-analysis of DNA methylation differences between 322 males and females in skeletal muscle, and integrated them with transcriptomic data. We 323 revealed that males display profound genome-wide hypomethylation compared with females, 324 and that hormone-related TFBSs and muscle fibre type proportions underlie the observed DNA 325 methylation sex differences. We also showed that many sex-biased genes found in GTEx 326 exhibit sex-biased DNA methylation, which was partially confirmed in the FUSION cohort. 327

We then validated the gene expression (qPCR) levels of three genes with large DNA 328 methylation and expression differences between the sexes across cohorts, and confirmed the 329 higher gene expression in males of GGT7 and ALDH1A1. Finally, we showed that the DMGs 330 are overwhelmingly involved in muscle contraction, as well as substrate metabolism and 331 anatomical structure-related pathways. 332

In the present study, the overwhelming majority (94%) of the DMPs were 333 hypomethylated in males. Interestingly, global autosomal hypomethylation in males has been 334 observed in various other tissues [54] , including blood [55, 56] and pancreatic islets [20] . We 335 aimed to uncover some of the biological mechanisms at the root of these epigenetic sex-336 differences by assessing the contributions of fibre type proportions and sex hormone levels to 337 the observed DNA methylation sex differences. We hypothesised that differences in fibre type 338

proportions between sexes may partly explain our findings [56] [57] [58] , as studies report that type 339 I fibres are hypermethylated compared with type II fibres [40] , and as females tend to have a 340 higher proportion of type I fibres than males [18] . Consistent with this, we observed that 341 females had higher proportions of type I muscle fibres than males and that type I fibre content 342 was mostly associated with DNA hypermethylation. Importantly, 16% of the loci exhibiting 343 sex-biased DNA methylation were also associated with fibre type proportions. This suggests 344 that at those CpGs, differences in DNA methylation between the sexes is due to the inherent 345 sex differences in fibre type proportions. Nonetheless, the vast majority of the loci that exhibit 346 sex-biased DNA methylation (84%; 48,008 CpGs) differ regardless of the sex differences in 347 fibre type proportions. A recent study on the FUSION cohort, adjusted for fibre type 348

proportions and found that it explains a substantial portion of the variability in DNA 349 methylation for many metabolic phenotypes of interest [59] . The enrichment for muscle contraction-related pathways among the DMGs across GO, 389 KEGG, and Reactome suggests that the sex-differential DNA methylation has functional 390 relevance in skeletal muscle function. Furthermore, enrichment of substrate metabolism 391 pathways among the DMGs suggests that the observed sex differences in carbohydrate, lipid, 392 and protein metabolism [16, 76] may have a molecular, more specifically an epigenetic, basis. 393 This is corroborated with results from transcriptomic studies, which report that skeletal muscle 394 female-biased genes are enriched for pathways involved in fatty acid metabolism while male-395 biased genes are enriched for pathways involved in protein catabolism (Lindholm et al., 2014a) . 396

Although not well-understood, the sex chromosome complement may also influence 397 autosomal DNA methylation patterns. In cultured fibroblasts, the presence of Sex-determining 398

Region Y (SRY) is associated with lower autosomal methylation levels [77] [78] [79] . Additionally, 399 a higher number of the X chromosomes, in the absence of SRY, leads to increased methylation 400 levels at a specific sex-differentially methylated autosomal region [79] . This could be attributed 401 to allele dosage compensation, a female-specific process that silences one of the X 402 chromosomes in a cell [80, 81] . Approximately one-third of genes 'escape' inactivation, remain 403 transcriptionally active in XX cells, [81] [82] [83] , and have been suggested to affect autosomal DNA 404 methylation via their histone marks [79, 84] . Moreover, females with Turner syndrome 405 (partially/fully missing one X) and monosomy X have lower global methylation than XX 406 females, but higher than XY males [85, 86] . The effect of sex chromosome complement on 407 autosomal DNA methylation in skeletal muscle has yet to be explored. Finally, genetic variants 408 (copy number variants (CNVs) and single nucleotide polymorphisms (SNPs)) may affect DNA 409 methylation status at specific loci [87], termed methylation-quantitative trait loci (meQTL). 410 However, with the exception of a female bias for large, rare CNVs, there are no large sex 411 differences in autosomal SNP minor allele frequencies [3] . 412

The relationship between DNA methylation and gene expression is complex; DNA 413 methylation at promoters, enhancers, and 1 st exons is generally believed to enhance gene 414 silencing, while DNA methylation at gene bodies can sometimes be associated with increased 415 gene expression (Rauch et al., 2009 , Lister et al., 2009 , Ball et al., 2009 , Aran et al., 2011 , 416 Jones, 2012 . Using a permutation test, we showed that DNA methylation differences between 417 the sexes at promoters and enhancers were more often associated with lower gene expression 418 than would be expected by chance alone. DNA methylation differences between the sexes were 419 also particularly prominent in chromatin states that are known differ between males and 420 females. This suggests that DNA methylation differences between males and females reflect 421 alterations in chromatin activity, and differential epigenetic states and expression are likely 422 functionally connected. In line with this, chromatin states that differ between the sexes have 423 been shown to be enriched for sex-biased genes across various tissues, including skeletal 424 muscle [46] . However, it is not yet possible to assess whether the relationship reflects 425 correlation or direct causality. There is still debate around whether epigenomic features drive 426 regulatory processes or are merely a consequence of transcription factor binding [46] . 427

We identified 326 genes with consistent differential skeletal muscle DNA methylation 428 and expression across 1,172 individuals altogether (369 individuals from three cohorts for 429 DNA methylation and 1,077 individuals from two cohorts for gene expression). Although we 430 found profound global DNA hypomethylation in males, of the overlapping genes there were 431 equivalent numbers of genes over-and under-expressed in males compared with females for 432 both GTEx and FUSION. Indeed, hypermethylation is not always associated with decreased 433 gene expression [88] . The substantial overlap between differentially methylated genes and 434 differentially expressed genes highlights many genes that may be of interest for their roles in 435 muscle-related processes. We focused on three of these genes that displayed a large DNA 436 methylation difference between males and females, are highly expressed in skeletal muscle, or 437 play a role in skeletal muscle function: HDAC4 given its role in neurogenic muscle atrophy 438 [89, 90] and GGT7 for its role in antioxidant activity [95] . Of the three genes which were validated 445 across GTEx, FUSION, and Gene SMART, two of them showed consistently higher male 446 expression levels (GGT7, ALDH1A1) while one showed opposite sex-biased expression 447 (FOXO3) in the young versus the old cohorts. FOXO3 expression was lower in males in the 448 older cohorts (GTEx and FUSION) , and higher in males in the younger cohort (Gene SMART). 449

Other studies have shown that males have higher FOXO3 expression in young skeletal muscle 450

[96] and that elderly females have higher skeletal muscle FOXO3 expression than younger 451 females [97] . While FOXO3 skeletal muscle gene expression differs between males and 452 females, it seems that the direction is opposite in young and old individuals, which emphasizes 453 the caution that should be used when interpreting sex differences across a large age range of 454 individuals. The promoter, 1 st exon, and gene body of GGT7 were hypomethylated in males 455 and males had higher GGT7 expression. GGT7 is highly expressed in skeletal muscle and 456 metabolises glutathione, which is a ubiquitous "master antioxidant" that contributes to cellular 457 homeostasis. Efficient glutathione synthesis and high levels of glutathione-dependent enzymes 458 are characteristic features of healthy skeletal muscle and are also involved in muscle 459 contraction regulation [98] . 460

In conclusion, we showed that the DNA methylation of hundreds of genes differs 461 between male and female human skeletal muscle. We uncovered important biological factors 462 underlying sex-specific skeletal muscle DNA methylation. Integration of the DNA methylome 463 and transcriptome, as well as gene expression validation, identify sex-specific genes associated 464 with muscle metabolism and function. Uncovering the molecular basis of sex differences across 465 different tissues will aid in the characterization of muscle phenotypes in health and disease. 466

The effects of other upstream drivers on sex differences in the muscle methylome, such as non-467 muscle cell type, the XY chromosomes, and genetic variants still need to be explored. 468 and biological confounding, such as population substructure, batch effects, and cellular 511 heterogeneity [105] . The inflation factor is higher when the expected number of true 512 associations is high (as it is for age); it is also greater for studies with higher statistical power 513

[1]. The results were consistent with the inflation factors and biases reported in an EWAS of 514 age in blood [1] . Results from the independent EWAS were combined using an inverse variance 515 weighted meta-analysis with METAL [106]. We used METAL since it does not require all 516 DNA methylation datasets to include every CpG site on the HumanMethylation arrays. For 517 robustness, we only included CpGs present in at least 2 of the 3 cohorts (633,645 CpGs). We 518 used a fixed effects (as opposed to random effects) meta-analysis, assuming one true effect size 519 of sex on DNA methylation, which is shared by all the included studies. Nevertheless, 520

Cochran's Q-test for heterogeneity was performed to test whether effect sizes were 521 homogeneous between studies (a heterogeneity index (I2 ) > 50% reflects heterogeneity 522 between studies). 523

To identify DMPs, we used linear models as implemented in the limma package in R 524

[107], using the participants' ID as a blocking variable to account for the repeated measures 525 design (for twin (GSE38291) and duplicate samples (Gene SMART), using 526 DuplicateCorrelation). The main sources of variability in methylation varied depending on the 527 cohort and were adjusted for in the linear model accordingly. For the Gene SMART study, we 528 adjusted the linear models for age, batch (2017 vs 2019), sex, timepoint and the interaction of 529 sex and timepoint (before and after four weeks of high-intensity interval training). For the 530 FUSION study, we adjusted the linear models for age, sex, BMI, smoking status, and OGTT 531 status. For the GSE38291 study, we adjusted the linear models for age, sex, and diabetes status. 532

All results were adjusted for multiple testing using the Benjamini and Hochberg correction 533

[108] and all CpGs showing an FDR < 0.005 were considered significant [109] . DMRs were 534 identified using the DMRcate package [110] . DMRs with Stouffer, Fisher, and harmonic mean 535 of the individual component FDRs (HMFDR) statistics < 0.005 were deemed significant. Effect 536 sizes are reported as mean differences in DNA methylation (%) between the sexes. 537

To assess whether fibre type proportions differ between males and females in the Gene 538 SMART and FUSION cohorts we used a beta regression model [111] using the betareg 539 package in R. We then included type I fibre ratio as a covariate in the linear models. To 540 investigate whether circulating hormones is associated with DNA methylation of the sex-541 DMPs, we included each sex hormone as a covariate in a linear model, separately. One male 542 from the Gene SMART cohort had missing circulating hormone values; and one two females 543 from the Gene SMART cohort had missing type I fibre proportions. For these three individuals 544 missing values were imputed with the mice package in R [112] . 

Myosin heavy chain is currently the best available marker for fibre typing [117] . Gene 556 SMART muscle sections were frozen in optimum-cutting temperature (OCT) medium by 557 holding the sample with OCT in liquid-nitrogen cooled Isopentane until frozen. Samples were 558 stored in -80°C until they were sectioned at 8µM with a cryostat. The IHC protocol was 559 performed as is described elsewhere [118] . Briefly, sections were blocked in 4% goat serum 560 (Invitrogen). Primary antibodies BA-F8 (MHCI), BF-35 (MHCIIA), and 6H1 (MHCIIX) were 561 purchased from DSHB, Iowa. Secondary antibodies goat anti-mouse IgG2b 350, goat anti-562 mouse IgG1 488, and goat anti-mouse IgM 555 were purchased from Invitrogen. Some samples 563

were fixed in paraformaldehyde for other analyses and for those samples an antigen-retrieval 564 protocol consisting of a 10 min incubation at 50° C of Proteinase K diluted in milliQ (1:1000) 565 and subsequent 1 min washes was performed before the IHC. Imaging was performed on the 566 Olympus BX51. 567

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

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

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

The hormone assays were completed in the accredited pathology laboratory at Monash 598 Health, Australia. Estradiol (E2) and Progesterone assays are competitive binding 599 immunoenzymatic assays performed on the Unicel DXI 800 system (Beckman Coulter). FSH 600 assay is based on Microparticle Enzyme Immunoassay (MEIA) and is carried out on the Unicel 601 DXI 800 system (Beckman Coulter). The LH and sex-hormone binding globulin (SHBG) 602 assays were performed using a sequential two-step immunoenzymatic ("sandwich") assay 603 carried out on a Unicel DXI 800 (Beckman Coulter). Testosterone was measured using the 604 HPLC-tandem mass spectrometry method using a liquid sample extraction (AB Sciex Triple 605

Quad 5500 liquid chromatography-tandem mass spectrometry). Free androgen index was 606 calculated as (total testosterone × 100)/SHBG. Free testosterone (fT) was calculated by the 607 Södergard fT calculation (36). 608 609

The Genotype-Tissue Expression (GTEx) Project sex-biased data was downloaded 611 from the GTEx Portal on 08/26/2020 and filtered for skeletal muscle samples. The enrichment 612 of DMG for GTEx DEGs was done by supplying the list of sex-biased genes to the gsameth 613 function in the missMethyl R package [115, 116] , which performs a hypergeometric test, taking 614 into account biases due to the number of CpG sites per gene and the number of genes per probe 615 on the EPIC array. Therefore, caution should be taken when interpreting the number of DMPs 616 reported per DMG. The analysis for direction of correlation between DNA methylation and 617 gene expression was performed by randomly shuffling DNA methylation effect sizes and 618 performing 10,000 permutations to assess how often a negative correlation occurs. This 619 analysis was performed for both GTEx and FUSION transcriptome data and yielded similar 620 results; data presented reflect results from the integration of differential methylation with 621 differential GTEx expression. Significance reported for GTEx sex-biased genes is represented 622 as the local false sign rate (lfsr) which is analogous to FDR [120] . 

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

The authors declare that they have no conflict of interest. 673

The dataset generated and analysed during the current study are available in the GEO 675 repository, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE171140. 676

Not applicable. (DMPs) with type I fibre proportion across all CpGs conducted with a meta-analysis of males and females, separately, from the Gene SMART and FUSION cohorts. Volcano plot of DNA methylation changes per percent increase in type I fibre content (expressed at percentage of beta value). Each point represents a tested CpG (665,904 in total) and those that appear in color are DMPs at a false discovery rate < 0.005; red DMPs are hypermethylated in type I fibres; blue DMPs are hypomethylated in type I fibres. The x-axis represents the amount of DNA methylation difference with increasing type I fibre content and the y-axis represents statistical significance (higher = more significant). (B) Percent of type I fibres in Gene SMART female and male muscle samples as determined by immunohistochemistry; p-value = 0.0001; 59% in females versus 47% in males. (C) Percent of type I fibres in FUSION female and male muscle samples as determined by RNA-seq; p-value = 2.06 x 10-7; 48.6% in females versus 39.9% in males. Each point is one of the 56,813 differentially methylated positions (DMPs) discovered in the full meta-analysis. To compare results from the full and partial meta-analysis we plotted the effect size (B value percentages) in the full metaanalysis (x-axis) against the effect size of the partial metaanalysis (y-axis). To show whether DMPs remained significant in the partial meta-analysis, we coloured points according to the false discovery rate (FDR) in the partial meta-analysis. 

Controlling bias and inflation in 764 epigenome-and transcriptome-wide association studies using the empirical null 765 distribution

A map of direct TF-DNA interactions in the human genome. 767 Nucleic acids research

The role of sex in the genomics of 769 human complex traits

Sex differences in immune responses that underlie COVID-19 771 disease outcomes

Sex bias and omission in neuroscience research is influenced 773 by research model and journal, but not reported NIH funding. Frontiers in 774 neuroendocrinology

Evidence for sex-specific genetic 776 architectures across a spectrum of human complex traits

The impact of sex on gene expression across human tissues

Sex differences in obesity, lipid metabolism, and 781 inflammation-A role for the sex chromosomes? Molecular metabolism

Y chromosome's roles in sex differences in disease

What a difference an X or Y makes: sex 786 chromosomes, gene dose, and epigenetics in sexual differentiation

Tissue-specific pioneer factors associate with androgen 792 receptor cistromes and transcription programs. The EMBO journal

Sex-specific differences in lipid and 795 glucose metabolism

Sex Differences in Gene Expression and Regulatory 797 Networks across 29 Human Tissues

Sex differences in global mRNA content of human skeletal muscle. 799 PLoS One

Sex differences in exercise metabolism and the role of 17-beta 801 estradiol

Genetic and epigenetic sex-specific adaptations to endurance 803 exercise

Sex-Based Differences in Skeletal 805

Muscle Kinetics and Fiber-Type Composition. Physiology (Bethesda)

Sex Differences in Human Fatigability: Mechanisms and Insight to

Sex differences in the genome-wide DNA methylation pattern and 810 impact on gene expression, microRNA levels and insulin secretion in human 811 pancreatic islets

Characterization of whole-genome autosomal differences of DNA 813 methylation between men and women

A study of the influence of sex on genome wide methylation

Sex influences DNA methylation and gene expression in human 817 skeletal muscle myoblasts and myotubes

The human skeletal muscle transcriptome: sex differences, 820 alternative splicing, and tissue homogeneity assessed with RNA sequencing

Sex-Related Differences in Gene Expression 823 in Human Skeletal Muscle

The landscape of sex-differential transcriptome and 825 its consequent selection in human adults

The effect of sex hormones on skeletal 827 muscle adaptation in females

Determinants of Receptor-and Tissue-829

Specific Actions in Androgen Signaling

Sex-specific differences in lipid and 831 glucose metabolism. Front Endocrinol (Lausanne)

JASPAR 2020: update of the open-access database of transcription 833 factor binding profiles

Sex Differences in Gene Expression and Regulatory 835 Networks across 29 Human Tissues

Identification of novel androgen-regulated pathways and mRNA 837 isoforms through genome-wide exon-specific profiling of the LNCaP transcriptome. 838 PLoS One

Identification of androgen response elements in the insulin-like growth 840 factor I upstream promoter. Endocrinology

The impact of ERα action on muscle metabolism and insulin 842 sensitivity -Strong enough for a man, made for a woman

Expression of both oestrogen receptor alpha and beta in human 845 skeletal muscle tissue

Effect of resistance exercise on muscle steroid receptor protein 847 content in strength-trained men and women

Large Scale Gene Expression Meta-Analysis Reveals Tissue-849

Sex-Biased Gene Expression in Humans

Strength and cross-sectional area of human 851 skeletal muscle. The Journal of physiology

Changes in skeletal muscle in males and females following 853 endurance training. Canadian journal of physiology and pharmacology

DNA methylation assessment from human slow-and fast-twitch 856 skeletal muscle fibers

Skeletal muscle fiber type: influence on contractile and 858 metabolic properties

Meta-analysis of genome-wide DNA methylation and integrative 860 OMICs in human skeletal muscle. bioRxiv

DNA methylation in the pathogenesis of type 2 diabetes in 862 humans

The Role of DNA Methylation in Mammalian Epigenetics

Integrative analysis of 111 reference human epigenomes

Systematic chromatin state comparison of epigenomes 868 associated with diverse properties including sex and tissue type

Dynamic regulation of transcriptional states by 871 chromatin and transcription factors

UniBind: maps of high-confidence direct TF-DNA interactions 873 across nine species. bioRxiv

Identification of methylation haplotype blocks aids in deconvolution of 875 heterogeneous tissue samples and tumor tissue-of-origin mapping from plasma DNA. 876 Nature genetics

Discovering high-resolution patterns of differential DNA 878 methylation that correlate with gene expression changes

Modeling complex patterns 881 of differential DNA methylation that associate with gene expression changes. Nucleic 882 acids research

Gender differences in strength and muscle fiber characteristics

A novel atlas of gene expression in human skeletal muscle reveals 887 molecular changes associated with aging. Skeletal muscle

Meta-analysis of human methylation data for evidence of sex-889 specific autosomal patterns

Sex differences in DNA methylation assessed by 450 K BeadChip in 891 newborns

Sex differences of leukocytes DNA methylation adjusted for 893 estimated cellular proportions. Biology of sex differences

Cell-type-specific disturbance of DNA methylation pattern: a chance to get 895 more benefit from and to minimize cohorts for epigenome-wide association studies

DNA methylation landscapes: 898 provocative insights from epigenomics

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

Single-cell transcriptional profiles in human skeletal muscle. 904 Scientific reports

A reference single-cell transcriptomic atlas of human skeletal 906 muscle tissue reveals bifurcated muscle stem cell populations. Skeletal muscle

Stem Cell-Based and Tissue 909 Engineering Approaches for Skeletal Muscle Repair

DNA methylation landscapes: provocative insights from 912 epigenomics

A genome-wide study of DNA methylation patterns and gene 914 expression levels in multiple human and chimpanzee tissues

Tissue-specific DNA methylation is conserved across human, mouse, 917 and rat, and driven by primary sequence conservation

Differential influence of peripheral and systemic sex steroids on 920 skeletal muscle quality in pre-and postmenopausal women

Intramuscular sex steroid hormones are associated with skeletal 923 muscle strength and power in women with different hormonal status

The effect of sex hormones on skeletal 926 muscle adaptation in females: Influence of sex hormones on female muscle 927 physiology

Androgen receptor in human skeletal muscle and cultured 929 muscle satellite cells: up-regulation by androgen treatment

Diminished androgen and estrogen receptors and aromatase levels 932 in hypogonadal diabetic men: reversal with testosterone

Time trajectories in the transcriptomic response to exercise-a meta-935 analysis

Let's Talk about Placental Sex

Mechanisms That Drive Female-and Male-Specific Fetal Growth and Developmental 938 Outcomes

A systematic review and meta-analysis examining whether 940 changing ovarian sex steroid hormone levels influence cerebrovascular function. 941 Frontiers in physiology

Risk for Venous Thromboembolism in Transgender Patients 943 Undergoing Cross-Sex Hormone Treatment: A Systematic Review. The Journal of 944 Sexual Medicine

Hormone replacement therapy improves contractile function and 946 myonuclear organization of single muscle fibres from postmenopausal monozygotic 947 female twin pairs

Gender differences in skeletal muscle substrate 949 metabolism-molecular mechanisms and insulin sensitivity. Frontiers in 950 endocrinology

Sexual Dimorphism in Mammalian Autosomal Gene Regulation 952 Is Determined Not Only by Sry but by Sex Chromosome Complement As Well

Minireview: Sex Differences in Adult and Developing Brains: 955 Compensation, Compensation, Compensation. Endocrinology

X chromosome dosage and presence of SRY shape sex-specific 958 differences in DNA methylation at an autosomal region in human cells

Requirement for Xist in X chromosome inactivation

X-inactivation profile reveals extensive variability in X-963 linked gene expression in females

The effect of X-linked dosage compensation on complex trait 965 variation

Landscape of X chromosome inactivation across human tissues

Multilocus loss of DNA methylation in individuals with 969 mutations in the histone H3 lysine 4 demethylase KDM5C

Widespread DNA hypomethylation and differential gene expression 972 in Turner syndrome

DNA methylation signature in peripheral blood reveals distinct 974 characteristics of human X chromosome numerical aberrations

Genetic impacts on DNA methylation: research findings 977 and future perspectives

Promoter DNA hypermethylation and paradoxical gene activation. 979 Trends in cancer

HDAC4 as a potential therapeutic target in neurodegenerative 981 diseases: a summary of recent achievements. Frontiers in cellular neuroscience

HDAC4 preserves skeletal muscle structure following long-term 984 denervation by mediating distinct cellular responses

Exercise-induced histone modifications in human skeletal muscle

Glucose sensing by skeletal myocytes couples nutrient signaling to 989 systemic homeostasis. Molecular cell

Exercise-associated DNA methylation change in skeletal muscle and 991 the importance of imprinted genes: a bioinformatics meta-analysis

Recent advances in understanding the role of FOXO3. 994 F1000Research

γ-Glutamyl transferase 7 is a novel regulator of glioblastoma growth. 996 BMC Cancer

Effect of sex on the acute skeletal muscle response to sprint 998 interval exercise. Experimental physiology

Skeletal muscle gene expression profiles in 20-29 year old and 65-1000 71 year old women

Glutathione and Nitric Oxide: Key Team Players in Use and 1002 Disuse of Skeletal Muscle

Genome-wide analysis of DNA methylation differences in 1006 muscle and fat from monozygotic twins discordant for type 2 diabetes

ChAMP: updated methylation analysis pipeline for Illumina 1009 BeadChips

Critical evaluation of the Illumina MethylationEPIC BeadChip 1011 microarray for whole-genome DNA methylation profiling

Cross-reactive DNA microarray probes lead to false discovery of 1014 autosomal sex-associated DNA methylation. The American Journal of Human 1015 Genetics

Abecasis, METAL: fast and efficient meta-analysis of 1020 genomewide association scans

Linear Models for Microarray Data

Controlling the false discovery rate: a practical and 1025 powerful approach to multiple testing

Redefine statistical significance

De novo identification of differentially methylated regions in the 1030 human genome

Beta regression for modelling rates and proportions

Multivariate imputation by chained 1034 equations in RJ Stat. Softw

Comprehensive characterization, annotation and 1036 innovative use of Infinium DNA methylation BeadChip probes. Nucleic acids research

missMethyl: an R package for analyzing 1041 data from Illumina's HumanMethylation450 platform

Gene set enrichment analysis for 1044 genome-wide DNA methylation data. bioRxiv

Fiber types in mammalian skeletal muscles

Rapid determination of myosin heavy chain 1048 expression in rat, mouse, and human skeletal muscle using multicolor 1049 immunofluorescence analysis

Flexible statistical methods for estimating and testing effects in 1054 genomic studies with multiple conditions

Supplementary