Comparing the epigenetic landscape in myonuclei purified with a PCM1 antibody from a fast/glycolytic and a slow/oxidative muscle

Muscle cells have different phenotypes adapted to different usage, and can be grossly divided into fast/glycolytic and slow/oxidative types. While most muscles contain a mixture of such fiber types, we aimed at providing a genome-wide analysis of the epigenetic landscape by ChIP-Seq in two muscle extremes, the fast/glycolytic extensor digitorum longus (EDL) and slow/oxidative soleus muscles. Muscle is a heterogeneous tissue where up to 60% of the nuclei can be of a different origin. Since cellular homogeneity is critical in epigenome-wide association studies we developed a new method for purifying skeletal muscle nuclei from whole tissue, based on the nuclear envelope protein Pericentriolar material 1 (PCM1) being a specific marker for myonuclei. Using antibody labelling and a magnetic-assisted sorting approach, we were able to sort out myonuclei with 95% purity in muscles from mice, rats and humans. The sorting eliminated influence from the other cell types in the tissue and improved the myo-specific signal. A genome-wide comparison of the epigenetic landscape in EDL and soleus reflected the differences in the functional properties of the two muscles, and revealed distinct regulatory programs involving distal enhancers, including a glycolytic super-enhancer in the EDL. The two muscles were also regulated by different sets of transcription factors; e.g. in soleus, binding sites for MEF2C, NFATC2 and PPARA were enriched, while in EDL MYOD1 and SIX1 binding sites were found to be overrepresented. In addition, more novel transcription factors for muscle regulation such as members of the MAF family, ZFX and ZBTB14 were identified.


Introduction
The phenotype of skeletal muscle fibers differs as an adaption to different tasks. Some fibers have short twitches and rapid shortening velocities, but low endurance. Such fibers are used for short, but explosive external work (sprinting, ballistic and corrective movements etc.), but are easily fatigued and not very energy efficient. Other fibers have slow twitches and shortening velocities, but are fatigue resistant and have a low energy expenditure. Such fibers are used for tasks like keeping body and limb posture [1][2][3].
Fibers are generally classified into "fiber types" related to the predominant myosin heavy chain (MyHC) isoenzyme expressed in the cell. MyHC is the strongest molecular determinant of shortening speed [2,4]. In rodents, there are four different MyHC genes expressed in adult limb muscles, namely Myosin Heavy Chain 7 (Myh7), Myh2, Myh1 and Myh4 giving rise to the fiber types 1 (slowest), 2A, 2X and 2B (fastest), respectively [2]. Partly coupled to the MyHC fiber type, the fibers display a spectrum of metabolic properties, from highly oxidative mitochondria-rich cells to cells that are mainly glycolytic [2,5]. At the molecular level, the different fiber types vary in isoform expression of various proteins such as calcium pumps, oxygen related proteins and also sarcomere components other than MyHC [2,3,5].
Skeletal muscle is the most important metabolic organ in the body, and it has been known for 40 years that differences in muscle fiber type composition are a strong individual predictor for developing metabolic syndrome [6-10]. Metabolic syndrome is a cluster of conditions increasing the risk of heart disease, stroke and type 2 diabetes. These conditions include increased blood pressure, high blood sugar levels, excess body fat around the waist, and abnormal cholesterol or triglyceride levels. Metabolic syndrome is on the rise, and in several countries, the prevalence is now over ¼ of the population [11]. Epigenetic mechanisms like histone modifications and chromatin structure have been suggested to play an important role in the development of and predisposition for metabolic syndrome [12], but data supporting this is currently scarce.
Muscle fibers are post-mitotic and represent an interesting example of a balance between long term phenotypic stability, yet malleability. Thus, phenotype can be changed by altering the pattern of electrical activity elicited in their sarcolemma by the motor neurons or electrical stimulation [3] and also by the mechanical stress created by contraction [13]. These two external factors are the major mechanistic foundations for the effects of exercise on muscle. However, the adult phenotype also depends on the embryonal cell lineage and this origin limits the adaptive ranges of training effects [14]. The transition process has been studied in detail in rats stimulated with electrical patterns mimicking extreme, but well-defined training [3]. It appears that some traits require very long-term treatment in order to be altered, e.g. type 2 to type 1 transformations might take more than 3 months [15]. In contrast, changes within type 2, i.e. 2B>2A>2X [16] and changes in metabolic and calcium sequestering enzymes can be altered more rapidly [17]. Epigenetic mechanisms may play an important role in shaping such differences [18][19][20][21].
Recently the existence of a long-term cellular memory was demonstrated in muscle cells related to hypertrophy [22,23], and in addition to a permanently elevated number of myonuclei [24,25], epigenetic mechanisms might be involved [26]. Epigenetic mechanisms might also be related to the observation that some individuals seem to have less malleable muscles than others, i.e. exercise resistance [27][28][29].
Since the overall chromatin environment and the modifications of the histones represent overreaching gene regulatory mechanisms with the potential of long-lasting stability [26,30,31], we set out to compare the chromatin environment in an extremely fast/glycolytic muscle (the m. extensor digitorum longus, EDL) and an extremely slow/oxidative muscle (the m. soleus) in mice.
The majority of studies on epigenetics are on tissue culture cells, and less has been done on tissue homogenates where results seem harder to interpret. In muscle, electron microscopy [32] and a specific marker for myonuclei [32,33] have revealed that only approximately 40-50% of the nuclei found in muscle tissue are myonuclei, and since accounting for cellular heterogeneity is critical in epigenome-wide association studies [34][35][36][37], we aimed at sorting out the myonuclei proper for our epigenetic analysis.
We recently reported that the nuclear envelope protein Pericentriolar material 1 (PCM1) can be used as a specific light microscopy marker to discern the skeletal muscle myonuclei in both rodents and humans [33], and we show here that this marker can be used to sort myonuclei to >95% purity from mice, rats as well as human samples and be used for subsequent epigenetic analyses. Furthermore, we show that the purification is necessary for the epigenetic landscape faithfully to reflect known features of muscle function, and the results indicate that purification should be used in studies aimed at elucidating the role of epigenetic mechanisms in muscle differentiation and plasticity.

Results
Purification of myonuclei allows a genome-wide analysis of muscle fiber specific chromatin PCM1 labelling is specific for myonuclei isolated from the muscle tissue. To remove interference from non-myofiber nuclei in the tissue, we took advantage of an antibody against Pericentriolar material 1 (PCM1) which we have previously reported to selectively label myonuclei in skeletal muscle fibers and not satellite cells or stromal cells [33,38] (S1A and S1B Fig).
To further prove the specificity of PCM1 for muscle myonuclei, we used a transgenic mouse model expressing Histone-2B coupled to GFP (H2B-GFP) under the control of the skeletal muscle Actin Alpha 1 (ACTA1) promoter, giving specific expression of the fusion construct in skeletal muscle fiber nuclei across fiber types [39]. Single fiber analysis showed that H2B-GFP and PCM1 positive nuclei invariably colocalized (Figs 1A and 1B and S1C). In some cases, DAPI labelled nuclei that were negative for both H2B-GFP and PCM1 were observed (arrows in Fig 1A), probably representing non-myonuclei adhering to the fibers after isolation.
PCM1 can be used to isolate myonuclei for subsequent epigenetic analysis in rodents and humans. Flow cytometry analysis of nuclei from skeletal muscle displayed a 98.0% overlap between GFP and PCM1 labelling of the nuclei in the transgenic mouse line (Fig 2A-2F).
Analysis by flow cytometry of the nuclear composition in the tibialis anterior (TA), EDL and soleus showed that the number of myonuclei varied from 40%-60% with the lowest fraction of myonuclei in the slow soleus (Fig 2G-2K). This confirms that without sorting, roughly half of the nuclei isolated from muscle homogenates are from other cell types than muscle fibers and will confound a genome-wide analysis aimed at this cell type.
To isolate the myonuclear fraction from the nuclei originating from the other cell types in the tissue, we used a magnetic-assisted sorting approach (Fig 2G-2J). Using this method, we were able to isolate the myonuclei with a very high purity averaging 95% for the three different muscle types (Fig 2G-2L). The same high sorting purity was obtained when the three same muscle types from rat and vastus lateralis biopsies from humans were analyzed (S2 Fig).
Purification of nuclei is required for a valid genome-wide analysis of the myo-specific epigenome. To investigate the importance of nuclear purification for epigenetic analysis of muscle cells, we performed chromatin immunoprecipitation on the PCM1 purified nuclei, followed by next-generation sequencing (ChIP-Seq) in soleus and EDL. An antibody against Histone H3 acetylated at Lysine 27 (H3K27ac) was used, since its enrichment at promoter regions follows the transcriptional activity [40,41], and compared with previously published H3K27ac data from whole soleus and EDL muscles [42]. In the non-sorted material, 12550 and 12243 promoters with H3K27ac peaks were identified in soleus and EDL respectively, and 27% and 14% of these did not have enrichment after purification (Fig 3A-3F). For both soleus and EDL, gene ontology analysis of the promoters with signal only in whole tissue showed enrichment for the immune system and angiogenesis (Fig 3D and 3E, and S1A and S1B Table). Extending the analysis to the full epigenomes for whole muscle tissue and the purified myonuclei, showed similar results with enrichment for angiogenesis and blood related functions for both muscles (S3A- S3D Fig and S2A-S2D Table).   The effects of purification were well illustrated by the observation that genes used to identify non-muscle cells were found to be enriched for H3K27ac in non-purified nuclei, but not after the purification, e.g. Fatty Acid Binding Protein 4 (FABP4) and Endothelial Cell-Specific Molecule 2 (ECSCR) which are both used as marks for endothelial cells [43][44][45][46][47] 65] was enriched in the total nuclei from the whole tissue, but not after purification of the myonuclei.
The loci for genes classically known to be a part of the skeletal muscle cell, such as the skeletal muscle specific version of actin (ACTA1) and Titin (TTN), displayed enrichment in both the full tissue and PCM1-sorted population ( Fig 3G).

The epigenetic landscape in purified myonuclei from EDL and soleus reflects the differences in functional properties
To explore how the epigenome differed between the slow/oxidative soleus and fast/glycolytic EDL muscle cells, we focused on differentially enriched (DE) peaks between the two muscles using the histone modifications H3K4me3, a mark associated with promoter regions of expressed genes and H3K27ac which is also associated with promoters of expressed genes as well as distal regulatory enhancers where the chromatin is in an open and accessible configuration [66][67][68]. A fold change (FC) >1.5 and a false discovery rate (FDR) < 1x10 -2 and 1x10 -7 for H3K4me3 and H3K27ac, respectively. We identified 719 H3K4me3 and 5309 H3K27ac DE peaks, corresponding to 7,6% and 22,2% of the sites, for each of these two histone modifications (Fig 4A and 4B). For H3K4me3, the majority of sites were localized around the promoter regions, while for H3K27ac, most of the sites were found at distal regulatory enhancer regions (S4A- S4F Fig). The higher number of DE peaks for H3K27ac suggest that distal regulatory enhancers play an important role in regulating the differences in function and phenotype, in line with previous observations of the epigenetic landscape between cell types from the same family [69,70].
The epigenome reflects fiber type specific differences in contractile properties. The most fundamental way of classifying muscle fiber types has been the particular MyHC gene expressed in the fiber, but the sarcomeres differ with respect to other contractile proteins as well [2,71]. These features were reflected in the epigenome ( Fig 4C). Thus, the slowest MyHC isoform Myh7 (coding for MyHC-1) on the mouse chromosome 14 had a larger enrichment in soleus compared to EDL, while the opposite was true for the fastest form Myh4 (MyHC-2B) that is clustered with the other type 2 heavy chains on the mouse chromosome 11 (Fig 4C). When considering all the four limb muscle MyHCs, the H3K27ac enrichment roughly resembled the histological fiber type distribution of the two muscles that we have found previously [72] (Fig 4E). Similar differences in the epigenetic enrichment were found at the loci coding for the other sarcomeric proteins, such as the myosin light chains (Myl1-3) (Fig 4D). The fast specific for whole muscle tissue (grey). See S1A and S1B Table for

PLOS GENETICS
isoform of the gene coding for the contraction regulator Myosin Binding Protein C2 (Mybpc2) was enriched in the EDL, while paradoxically, the gene coding for the slow isoform Myosin Binding Protein C1 (Mybpc1), was not enriched in the soleus, however, this is in agreement with RNA and protein levels for the isoform [73][74][75].
The epigenome reflects fiber type specific differences in twitch duration. One of the major differences between fiber types is the twitch duration, hence the EDL is a so called fasttwitch muscle and soleus a slow-twitch muscle. A major determinant of the twitch speed is the handling characteristics of Ca 2+ [2,3], and the expression level of relevant isoenzymes differs between different types of muscles [17,50,75,76]. Our data showed that this was also reflected at the epigenetic level, such as for the ATPase Sarcoplasmic/Endoplasmic Reticulum Ca2 + Transporting (Atp2a)-1 and -2 (SERCA1 and -2), Casq1 and -2, coding for the two calcium buffer proteins, Calsequestrin 1 and -2, responsible for regulating the calcium levels in fast and slow muscle fibers respectively, and Pvalb, coding for the calcium binding protein Parvalbumin, highly specific for fast fibers [17] (Fig 5).
For the genes coding for contractile and calcium handling proteins, the differences in the epigenome were not only restricted to the promoter regions, but also the inter-and intragenic areas with H3K27ac enrichment, such as the enhancer located inside the MyHC cluster at chromosome 11 between the long non-coding RNA Linc-Myh and Myh2 (MyHC-2A) promoting fast glycolytic phenotype [77], and upstream enhancer regions at the slow isoforms of Myh7 (MyHC-1) and Atp2a2 (SERCA2) genes [78] (Figs 4B and 5). Furthermore, we also identified the recently discovered fast phenotype specific enhancer located in the locus coding for the Pvalb (Parvalbumin) gene [42]. This indicates the complexity and detailed regulation of the genes involved in contraction and calcium handling in the different muscle types.
The epigenome reflects fiber type specific differences in metabolic properties. Ontology analysis of the H3K27ac DE peaks showed that the regions are associated with genes involved in defining the metabolic properties for the two muscle types, such as lipid metabolism and mitochondria for soleus and muscle contraction and metabolism of hexose (e.g. sugar) for EDL (Fig 6A and 6B, and S3A and S3B Table). Interestingly, a similar analysis carried out on whole tissue identified fewer DE regions and only non-specific general functions (S4G and S4H Fig, and S3C and S3D Table), again, underscoring the importance of sorting out the relevant nuclei from heterogeneous tissues.
To gain further insight into how the epigenome adds to differences between the two types of muscles, we examined the differences in enrichment using the dynamic H3K27ac mark. The analysis identified differences in H3K27ac enrichment, both at promoters and distal regulatory regions of genes, such as Peroxisome proliferator-activated receptors (PPAR)-and the Fatty acid degradation pathways in soleus and the glycolysis and mitogen-activated protein kinase (MAPK) pathways in EDL (S4 Table). Detailed analysis of the epigenetic landscape for two of the primary energy generating pathways, the fatty acid oxidation in soleus and glycolysis in EDL, identified a complex epigenetic environment with several promoters and enhancers with a different enrichment for the genes involved in the pathways (Fig 6C-6E). For both the fatty acid oxidation pathway and glycolysis, none of the involved genes were associated with DE regions for the opposite muscle type, e.g. glycolysis enriched in soleus and fatty acid oxidation in EDL.

Differences in the epigenome indicates that the phenotype of fast/glycolytic and slow/oxidative muscles are determined by different regulatory networks
We sought to unravel the gene regulatory programs leading to the fast/glycolytic and slow/ oxidative phenotypes respectively, by using the myonuclei-specific DE H3K27ac regions in combination with DNase I hypersensitivity data for whole skeletal muscle [79,80], which with a high resolution identifies where transcription factors and other proteins are bound to the chromatin [81]. Motif enrichment analysis of the DE regions identified two distinct groups of transcriptional regulators for the two muscles. In soleus, the SRY-Box Transcription Factor (SOX) and Nuclear Factor of Activated T Cells (NFAT) motifs were the most overrepresented, while in EDL, it was the E-boxes and Sine Oculis Homeobox Homolog (SIX) motifs (S5 Table).
To consolidate and specify the analysis, we used the transcription factors from the motif enrichment analysis with the most significantly enriched promoter region and high-quality motifs for each muscle. This revealed two distinct groups of transcription factors, one for each muscle type. For soleus, factors such as SOX6, NFATC2, Myocyte Enhancer Factor 2C (MEF2C), MAF BZIP Transcription Factor K (MAFK), Peroxisome Proliferator-Activated Receptor Alpha (PPARA) and Transcription Factor 7 Like (TCF7L2) were identified, while for EDL, factors such as Myogenic Differentiation1 (MYOD1), SIX1, MAFG, Zinc Finger and BTB Domain Containing 14 (ZBTB14), and Zinc Finger Protein X-Linked (ZFX) (S5 Table). The transcription factors were identified to have 4530 and 4371 predicted binding sites connected to 464 and 392 DE promoter regions for soleus and EDL, respectively. Functional analysis of the genes connected to the transcription factors revealed that they are associated with central differences in the functional phenotype between the two muscles (S6 Table). Combining the information about factors, genes and functions revealed distinct regulatory networks for the two muscles (Fig 7). Though different, they also share similarities, such as regulation of muscle contraction, but with a different set of factors and genes involved (S6 Fig).

Discussion
We show that PCM1 selectively marks mature myonuclei in skeletal muscle, and that this marker can be used to isolate the myonuclei from fresh and frozen muscle tissue collected from mice, rats and humans. The method can be ubiquitously applied, also to experimental  Table  for full list. Ontologies were identified with GREAT [143] using closest promoter. (C-D) ChIP-Seq profiles of the myo-specific H3K4me3 and H3K27ac enrichment in soleus and EDL, for two loci with different enrichment, encoding Fatty acid translocase (Cd36/Fat) at chromosome 5 and Lactate dehydrogenase A (LDHA) at chromosome 7, respectively. (E-F) Schematic representation showing the genes encoding the enzymes involved in the fatty-acid oxidation and glycolysis, respectively. Genes with a differently enriched H3K27ac promoter colored in red and blue for soleus and EDL, respectively. Grey denoted no difference in enrichment. Differently enriched alternative promoters and enhancer regions are shown as nodes with relative distance to the primary promoter region. Color code denotes the FDR significance level of the differential enrichment between the two muscle types.
https://doi.org/10.1371/journal.pgen.1009907.g006 . The last row shows biological functions assigned to the genes. Edges between transcription factor and genes represent identification of predicted binding sites for the respective transcription factor and the gene it is connected to. Edges between genes and biological function denote the functional properties assigned to the gene. Transcription factor motifs were identified with AME [151] using high-quality motifs from human and mouse FDR < 1x10 -5 . For factors with similar motifs, the factor with the most significantly enriched promoter in the respective muscle type was used. Biological functions represent functional annotation clusters associated with the differential enriched genes (Clusters identified with David [144,145] using an enrichment score >2.2 for the clusters and a p-value and fold change for the ontologies at < 5x10 -2 and >3, respectively). Network visualized with Cytoscape. For the sake of clarity, only selected genes, with their most significant edge to each transcription factor, and terms are shown. See S5 conditions where transgenic markers are not feasible, such as human studies. The isolation method does not involve any enzymatic digestions or other treatments at elevated temperatures, thereby preserving the native state of the epigenome.
Propper isolation of the myonuclei is important since up to 60% of the cell-nuclei in muscle tissue can reside in other cell types such as satellite-, stroma-, and vascular cells, and in connective tissue [32,33]. Based on our FACS data; in mice and rats these other nuclei constituted about 40% in the fast EDL and TA and 20% in human biopsies from vastus lateralis. The corresponding numbers for the slow soleus were 60% in mice and 40% in rats. We don't know if this difference is a due to muscles or species variability.
The epigenetic signature of individual genes is highly cell-specific. Studies of other heterogeneous tissues have concluded that accounting for cellular heterogeneity is critical in epigenome-wide association studies [34][35][36]. Mixed cell populations are prone to increase the number of false negatives and positives [41,82]. This point has also been addressed in skeletal muscle by mechanically isolating single fibers preceding transcriptome [50] and methylome [83] analyses; reducing the signal arising from non-myogenic cells residing in the tissue.
Recent studies have mapped the genome-wide landscape of whole skeletal muscles, but without discriminating between the different cell types in the tissue [42,84]. Compared to analysis of the whole population of nuclei, we show that analysis of purified myonuclei gives a more detailed profile of the myo-specific regulatory environment, also confirmed by a recent analysis of single nuclei [85]. Purification of nuclei affected both the signal from specific gene loci as well as the general functional ontological analysis, which agreed better with the function of the muscle and the two muscle types.
We compared myonuclei from the soleus and EDL typifying a slow/oxidative and a fast/glycolytic muscle in mice, using the markers H3K4me3 and H3K27ac. The first is found in the promoter regions of expressed genes, the latter is associated with active promoters and distal regulatory enhancers where the chromatin is in an open and accessible configuration [66][67][68]. The two muscles had a large part of the epigenome in common, but important differences were found in connection to genes coding for proteins involved in contraction, calcium handling and metabolic pathways, reflecting differences at the transcriptional and protein level between different fiber types [5,50,75,86,87]. At the genome-wide level, we revealed differences in distal enhancers at the loci reflecting different regulatory networks in the fast and the slow muscle. The two MyHC loci provide interesting examples of differences in such regional regulation. In both mice and humans, the slow type 1 Myh7 (MyHC-1) gene is located on chromosome 14 head to tail with the cardiac isoform Myh6 (MyHC-1alpha), while the fast type 2 isoforms are clustered on a 300-600 kb segment in the order Myh2, Myh1 and Myh4, reflecting the order of increasing shortening velocity, (MyHC-2A-2X-2B) on the mouse chromosome 11 (in humans 17) [88]. Interestingly, by exercise, changes are relatively easily obtained within the different type 2 MyHCs of the gene cluster, while 1,2 transformation is less common, and cannot be achieved, except under extreme conditions [15].
In the fast locus, an intergenic SIX1 responsive enhancer is located between Myh3 (MyHCemb) and non-coding RNA Linc-Myh. It drives the expression of the fast myosin through higher order chromatin interactions and locks the muscle fiber in the fast-glycolytic phenotype [77]. Furthermore, the enhancer regulates the Linc-Myh, which might play a role in maintaining the fast phenotype, presumably through trans-induced mechanisms [77]. This area, which is one of the most differently enriched in our data for the glycolytic EDL, revealed several additional enhancer regions upstream and downstream of the Linc-Myh, further adding to the complexity of the locus. In muscle cell culture, this area has been identified to be part of a special group of enhancers known as super-enhancers [89]. This group of enhancers is characterized by spanning several kb and being important for defining the identity of the cell [90][91][92]. Such areas were among the most enriched in EDL, and it is possible that the chromatin configuration stabilizes the fast fiber type and provides an obstacle for transformation in the slow direction. As discussed in the introduction, it appears that some traits require very long-term treatment, while changes between different type 2 MyHC are readily altered by endurance exercise or high amounts of electrical stimulation. Type 2 to 1 transformations might take more than 3 months [15]. We hypothesize that this delay is related to the epigenetic landscape.
Several of the transcription factors identified in our analysis of regulatory networks are well known regulators of skeletal muscle phenotype and function, such as NFAT, SOX, SIX and MYOD1 [2,3,93]. Our enriched factors identified by isolating the myonuclei prior to epigenetic analysis are largely in agreement with a recent study using an alternative bioinformatic approach to identify and investigate the epigenetic regulatory environment in slow and fast muscles at the single nuclei level [85].
Other factors we identified, such as MAF transcription factors and ZBTB14, are just recently described to have a role in muscle regulation at the gene level [85,94]. The latter is identified to be involved in regulation of metabolism; a function also assigned to the factor in our regulatory network. More specifically, it is found to be a negative regulator of the myokines Interleukin 6 (IL6) and Leukemia inhibitory factor (LIF) by binding to the promoter region of the two genes [94].
We identified two additional factors that are less explored in skeletal muscle, namely the ZFX enriched in EDL and TCF7L2 in soleus, but both have previously been suggested to be a part of skeletal muscle specific regulatory networks [75,95]. On the molecular level, the ZFX is identified to facilitate transcription through binding in the promoter regions just downstream of the transcription start sites of target genes [96]. In the case of TCF7L2, it is found to be involved in cardiac hypertrophy through the Wnt/β-catenin signaling pathway, where it modulates the chromatin environment and regulates the expression of the slow Myh7 (MyHC-1) and the Wnt/β-catenin target, the Myc gene [97]. Interestingly, in our regulatory network, the factor is also associated with Myh7, hence pointing towards a conserved role in oxidative muscle tissues for TCF7L2 in regulating the slow-twitch myosin. To shed further light on the two factors´role in skeletal muscle biology with a focus on humans, we performed a meta-analysis of their expression in muscle using the MetaMEx [98] database which focuses on the response of human skeletal muscle to interactivity and exercise. The analysis showed that ZFX is upregulated by prolonged inactivity, while TCF7L2 is upregulated during acute aerobic exercise training. Furthermore, both factors are found to be expressed at the protein level in skeletal muscle [99,100]. For both factors, more research is needed to elucidate their function in skeletal muscle, such as if ZFX is having a similar role at promotor regions in enhancing the transcription of target genes during prolonged inactivity. While in the case of TCF7L2, whether it plays a role in Wnt/β-catenin induced hypertrophy in skeletal muscle [101] and whether the transcription factor regulates the slow-twitch myosin expression and Wnt/β-catenin targets such as Myc.
For the soleus, it is notable that we identified the transcription factors MEF2C, PPARA, SOX6 and NFATC2 as they are known to influence muscle in the slow/oxidative direction [2,[102][103][104][105]. MEF2C is previously identified to be important for slow fiber development and for their energy homeostasis [106,107]. MEF2C and its family members are known to interact with co-regulators that are capable of altering the epigenetic environment, such as histone acetyltransferase P300 [108,109] and lysine methyltransferase Myeloid/lymphoid or Mixed-Lineage Leukemia 4 (MLL4). In the latter case, the transcription factor, together with its coregulator, establishes an open chromatin environment at active enhancers that enforce an oxidative/slow muscle phenotype, e.g. by regulating sarcomeric slow genes such as Myh7 (MyHC-1) and Tnnc1 (Troponin C1) through enhancers located upstream of their respective promoters [78]. MEF2C is also participating in activity-related calcium regulation, both by being activated by upstream regulators, such as MAPK and calmodulin-dependent protein kinases (CaMK) [110][111][112], but also through regulation of the genes participating in the signaling cascades [113][114][115]. In the cardiac tissue, MEF2C is found to regulate the expression of the slow sarcoplasmic calcium pump, Atp2a2 (SERCA2), through an enhancer region [116], a connection also identified in our soleus network.
The transcription factors PPARA and Estrogen Related Receptor Beta (ESRRB), enriched in soleus, are mediating fatty acid degradation and oxidative metabolism [86,117,118]. Similar roles are known for additional members of their respective families [119,120]. In the case of PPARA, it is furthermore reported to control the expression of metabolic genes in response to changes in the environment through distal enhancers [121] and to be involved in regulation of contraction [122,123].
For EDL, it was notable that SIX1 and MYOD1 were enriched since there is evidence that they drive the phenotype in the fast/glycolytic direction [3,124,125]. The SIX1 transcription factor is reported to be an important regulator of key genes defining the fast-glycolytic phenotype, where it regulates genes mediating glycolysis (e.g. Aldoa (Aldolase A), Pfkm (Phosphofructokinase) and Eno3 (Enolase 3) and calcium homeostasis (Pvalb (Parvalbumin)) [126,127] in addition to regulating the non-coding RNA Linc-Myh and its enhancer region [127]. MYOD1 is found to be important for correct glycolytic fiber type specification and composition [125,128], and is thought to be a master regulator working as an organizer of enhancers [129][130][131]. In line with this, it is known to recruit histone modifiers and co-regulators, such as p300 and its homologue CBP, to target sites [108,132,133].
In a recent muscle cell differentiation study in cell culture, MYOD1 was identified to take part in the regulation of the super-enhancer inside the fast MyHC locus covering the Linc-Myh and its enhancer [134]. Pointing towards an interesting collaboration at the locus between MYOD1 as a master regulator and the more specific role of SIX1 in locking the expression of the fast myosin genes [77], thereby promoting the glycolytic phenotype.
In summary, we have shown that PCM1 can be used to isolate the myonuclei from the complex whole muscle tissue and that in order to get a clear and faithful representation of the whole epigenome, such sorting should be performed prior to analysis of nuclei from the complex muscle tissue. Comparison of the myo-specific epigenomes of the two muscle extremes soleus and EDL show, as expected, that the majority of the epigenome was common. The parts that differed were regulatory elements and genes related to contraction speed, twitch duration, and metabolism. In addition, our epigenetic analysis revealed that the two muscles have distinct regulatory networks associated with genes defining the disparate phenotype. All procedures with the ACTA1 rtTA ; TRE H2B-GFP mice were approved by the Cincinnati Children's Hospital Medical Center's Institutional Animal Care and Use Committee.

Ethics statement
Human vastus lateralis biopsies were a kind gift from Dr. Tim Snijders, Maastricht University. the samples were obtained as a part of a study approved by the Medical Ethical Committee of the Maastricht University Medical Centre+ (METC 15-3-003) which conformed to the standards for the use of human subjects in research as outlined in the most recent version of the Helsinki Declaration.

Animal procedures
Female NMRI mice (28-31 gram) and male Sprague Dawley rats weighing 210-260 g were used. Animals were kept at the animal facility at the Department of Biosciences or Department of Medicine, University of Oslo. Animals were housed with a 12 h light/dark cycle with ad libitum access to food and water. Before surgery, animals were sedated with 2% isoflurane (506949, Forene, Abbot) in the air. Following deep anesthesia, the hind limb muscles extensor digitorum longus (EDL), soleus and tibialis anterior (TA) were surgically removed, directly frozen in liquid nitrogen and transferred to cryotubes before being stored at -80˚C. Animals were sacrificed by cervical dislocation while under deep anesthesia. The research was conducted in accordance with the Norwegian Animal Welfare Act of 20th December 1974. The

Immunohistochemistry and imaging
For cross section tissue sectioning, immunohistochemical staining and subsequent image processing were performed as previously described [33].
For single fiber imaging; EDL muscles from the ACTA1 rtTA ; TRE H2B-GFP mice were collected and incubated in high-glucose DMEM (HyClone Laboratories) with 0.3% collagenase type I (Sigma-Aldrich) at 37˚C for 40 minutes, then washed with PBS. Muscles were gently triturated using a glass pipette to loosen digested myofibers until they shed from muscle and fixed overnight in 1% paraformaldehyde. Fibers were then put in staining buffer (0.6% Igepal-CA-630, 5% BSA and 1% goat serum in PBS) for 1h at room temperature, followed by staining with an antibody against PCM1 (1:1000, HPA023370, Sigma Aldrich) in staining buffer overnight at 4˚C with gentle agitation. Next day, fibers were washed (3x 30 min in staining buffer), incubated for 3 hours with secondary antibody (Abcam #150083) at 4˚C with gentle agitation, washed 3x 30 min in staining buffer and mounted with Fluoromount-G with DAPI (Southern-Biotech #0100-20). For Figs 1A and 1B and S1C, fibers were imaged on a Nikon A1 R confocal system through NIS-Elements AR software (ver.5.10.01) with a 60x oil immersion objective with NA of 1.4, using a Galvano scanner. For S1A and S1B Fig, fibers were imaged on a Nikon Eclipse Ti inverted microscope with a 20x air immersion objective.

Myonuclear isolation
All samples and buffers were kept on ice during the procedure. The muscles were removed from the -80˚C freezer and transferred to -20˚C and minced into pieces of 1-2 mm before being transferred to gentleMACS M-tubes (Miltenyi Biotec) equilibrated to the same temperature. Samples were transferred to ice and stored for 5 minutes before homogenized in 5 ml of lysis buffer (10 mM Tris-HCl pH 8.0 (cold adjusted), 5 mM CaCl 2 , 3 mM MgAc, 2 mM EDTA, 0,5 mM EGTA) on a gentleMACS Dissociator (Miltenyi Biotec) using the default homogenization program protein_01. All buffers were supplemented with 5 mM Na-butyrate (Sigma Aldrich), 5 mM PMSF (Sigma Aldrich) and 1x protease inhibitor (Complete Protease Inhibitor Cocktail, Roche) immediately before use. Following homogenization, the lysate was diluted in 5 ml lysis buffer containing 0.4% Triton-X100. Samples were mixed 10 times with a 1 ml pipette before being filtered through 100 μm and 30 μm strainers (Falcon, Sigma Aldrich) to remove large aggregates. The homogenates were transferred to a 15 ml tube and centrifuged at 1000 g for 5 minutes in a centrifuge with a swing-out rotor at 5˚C. The pellets were resuspended in 1 ml nuclei staining buffer (5% BSA wt/vol, 0.2% IGEPAL-CA630, 1 mM EDTA in Dulbecco's Phosphate Buffered Saline (DPBS) pH 7.4) [41] containing an antibody against PCM1 (1:1000, HPA023370, Sigma Aldrich). Samples were incubated using a tumbler (40 rpm) for 1 hour at 6˚C. Then centrifuged for 5 minutes at 800 g at 5˚C and the pellet was resuspended in 500 μl sorting buffer (1% BSA wt/vol, 2% skimmed milk powder (Sigma Aldrich), 1 mM EDTA in DPBS (pH 7.4) and spun down at 600 g at 5˚C for 5 minutes at slow acceleration settings. The pellets were resuspended in 100 μl beads buffer, consisting of 80 μl sorting buffer with 20 μl secondary Anti-rabbit IgG MicroBeads (Miltenyi Biotec), and incubated for 15 minutes at 4˚C. After incubation, the samples were centrifuged as above and washed once with Sorting buffer before resuspended in 1 ml of the same buffer and incubated for 5 minutes at a tumbler at 40 rpm at 6˚C. After incubation, the samples were carefully applied to a M column (Miltenyi Biotec), washed 3 times with 1 ml sorting buffer and eluted with 1 ml Sorting buffer after removal from the magnet. The eluates were applied to another M column, washed again and eluted in 1 ml elution buffer (DPBS with 1 mM EDTA). The efficiency of sorting was determined by staining aliquots before and after myonuclear enrichment.

Chromatin immunoprecipitation
For the chromatin immunoprecipitation, a modified version of [135,136] was used. All buffers were supplemented with 5 mM Na-butyrate (Sigma Aldrich), 5 mM PMSF (Sigma Aldrich) and 1x protease inhibitor (Complete Protease Inhibitor Cocktail, Roche) immediately before use. Nuclei, approximately 300 000, were fixated with 1% freshly prepared formaldehyde (Pierce Cat.No.28906) in crosslinking buffer (DPBS supplemented with 50 mM HEPES, pH 7.5, 1 mM EDTA, 0.5 mM EGTA) for 2 minutes at room temperature. Reaction was stopped by adding glycine to a final concentration of 125 mM. After 5 minutes of incubation at RT, nuclei were transferred to ice and washed twice with ice cold DPBS and lysed in 130 μl sonication buffer (50 mM Tris-HCl pH 8.0 (cold adjusted), 1 mM EDTA, 0.1% wt/vol SDS). Samples were sonicated for 8 minutes (30 sec on/off cycles) using a Bioruptor Pico (Diagenode) yielding average DNA fragments of 200-300 bp. Lysates were diluted in an equal volume of 2x ChIP buffer (20 mM Tris-HCl pH 8.0 (cold adjusted), 280 mM NaCl, 1 mM EDTA, 1 mM EGTA, 2% Triton X-100, 0.2% Na-deoxycholate, 1% Igepal CA630, 0.1, % SDS), and adjusted to 800 μl with ChIP buffer (10 mM Tris-HCl pH 7.5 (cold adjusted), 140 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, 1% Triton X-100, 0.1% Na-deoxycholate, 0.1% SDS). Samples were centrifuged at 12 000 g for 10 minutes and the supernatant was transferred to new tubes and precleared with 10 μl Protein A Dynabeads (Invitrogen) for 1 hour. Before immunoprecipitation, a 10% input sample was removed. For each sample, ChIP assays were performed in parallel with antibodies directed against H3K27ac (C15410196, Diagenode) and H3K4me3 (C15410003, Diagenode). Antibodies (2 μg per sample) were incubated with 10 μl Protein A Dynabeads for two hours on a rotator at 40 rpm at 6˚C. Beads were transferred to PCR-tubes and captured with a magnetic rack followed by addition of 150 μl ChIP ready chromatin. Samples were incubated at 40 rpm on a tube rotator overnight at 40 rpm at 6˚C. Next morning, the supernatant was removed and the beads were washed with 100 μl ChIP-buffer, followed by a washing step with Wash buffer A (20 mM Tris-HCl pH 8.0 (cold adjusted), 500 mM NaCl, 1 mM EDTA, 1 mM EGTA, 2% Triton X-100, 0.2% Na-deoxycholate, 0.1% SDS) and a washing step with Wash buffer B´(20 mM Tris-HCl pH 8.0 (cold adjusted), 250 mM LiCl, 1 mM EDTA, 1 mM EGTA, 2% Triton X-100, 0.2% Na-deoxycholate, 0.1% SDS) and once with TEbuffer (20 mM Tris-HCl pH 8.0). All washing steps were being performed for 5 minutes using a rotor at 40 rpm and at 6˚C, with the exception of the TE-buffer washing steps which were carried out at room temperature. After the final washing step, the beads were resuspended in 96 μl Elution buffer (20 mM Tris-HCl pH 7.5, 5 mM EDTA, 50 mM NaCl, 1% SDS) and incubated on a heating block at 37˚C with 3 μl RNaseA (20 μg/μl) (Roche). After 1 hour, 1 μl Proteinase K (20 μg/μl) was added and samples were incubated for 1 hour at 50˚C followed by 4 hours at 68˚C on thermoshaker at 1200 rpm. ChIP DNA was purified with Zymo ChIP DNA Clean & Concentrator kit in 10 μl (Zymo Research).

Flow cytometry
For flow cytometry analysis, nuclei were isolated as described above. For the ACTA1 rtTA ; TRE H2B-GFP , TA muscles were excised, minced with a razor, and homogenized in PBS containing 0.25 M sucrose and 1% BSA using an Ultra-Turrax T25. The homogenate was then incubated for 5 minutes with addition of Triton-X100 to a final concentration of 0.36% at 4˚C. Samples were filtered through a 100 μm strainer, then filtered again through a 40 μm strainer. The nuclei pellet was collected after centrifugation (3000 g for 5 minutes at 4˚C), and resuspended in sorting buffer (2% BSA/ PBS). PCM1 (HPA023370, Sigma Aldrich). Staining was performed on ice for 1 hour (1:1000). After washing, samples were stained with the secondary antibody Alexa Fluro 647 anti-rabbit IgG (A32795, Invitrogen) (1:1000) for 30 minutes in combination with Hoechst 33342 (Life Technologies) added 1:5000. Stained samples were analysed on a BD FACSCanto II or a BD LSRII. Data analysis was performed utilizing BD FACS Diva, FlowJo and Kaluza.

Next-generation sequencing
Libraries for next-generation sequencing (NGS) were created with the Swift Accel-NGS 2S Plus DNA Library Kit (Swift Biosciences) following the recommendations from the manufacturer with the exception that the ratio of beads and PEG was 1.5 and 1.3 in step Repair step I and II, respectively, and amplified with 12 PCR cycles. The libraries were sequenced on an Illumina HiSeq 2500 40 bp paired-end, at the Oslo University Hospital Genomics Core Facility (oslo.genomics.no).

Bioinformatics analysis
ChIP-Seq reads were trimmed with TrimGalore [137] and mapped to the mouse genome (mm9) using BWA [138]. Duplicated and low-quality reads were removed (q >10) using Samtools [139]. Peaks were identified with PePr for each histone modification over the respective input samples using the following parameters-sharp for peak calling and a p-value threshold of 0.05 [140]. Differently enriched regions were identified with PePr-diff with the histone modification over the respective input for each sample and using threshold on foldchange 1.5 or above and an FDR-value on 1x10 -5 for H3K4me3 and 1x10 -7 for H3K27ac or below in addition to intersect with peaks identified in their respective sample.
For gene annotation, UCSC Known Genes (mm9) with an ensemble id were used. Promoters were defined as 750 bp upstream and 250 bp downstream from transcription start site (TSS). For visualization, samples were normalized to input and BigWig files were created as described in [141] and visualized using the IGV browser. Heatmaps of ChIP-Seq signals were created with deepTools version 3.3.0 [142].
Gene ontology analyzes of DE peaks were performed with GREAT v. 4.04 [143] using single closets gene while pathway analysis was conducted with David v. 6.8 [144,145]. Redundancy reduction was performed with ReViGO [146].
Regulatory networks were visualized with Cytoscape [147]. DNase I hypersensitive enrichment data for mouse skeletal muscle were obtained from [79,80]. For the analysis, only middle point of DNase I enriched areas -/+ 100 bp inside the myo-specific H3K27ac differently enriched areas for soleus and EDL were used. Motif enrichment analysis was conducted with AME Version 5.0.4 [148] using high quality motifs from human and mouse and threshold on FDR < 1x10− 5 . For factors with similar motifs, the factor with the most significant H3K27ac enrichment at the promoter region (FDR < 1x10 -15 ) was included. Functional annotation clustering of enriched promoters was conducted with David v. 6.8 [144,145]. For external data used see S7 Table. Data have been made publicly available at the GEO database [149,150] under accession number GSE182667.  Table. GO analyses of promoter regions. Gene ontology analyses of unique promoters with H3K27ac enrichment in whole soleus and EDL muscles compared to myonuclei from same muscle type. (XLSX) S2 Table. GO analyses of unique peaks in whole muscle and myonuclei. Gene ontology analyses of unique peaks with H3K27ac enrichment in whole soleus and EDL muscles and myonuclei from same muscle type. (XLSX) S3 Table. GO analyses of differently enriched peaks in myonuclei and whole muscle. Gene ontology analyses of differently enriched H3K27ac peaks between soleus and EDL in myonuclei and whole muscle.