Application of ATAC-Seq for genome-wide analysis of the chromatin state at single myofiber resolution
Abstract
Myofibers are the main components of skeletal muscle, which is the largest tissue in the body. Myofibers are highly adaptive and can be altered under different biological and disease conditions. Therefore, transcriptional and epigenetic studies on myofibers are crucial to discover how chromatin alterations occur in the skeletal muscle under different conditions. However, due to the heterogenous nature of skeletal muscle, studying myofibers in isolation proves to be a challenging task. Single-cell sequencing has permitted the study of the epigenome of isolated myonuclei. While this provides sequencing with high dimensionality, the sequencing depth is lacking, which makes comparisons between different biological conditions difficult. Here, we report the first implementation of single myofiber ATAC-Seq, which allows for the sequencing of an individual myofiber at a depth sufficient for peak calling and for comparative analysis of chromatin accessibility under various physiological and disease conditions. Application of this technique revealed significant differences in chromatin accessibility between resting and regenerating myofibers, as well as between myofibers from a mouse model of Duchenne Muscular Dystrophy (mdx) and wild-type (WT) counterparts. This technique can lead to a wide application in the identification of chromatin regulatory elements and epigenetic mechanisms in muscle fibers during development and in muscle-wasting diseases.
Editor's evaluation
The authors have described an innovative application of ATAC-Seq for genome-wide analysis of the chromatin state at single myofiber resolution.
https://doi.org/10.7554/eLife.72792.sa0Introduction
Skeletal muscle evolved for contraction and the production of force. The main component of skeletal muscle are myofibers which are formed from the fusion of myogenic precursor cells (Buckingham et al., 2003) resulting in large postmitotic syncytia that are composed of repeating contractile units, called sarcomeres (Huxley and Hanson, 1954). Myofibers exhibit wide variations in their metabolic activity and contractile properties (Zierath and Hawley, 2004). In addition, they have a highly adaptive nature where their size, myosin heavy chain isoform, energy metabolism, and the overall skeletal muscle mass, among other characteristics, are regulated by complex processes involving rates of protein turnover (Bohé et al., 2001; Mittendorfer et al., 2005), as well as transcriptional (Quiat et al., 2011) and posttranscriptional (Weskamp et al., 2021) control of gene expression. Due to their adaptive nature, myofibers can change in response to exercise (Zierath and Hawley, 2004; Wilson et al., 2012; Dons et al., 1979), aging (Deschenes, 2004) and diseases, such as sarcopenia (Thompson, 2002; Nilwik et al., 2013) and cachexia (Roberts et al., 2013). Therefore, the study of the myofiber transcriptome and epigenome can provide key insights into how skeletal muscle adapts and changes under various stimuli, and it can potentially lead to the discovery of novel therapeutic venues for muscle related diseases.
Myofibers also act as a key signaling component of muscle stem cells (MuSCs) (Lazure et al., 2020), which are in turn required for the regeneration of muscle fibers after injury (Snow, 1977; Shadrach and Wagers, 2011; Dumont et al., 2015). Skeletal muscle is a very heterogenous tissue composed not only of myofibers and their associated MuSCs, but also numerous different non-myogenic cell types (Giordani et al., 2019). Previous studies using whole muscle sequencing captures not only the myofibers but also the other resident cell types in the muscle, making it challenging to attribute any changes in the transcriptome and epigenome specifically to myofibers as they could be due to changes in these other cell types. Recent advances in Next Generation Sequencing (NGS) now allow for high-dimensional analysis at a single-cell level. Recent studies using these technologies to study muscle tissue, such as single nucleus RNA-Seq and single nucleus ATAC-Seq have analyzed the transcriptome and epigenome of the myonuclei within the muscle fiber (Petrany et al., 2020; Dos Santos et al., 2020; Kim et al., 2020). However, they present certain limitations where they sequence all myonuclei present in the muscle and cannot distinguish between different myofibers, as well as having low sequencing depth with a limited capacity for downstream analyses.
The chromatin state plays a key role in transcriptional regulation and the determination of cellular identity (Zhu et al., 2018). Although the accessible regions make up only 3% of the total genome, it represents over 90% of known transcription factor binding sites (Thurman et al., 2012). Chromatin accessibility is a determinant of gene expression, and changes in chromatin accessibility have been identified in different biological and disease conditions such as during development (Liu et al., 2019a; Trevino et al., 2020), cancers (Corces et al., 2018; Rendeiro et al., 2016), and neurological disorders (Wang et al., 2020; Bastle and Maze, 2019). Thus, in recent years, the study of epigenetics and chromatin accessibility has become a promising field for the development of novel therapeutics. Today, ATAC-Seq is a widely used method that allows for the mapping of the accessible chromatin regions in the genome. ATAC-Seq relies on the hyperactive Tn5 transposase that fragments the accessible regions in the genome while simultaneously ligating sequencing compatible adaptors (Corces et al., 2017; Buenrostro et al., 2015). Over the years, ATAC-Seq has been applied to many different cell types and tissues (Liu et al., 2019b; Yan et al., 2020; Rocks et al., 2022). However, to our knowledge, it has not been performed on a single myofiber, possibly due to the rigidity of their membrane, high levels of mitochondria (Janssen et al., 2000; Ortenblad et al., 2018; Mishra et al., 2015), and the low number of myonuclei that are present in a single myofiber (Neal et al., 2012; Cramer et al., 2020).
Here, we have adapted OMNI-ATAC-seq to determine the genome-wide chromatin accessibility of myonuclei contained within a single Extensor Digitorum Longus (EDL) muscle fiber of a mouse. The single myofiber ATAC-Seq (smfATAC-Seq) method that we applied in this study allows for the investigation of the accessible chromatin regions of a single myofiber, without the presence of other confounding cell types. The smfATAC-Seq has a sequencing depth of approximately 6 million final reads aligned which provides approximately 30,000 peaks called. Using this method, we provide comparative analysis of chromatin accessibility between resting and regenerating myofibers, as well as their MuSC progenitors. In addition, application of this method to the study of myofibers isolated from a mouse model of Duchenne Muscular Dystrophy (mdx) (McGreevy et al., 2015) and their WT counterparts provide a genome-wide assessment of changes in chromatin accessibility in Duchenne Muscular Dystrophy (DMD). This method can be used in the future to profile the epigenetic state of myofibers in different disease conditions, and under various physiological and physical stimuli, and to identify active cis-regulatory elements in muscle fibers.
Results
Generation of ATAC-Seq libraries from a single myofiber
A single EDL myofiber of a mouse contains an average of 200–300 myonuclei (Neal et al., 2012; Cramer et al., 2020), making genome-wide analyses of the chromatin state difficult. With the advancements in next generation sequencing (NGS) and the development of the OMNI ATAC-Seq protocol (Corces et al., 2017), analysis of chromatin accessibility of samples with an input of as low as 500 cells is now possible (Corces et al., 2017). However, myofibers present additional challenges with their rigid membrane and high levels of mitochondria (Janssen et al., 2000; Ortenblad et al., 2018; Mishra et al., 2015). Here, we report a robust protocol for the successful application of ATAC-Seq on a single myofiber isolated from the EDL muscle. Our method relies on the lysis and permeabilization of a single myofiber followed by transposition with a hyperactive Tn5 transposase (Corces et al., 2017) (Figure 1). DNA fragment sizes obtained from the smfATAC-Seq were of a similar range in size as those obtained from conventional OMNI ATAC-Seq which we have performed on 5000 MuSCs that were freshly isolated by Fluorescence Activated Cell Sorting (FACS) (Figure 1—figure supplement 1). Furthermore, analysis showed that only 0.9–2.09% of reads were derived from the mitochondria in smfATAC-Seq (Table 1), suggesting that this method is highly efficient for the removal of mitochondria from mitochondria-rich myofibers. Following the removal of mitochondrial reads, there were approximately 6 million final reads aligned and 30,000 peaks called, demonstrating a sufficient sequencing depth for downstream analysis (Table 1).
smfATAC-Seq can be used to study chromatin accessibility of myofibers under different physiological conditions
In addition to adapting the OMNI ATAC-Seq method to study the chromatin accessibility of a single myofiber, we also demonstrate the application of this technique for comparative analysis of chromatin accessibility between myofibers under different conditions. For that purpose, we performed ATAC-Seq on myofibers that were in a resting (uninjured) or regenerating (injured) state. Uninjured and injured (7 days post cardiotoxin (CTX) induced injury) myofibers were isolated from wild-type C57BL/6 mice and smfATAC-Seq was performed to compare the changes in chromatin accessibility during regeneration. In addition, as a further quality control, we compared the chromatin accessibility between myonuclei within a single myofiber and 5000 freshly isolated MuSCs. This analysis not only identified accessible regions of chromatin in myofibers and MuSCs, but it also revealed a repertoire of active cis-regulatory elements in each sample.
Apart from the myofibers and their associated MuSCs, skeletal muscle also contains many non-myogenic cells such as endothelial cells, adipocytes, hematopoietic cells, fibroblasts, fibro/adipogenic progenitors (FAPs), and macrophages (Giordani et al., 2019; De Micheli et al., 2020a; De Micheli et al., 2020b). Our smfATAC-Seq method allows for the analysis of chromatin accessibility of a single myofiber without the confounding effect of these contaminating cell types. Given that the whole muscle contains non-myogenic cell types, we first compared smfATAC-Seq to an ATAC-Seq performed on whole EDL muscle by Ramachandran et al., 2019 (GSM3981673) (Ramachandran et al., 2019) for the enrichment of ATAC-Seq peaks on the genes of non-myogenic cells. We obtained the list of genes that are solely expressed in the whole muscle (RPM of at least 10) but not in the myofibers (RPM of 0) by using an RNA-Seq dataset performed on whole muscle and a single myofiber by Blackburn et al., 2019 (GSE138591) (Blackburn et al., 2019). This list represents genes that are only expressed by the muscle resident non-myogenic cell types, designated as “non-fiber muscle genes”. We determined the number of peaks in smfATAC-Seq that overlap with non-fiber muscle genes, which revealed that only 0.1% of the peaks overlapped with the top 100 non-fiber muscle genes (Table 2). In comparison, 0.33% of the peaks in the whole EDL muscle ATAC-Seq (GSM3981673) (Ramachandran et al., 2019) overlapped with the top 100 non-fiber muscle genes (Table 2). The significant difference in the overlap with the non-fiber muscle genes between the whole muscle ATAC-Seq and smfATAC-Seq suggest that the whole muscle ATAC-Seq has enrichment of peaks associated with non-myogenic genes when compared to the smfATAC-Seq, which implies that smfATAC-Seq can successfully exclude the non-myogenic cell types. In contrast, the number of overlapping peaks with all the genes expressed in whole muscle in the EDL ATAC-Seq and smfATAC-Seq were similar (Table 2). To further illustrate the absence of the non-myogenic cell types in the smfATAC-Seq samples, peaks at the promoter regions of marker genes of muscle resident cells were searched for. Specifically, Platelet and Endothelial Cell Adhesion Molecule 1 (Pecam1) was used to determine whether endothelial cells were present (Khan et al., 2005). Similarly, Resistin (Retn) and Cd45 were used as markers for adipocytes and hematopoietic cells, respectively (Steppan et al., 2001; McKinney-Freeman et al., 2009). The cell Surface Antigen Thy1 was the marker selected for fibroblasts (Agorku et al., 2019). In addition, Lymphocyte antigen 6A (Ly6a) and Adhesion G-protein-coupled receptor E1 (Adgre1) were selected for fibro/adipogenic progenitors (FAPs) and macrophages, respectively (Waddell et al., 2018; Joe et al., 2010). None of these marker genes had ATAC-Seq peaks at their promoters, indicating that only a single myofiber is processed without any other contaminating cell types (Figure 2—figure supplement 1).
smfATAC-Seq can be successfully applied to myofibers under different conditions. For instance, we applied smf-ATAC-Seq to analyze the chromatin accessibility of myofibers under resting and CTX-mediated injury conditions. In a disease condition or in injury, not all myofibers undergo damage or regenerate to the same degree. Therefore, individual myofibers within a muscle can be in different physiological and disease conditions asynchronously (Folker and Baylies, 2013). Damaged or regenerating myofibers can be visualized by their characteristic feature of centrally located nuclei (Roman and Gomes, 2018). Injured myofibers in this method were visually selected for the presence of the centrally located myonuclei by Hoechst staining and the selected myofiber was used for downstream processing with smfATAC-Seq (Figure 2A and B). The selection of a specific myofiber that our smfATAC-Seq allows for, as well as the application of trypsin to remove any associated cells that may be present, results in the sequencing of DNA fragments corresponding purely to the myonuclei within a specific myofiber.
smfATAC-Seq can identify the accessible chromatin regions of a single myofiber
To validate the quality of the ATAC-Seq data generated from a single EDL myofiber, we first investigated the profiles of the ATAC-seq samples from both injured and uninjured myofibers as well as freshly sorted MuSCs. We investigated the similarity between biological replicates for each condition (i.e. uninjured and injured myofibers and MuSCs) by visualization of ATAC-Seq peaks for the muscle-specific gene muscle creatine kinase (Ckm) (Tai et al., 2011), the housekeeping gene Gapdh and the MuSC specific gene myogenic factor 5 (Myf 5) (Rudnicki et al., 1993) (Figure 2—figure supplement 2A-C). This analysis not only indicates that smfATAC-Seq can reliably detect chromatin accessibility in a single myofiber but also the presence of comparable peaks in each biological replicate within the specific condition shows the similarity between the samples (Figure 2—figure supplement 2A-C). In addition, Pearson correlation analysis between the biological replicates showed a high correlation of ATAC-seq reads between the replicates, indicating consistency within the samples (Figure 2—figure supplement 2D-J). Furthermore, we mapped ATAC-Seq peaks with DNase-Seq from skeletal muscle (Sequence Read Archive, accession # SRX191047) and the similarity between our ATAC-Seq and previous DNase-Seq was confirmed through the common peaks present for representative genes, as visualized on the IGV (Figure 2—figure supplement 2A-C). We also analyzed the overlap between the smfATAC-Seq on EDL myofibers with the ATAC-Seq performed on the whole EDL muscle by Ramachandran et al., 2019 (GSM3981673) (Ramachandran et al., 2019). This analysis revealed that 65% of the smfATAC-Seq peaks in the uninjured myofibers overlap with the whole EDL muscle ATAC-Seq (Table 3).
Accessible chromatin regions are associated with various histone marks such as H3K27ac and H3K4me3 (Zhang et al., 2015; Berger, 2007; Barrera et al., 2008). Thus, we compared the smfATAC-Seq to publicly available datasets of ChIP-Seq on H3K27ac in EDL muscle that was previously performed by Ramachandran et al., 2019 (GSM3515022, GSM3515023) (Ramachandran et al., 2019). The comparative analysis has revealed that there were only 97 peaks in the smfATAC-Seq that did not overlap with the H3K27ac peaks, while the majority of the peaks, 6090 peaks, were common to the H3K27ac peaks present in the entire EDL muscle (Figure 2—figure supplement 2K). This demonstrates that the accessible regions that are assessed by smfATAC-Seq correspond to the regions of the chromatin marked by histones that are associated with open chromatin such as H3K27ac. Overall, these analyses suggest that smfATAC-Seq can robustly measure chromatin accessibility and identify active cis-regulatory elements in a single EDL myofiber.
Following the initial quality control and the correlation analysis, the biological replicates from the same condition were pooled for further analysis. Peak annotation analysis for MuSCs revealed that more than half of the peaks were in the intron/distal intergenic regions (i.e. enhancer regions) and about 25% of the peaks were in the promoter region (Figure 2C, Table 4). Peak annotations for the uninjured and injured single myofibers also showed a great proportion of peaks in the enhancer and promoter regions (Figure 2D and E, Table 4). In addition, an enrichment of ATAC-seq reads around Transcription Start Sites (TSS) (±1 kb) from all datasets was observed, which is a typical result that is expected from ATAC-Seq (Yan et al., 2020) (Figure 2F–H).
To further assess the quality of the ATAC-Seq data, we analyzed select genes that are expressed by either MuSCs or myofibers. For instance, in the myofiber samples, we confirmed the presence of ATAC-seq peaks in the promoter regions of Ckm, Actin alpha 1 (Acta1), Myogenic factor 6 (Myf6), and Myosin heavy chain 4 (Myh4), all of which are expressed by myofibers but not MuSCs (Lazure et al., 2020; Tai et al., 2011; Nowak et al., 1999; Stuart et al., 2016) (Figure 2I–L). On the other hand, in MuSCs we observed peaks in the promoter regions of Paired Box 7 (Pax7), and Myf5, genes that are known to be expressed in MuSCs (Rudnicki et al., 1993; Seale et al., 2000) (Figure 2M–N). Gapdh was used a housekeeping gene for all samples (Figure 2O) and Pou5f1, a marker of pluripotency, was used as a negative control (Figure 2P). These observed peaks for known expressed genes demonstrate that our method, smf-ATAC-Seq, can reliably analyse chromatin accessibility in a single myofiber.
Muscle regeneration and repair rely on the temporal expression of Myogenic Regulatory Factors (MRFs), Myf5, MyoD, Myog, and Myf6/MRF4 (Hernández-Hernández et al., 2017; Montarras et al., 1991; Asfour et al., 2018). Therefore, we assessed the chromatin accessibility of the MRFs in MuSCs and in the myofibers under homeostasis and regeneration (Figure 2—figure supplement 3). We observed peaks in the promoter regions of Myf5 only in the MuSCs but not in the myofibers and peaks in the promoters of Myog and Myf6/MRF4 were solely observed in the myofibers (Figure 2—figure supplement 3). However, we observed peaks in the promoter regions of MyoD in both the MuSCs and myofibers (Figure 2—figure supplement 3).
Uninjured and injured myofibers and MuSCs display distinct chromatin states
To show the global differences in chromatin accessibility between MuSCs, uninjured and injured myofibers, we first performed heatmap clustering of Pearson correlation coefficients on all the replicates/samples, which shows that the biological replicates within conditions are more similar to one another than to those from the other conditions (Figure 3A). This can also be observed through Principal Component Analysis (PCA) where each condition clusters separately, with the injured and uninjured myofibers being more similar to one another than to MuSCs (Figure 3B). To test whether the differences between regenerating and resting myofibers are overshadowed by their differences with MuSCs, we performed heatmap clustering of Pearson correlation coefficients and PCA analysis for injured and uninjured myofibers only, without MuSCs (Figure 3—figure supplement 1). This further highlighted how the uninjured and injured myofibers cluster separately (Figure 3—figure supplement 1).
To ensure that the differences seen between myofibers were not due to differences in fiber types, we investigated the chromatin accessibility of known marker genes for slow and fast fiber types. Troponin I2 (Tnni2) and Troponin T3 (Tnnt3), markers of fast fiber types, (Mullen and Barton, 2000; Wei and Jin, 2016) had a high level of chromatin accessibility while Troponin T1 (Tnnt1) and Myosin heavy chain 7 (Myh7), which are expressed in slow fiber types displayed no chromatin accessibility (Wei and Jin, 2016; Meredith et al., 2004) (Figure 3—figure supplement 2). This data indicates that only fast fiber types were analyzed in this study and that the differences in the chromatin state between the injured and uninjured myofibers were not due to the differences in the fiber types.
Differential analysis was performed on the ATAC-Seq peaks based on the regions defined by the consensus peak sets derived from the uninjured and injured myofibers and MuSCs conditions. The clustering analysis based on the consensus peak set shows the overall chromatin state differences between the MuSCs and injured and unjured myofibers where the replicates within a condition are more similar to one another than to those of the other conditions (Figure 3C).
The unique chromatin state in each condition can be observed through the pile-up analysis that we have performed for the peaks identified as more accessible (LFC >2) and less accessible (LFC <0.5) between uninjured myofibers and MuSCs as well as between uninjured and injured myofibers (Figure 3D). In addition, the proportion of the differential peaks corresponding to various genomic regions, such as promoters and enhancers, differs depending on whether we compare MuSCs to myofibers or compare the myofibers during regeneration and homeostasis. For instance, the differential peaks between uninjured myofibers and MuSCs were mostly found close to the promoter region ( ≤ 1 kb), whereas in the uninjured and injured myofibers comparison, a greater proportion of differential peaks were found in the intron/distal intergenic regions (i.e. enhancer regions) (Figure 3E, Table 5). This implies that MuSCs and myofibers mostly differ in their promoter accessibility, whereas myofibers during homeostasis and regeneration differ mostly at the level of distal regulatory elements.
Furthermore, we performed occupancy analysis (using DiffBind) in order to determine the unique and common peaks between the conditions. Since the occupancy analysis relies on the peak score, the distribution pattern of the peak scores for all the conditions was assessed and were found to be similar despite the observed differences in the total number of peaks between the conditions (Figure 3—figure supplement 3A). Occupancy analysis between uninjured myofibers and MuSCs revealed that MuSCs have 45,533 unique peaks while myofibers contain only 535 unique peaks which are not present in MuSCs (Figure 3F). There are also many common accessible regions as seen from 5,930 peaks that are common to both uninjured myofibers and MuSCs. This analysis suggests that myonuclei share a large number of open chromatin regions with their parental stem cells. On the other hand, the occupancy analysis between uninjured and injured myofibers revealed that there are 6,352 overlapping peaks between the regenerating and resting myofibers. However, this analysis also revealed that there are 15,352 unique peaks in the injured myofibers and only 83 peaks that are unique to the resting myofiber (Figure 3F). Furthermore, when comparing the read count between conditions around the center of unique peaks, it can be observed that each condition displays a unique open chromatin signature (Figure 3—figure supplement 3B,C).
Comparative analysis of the chromatin state between MuSCs and myofibers
To get a better understanding of the functional differences in chromatin accessibility between MuSCs and myofibers, we first performed Gene Ontology (GO Biological Process) analysis on the genes associated with the nearest peaks from the uninjured myofiber and on the genes associated with the nearest unique peaks in the myofiber compared to MuSCs. As expected, this revealed myofiber-specific biological processes such as myofiber structure and organization (Figure 4A and Figure 4—figure supplement 1A). On the other hand, GO term analysis on all the genes nearest to each peak and on the genes associated with the nearest unique peaks to MuSCs revealed biological processes such as adherens junction organization, membrane permeability, and regulation of notch signaling which play key roles in MuSCs quiescence and function (Bjornson et al., 2012; Mourikis and Tajbakhsh, 2014) (Figure 4B and Figure 4—figure supplement 1B). The analysis above also revealed that genomic regions that remain in an open chromatin state when MuSCs fully differentiate into myofibers correspond to genes that are involved in processes such as mitochondrial transport, regulation of transcription, and regulation of metabolites and energy (Figure 4—figure supplement 1C). The changes in the chromatin state between MuSCs and myofibers can also be observed from the volcano plots showing differential peaks between conditions labeled by their nearest gene (Figure 4D). For instance, genomic regions associated with genes such as muscle-specific titin-capping protein (Tcap), a component of the skeletal muscle z-disc, as well as genes that are involved in regulatory and structural functions in skeletal muscle such as titin gene (Ttn) are associated with more accessible chromatin regions in the myofiber compared to MuSCs (Markert et al., 2010; Hackman et al., 2002) (Figure 4D).
Additionally, we performed GO term analysis between uninjured and injured myofibers, which revealed that globally, the accessible regions in the resting and regenerating myofibers corresponded to genes involved in similar processes. GO term analysis on the genes associated with the nearest peaks from uninjured myofibers and from the injured myofibers as well as the genes associated with the nearest peaks which are common between uninjured and injured myofibers revealed biological processes involved in striated muscle cell development, actomyosin structure, and sarcomere organization, which are important for myofiber structural formation and for the proper function of myofibers (Figure 4A and C, Figure 4—figure supplement 1D). On the other hand, genes associated with the nearest unique peaks from the injured myofibers mostly belong to processes involved in structural components of the myofiber while the genes associated with the nearest unique peaks from uninjured myofibers correspond to genes involved in ion transport and metabolism (Figure 4—figure supplement 1E-F).
Furthermore, we analyzed the enrichment of transcription factor binding motifs in the sequences under peaks common between the injured and uninjured myofibers overlapping the promoters (±5 kb of TSS) (Figure 4—figure supplement 2A) as well as in the peaks that are unique to injured and uninjured myofibers overlapping the promoters (±5 kb of TSS) (Figure 4—figure supplement 2B). The top motifs that were enriched in the sequences under peaks common to injured and uninjured myofibers include binding site for Mef2a (Figure 4—figure supplement 2A). On the other hand, the top motifs that were enriched in the sequences under peaks unique to injured myofibers included binding sites for JUN and Stat3 (Figure 4—figure supplement 2B). However, due to the low number of unique peaks in the uninjured myofibers (Figure 3F), there was no significant motif that enriched in that peak set.
Identification of cell-type-specific pathways by global analysis of chromatin accessibility
To further understand the functional differences in chromatin accessibility between different cell types, we investigated the cell-type-specific pathways. To accomplish this, Gene Set Enrichment Analysis (GSEA) was performed on genes associated with differentially accessible peaks between the conditions. Importantly, the GSEA between uninjured and injured myofibers revealed that inflammatory response and Il2-Stat5 signaling, and injury related pathways are still operational even after 7 days of CTX-mediated injury to muscle (Laurence et al., 2007) (Figure 5A).
The GSEA between uninjured myofibers and MuSCs revealed that one of the significantly enriched pathways is myogenesis, where we can observe that genes associated with differentiation and myofiber function such as Myf6, Ckm,and Tropomyosin 2 (Tpm2) (Lazure et al., 2020; Tai et al., 2011; Buckingham, 1994; Jin et al., 2016) have higher accessibility in the myofiber compared to MuSCs (Figure 5B-D). On the other hand, genes associated with quiescence and MuSCs such as alpha seven integrin (Itga7) and Gpx3 are more accessible in MuSCs compared to the myofiber, as expected (El Haddad et al., 2012; Pasut et al., 2012) (Figure 5C and E).
Moreover, differential chromatin accessibility between the MuSCs and their myofiber derivatives show differences in pathways that are known to be important for muscle, such as Notch and TGFβ signalling (Bjornson et al., 2012; Mourikis and Tajbakhsh, 2014; Carlson et al., 2009; Girardi et al., 2021) (Figure 5—figure supplement 1). For example, increased accessibility of Notch1 is seen in MuSCs while increased accessibility of Jagged-2 (Jag2) is observed in the myofibers regardless of whether they are regenerating or homeostatic as seen by the height of the peaks at their promoters (Figure 5—figure supplement 1A-C). On the other hand, for TGFβ signalling, Noggin (Nog) shows more accessibility in the myofibers while bone morphogenetic protein-4 (Bmp4) has increased accessibility in MuSCs (Figure 5—figure supplement 1D-F). Taken together, this data shows that smfATAC-Seq is an effective method to analyze chromatin accessibility and to identify active cis-regulatory elements in a single muscle fiber as well as to compare muscle fibers under different physiological conditions.
Comparative analysis of the cromatin state between WT and MDX myofibers
To demonstrate the applicability of our method to a disease condition, we performed smfATAC-Seq on myofibers isolated from mdx mice, a model for Duchenne’s muscular dystrophy (DMD), and their WT C57BL/10ScSn counterparts. In order to solely assess the effect of the loss of dystrophin without the effect of regeneration, myofibers that were not actively regenerating were selected for processing from both conditions. As was performed in the previous cohort, we began by confirming the similarity between biological replicates by visualization of ATAC-Seq peaks for Ckm and housekeeping gene Rps2, as well by the Pearson correlation analysis between the biological replicates (Figure 6—figure supplement 1). After consistency of the samples within the conditions was established, the biological replicates from the same condition were pooled for further analysis. First, enrichment of the ATAC-Seq reads around the TSS from both mdx and WT myofibers was confirmed (Figure 6A). Peak annotation analysis revealed that most of the peaks were in the promoter and enhancer regions for both data sets (Figure 6B, Table 4). To further show that smf-ATAC-Seq can successfully assess the chromatin accessibility in a single myofiber of an mdx and WT EDL muscle, we looked at the presence of ATAC-seq peaks in the promoter regions of Ckm, Acta1, Myh4, and the housekeeping genes Gapdh and Rps2, whereas Pou5f1 was used as a negative control (Figure 6—figure supplement 2A-F). We also confirmed that the myofibers from mdx and WT conditions were fast type (Figure 6—figure supplement 2G-H) and that they exclusively represent the myonuclei, without the presence of confounding cell types in the muscle (Figure 6—figure supplement 2I-N).
To investigate the differences in the chromatin state between mdx and WT, we first performed heatmap clustering of Pearson correlation coefficients and PCA (Figure 6C and D). These analyses showed that, mdx and WT myofibers cluster separately and the biological replicates for each condition generally cluster together (Figure 6C and D). We then performed differential analysis of ATAC-Seq peaks followed by pile-up analysis for the less accessible (LFC <0.5) and more accessible (LFC >2) peaks between WT and mdx myofibers (Figure 6E). The results clearly showed that mdx and WT myofibers exhibit extensive differences in their chromatin states (Figure 6E). Peak annotation analysis on the differential peaks between mdx and WT myofibers revealed that more than half of the differential peaks were found in enhancer regions, indicating that they differ mostly at the level of distal regulatory elements (Figure 6F, Table 5). In addition, occupancy analysis revealed that there are 22,967 overlapping peaks between mdx and WT myofibers, and that mdx myofibers possess 7599 unique peaks, while WT myofibers contain 2895 unique peaks (Figure 6G).
Furthermore, to understand the functional differences in chromatin accessibility, we performed Gene Ontology (GO Biological Process) analysis on the genes associated with the nearest peak for all peaks in the mdx and WT myofibers, as well as on the genes associated with the nearest common peaks between them (Figure 6—figure supplement 3). These analyses revealed processes involved in mitochondrial transport, myofibril assembly, sarcomere organization, and striated muscle cell development, which are important for myofiber structure, organization, and function (Figure 6—figure supplement 3). However, GO term analysis on the unique peaks of each condition revealed that different processes are affected between mdx and WT. GO term analysis on the genes associated with the nearest unique peaks from mdx myofibers revealed processes that are important for myofiber structure and organization such as actin cytoskeleton organization and striated muscle cell differentiation (Figure 6H). On the other hand, GO term analysis on the genes associated with the nearest unique peaks from WT myofibers revealed processes mostly involved in metabolism (Figure 6I). Since the observed differential biological processes between WT and mdx myofibers were similar to those seen between injured vs uninjured myofibers, we then compared the overall differences in chromatin accessibility between mdx, WT, and injured myofibers. We performed heatmap clustering of Pearson correlation and PCA analysis between WT, mdx and injured myofibers (Figure 6—figure supplement 4). These analyses have revealed that injured myofibers were more similar to the mdx than they are to the WT C57BL/10ScSn myofibers. However, as expected due to the different genetic backgrounds of the mice between injured and the WT and mdx mice, WT and mdx myofibers were more similar to each other than they are to the injured C57BL/6 myofibers (Figure 6—figure supplement 4).
Finally, we determined the top motifs that are enriched in the sequences under the peaks that are common between the mdx and WT myofibers overlapping the promoters (±5 kb of TSS) (Figure 6—figure supplement 5A) as well as the sequences under peaks that are unique to mdx and WT overlapping the promoters (±5 kb of TSS) (Figure 6—figure supplement 5B-C). The top significantly enriched motifs in the peaks common between mdx and WT included Mef2a and JUN (Figure 6—figure supplement 5A) while the top motifs enriched in the peaks unique to mdx included transcription factors such as Foxo1 (Figure 6—figure supplement 5B).
Overall, this data shows that smfATAC-Seq can be reliably used to study myofibers in disease conditions, revealing that there are substantial differences in the chromatin accessibility of myofibers in the mdx mice compared to their WT counterparts.
Discussion
Analysis of the myofiber-specific chromatin state and gene expression profile is very limited due to the heterogenous nature of muscle with the presence of numerous non-myogenic cells in the tissue. Whole muscle or muscle biopsies represent a pooled result of numerous cell types which are present in the muscle tissue. To overcome this limitation, single nucleus RNA-Seq (snRNA-Seq) and single nucleus ATAC-Seq (snATAC-Seq) have been developed and performed on myonuclei to allow for the computational removal of other cell types (Petrany et al., 2020; Dos Santos et al., 2020; Kim et al., 2020). However, these methods still sequence all myonuclei present in the muscle and cannot distinguish between different myofibers within a muscle. Although snATAC-Seq provides high dimensionality, it is limited in sequencing depth due to the generation of sparse reads. Although computational pseudo bulking of snATAC-Seq can increase read numbers for comparative analysis between samples and conditions, the pooled reads represent the average of all myonuclei within the sample.
In this study, we have introduced a highly effective protocol based on the adaption of OMNI ATAC-Seq (Corces et al., 2017) to quantify chromatin accessibility of a single EDL myofiber with high resolution and sequencing depth. This method allows for comparative analysis of chromatin accessibility within and between muscle types with a potential for wide-spread use in future studies to investigate myofiber-specific epigenetic alterations in skeletal muscle.
The smfATAC-Seq protocol that we introduce in this study investigates the open chromatin state of myofibers at a single myofiber resolution, and with a high sequencing depth that allows for peak calling and differential peak analysis. Using this method, we have demonstrated that accessible chromatin regions of myonuclei contained within a single EDL myofiber can be tagmented and that high-quality sequencing ready libraries can be generated from these fragments. Sequencing of these libraries allow for sufficient depth and peak calling that can be used for genome-wide analysis of chromatin accessibility between myofibers. In this study, we have also demonstrated that smfATAC-Seq strictly investigates a single myofiber without the confounding presence of muscle resident non-myogenic cell types. Additionally, application of trypsin to the isolated myofibers effectively removes MuSCs that are associated with myofibers (Blackburn et al., 2019) and was confirmed by the absence of peaks at the promoters of known genes associated with muscle stem and niche cells. Although all the smfATAC-Seq samples sequenced in this study were fast type myofibers, this protocol can be used to distinguish between different fiber types which could be applied to study myofiber heterogeneity under various physiological conditions. smfATAC-Seq peaks associated with genes involved in muscle structure and function such as Acta1, Ckm and the myosin heavy chain cluster, indicates that our smfATAC-Seq is a robust technique to investigate genome-wide chromatin accessibility of a single myofiber.
A key implication of this technique is its applicability to the study of changes in chromatin accessibility between myofibers in different contexts. As a demonstration of this, we have performed smfATAC-Seq on uninjured myofibers as well as injured myofibers isolated 7 days post-injury to investigate the changes in chromatin accessibility that occurs during regeneration. Through Pearson correlations and PCA analysis, we showed resting and injured myofibers cluster separately, indicating the power of smfATAC-Seq to determine chromatin signature from even the minute starting material of a single myofiber. In addition, through occupancy analysis, we showed that there is a large difference in the number of unique peaks present in the injured myofibers compared to uninjured myofibers. This indicates that there are major modifications to chromatin accessibility in the context of regeneration. However, GO term analysis of genes associated with the accessible chromatin regions in both injured and uninjured fibers are very similar in the biological processes and pathways that are enriched such as striated muscle cell development, actomyosin structure and sarcomere organization, which are key factors for the proper structure and function of muscle. Despite these similarities, there are certain trends in which uninjured myofibers have increased accessibility in genes involved in energy metabolism, while injured myofibers have greater accessibility in genes involved in myogenesis and inflammatory response which is what would be expected in the case of an injury and regeneration (Tidball, 2011). Despite the increase in chromatin accessibility during injury, the accessible chromatin regions in both injured and uninjured fibers are associated with genes involved in similar biological processes. This similarity at the gene network despite differences in chromatin profile may suggest activation of multiple enhancers on core muscle structural genes in the case of injury. Another possible reason could be the length of the recovery time where at 7 days post injury, a number of genes activated early in the regeneration process may have returned to levels seen in the steady state. Previously, a study investigating the changes in the transcriptional profile of MuSCs and various muscle resident cells throughout different time points of muscle injury using single cell RNA-Seq, revealed that after seven days of regeneration most cell types returned to a state that was similar to homeostasis (De Micheli et al., 2020a). Therefore, it is possible that harvesting the injured EDL myofibers 7 days post injury allowed these myofibers to return to a state reminiscent of homeostatic myofibers. Further, our analyses of the myofibers in this study indicates that this technique can effectively compare samples between conditions and could see future use in the study of chromatin accessibility of myofibers under different biologically relevant conditions.
We have also used smfATAC-Seq to compare changes in chromatin accessibility between MuSCs and myofibers. Our data shows that the regions of open chromatin in the myofibers correspond to genes involved in structural components of the muscle, such as the z-disc, which are important for the proper functioning of the muscle. On the other hand, open regions of chromatin in the MuSCs mostly correspond to genes involved in membrane permeability, adherens junction organization and signalling pathways implicated in the regulation of MuSC function (Relaix et al., 2021). The analysis also revealed that chromatin regions that are accessible in both MuSCs and myofibers correspond to genes that are crucial for the general function of cells such as those involved in mitochondrial transport, regulation of transcription, and regulation of metabolites and energy.
Lastly, we performed smfATAC-Seq on non-regenerating myofibers isolated from the mouse model of DMD (McGreevy et al., 2015) and their WT counterparts. DMD is a type of muscular dystrophy caused by a loss-of-function mutation in the Dystrophin (DMD) gene that encodes for a protein that has a crucial role in muscle structure (Gao and McNally, 2015). Lack of functional dystrophin in DMD leads to unstable and fragile myofibers that continuously need to be regenerated, in turn leading to progressive muscle degeneration (Gao and McNally, 2015). We have not only shown that smfATAC-Seq can reliably investigate the chromatin accessibility from a single myofiber of mdx and WT EDL muscle, but through Pearson correlations, PCA and differential peak analysis we have also shown that DMD is associated with substantial alterations in chromatin accessibility in myonuclei. Through occupancy and GO term analyses, we have shown that a great proportion of peaks that are common between the mdx and WT myofibers are mostly associated with processes involved in muscle structure and organization. However, we have shown that mdx myofibers have more unique peaks compared to the WT, suggesting that chromatin accessibility of myonuclei is increased in the mdx disease model. Our analyses have shown that the unique peaks in mdx are associated with biological processes involved in myofiber structure and organization. On the other hand, unique peaks in the WT myofibers are associated with processes mostly involved in energy and metabolism. It is possible that the progressive muscle degeneration due to loss of muscle fiber integrity and stability causes mdx myofibers to compensate by increasing the activity of processes involved in muscle structure and organization while WT myofibers retain their activity in metabolism. It should be noted that the differences in chromatin accessibility that we have observed between resting and regenerating myofibers and between the mdx and WT myofibers are similar, which could be explained by the degeneration and continuous round of regeneration in the mdx myofibers.
Studies in the future could utilize smfATAC-Seq to further investigate the changes in chromatin accessibility in the mdx myofibers and investigate the changes associated with DMD at the level of chromatin to potentially investigate therapeutic avenues for this disease, as well as other muscle wasting diseases.
Overall, smfATAC-Seq is a robust molecular tool that can be used to analyze genome-wide chromatin accessibility of a single myofiber. The sequencing depth from this approach, allows for in-depth analysis, peak calling, quantitative analysis of chromatin accessibility and to identify active enhancers and promoters in a single muscle fiber. smfATAC-Seq can be used to study the epigenetic alterations that occur in muscle fibers during development, diseases, and in response to exercise.
Materials and methods
ATAC-Seq on a single myofiber
Isolation of Extensor Digitorum Longus (EDL) from cardiotoxin-induced injured muscle
View detailed protocolThe Extensor Digitorum Longus (EDL) muscle was injured by intramuscular injection of 50 µL of 5 µM cardiotoxin (CTX) (Sigma, 11061-96-4). Mice were treated with carprofen 20 minutes prior to CTX injection and were injected with CTX under anesthesia by isoflurane. Mice were sacrificed 7 days post injury and the EDL was collected from the hind limb of each mouse with the contra lateral EDL being used for the isolation of uninjured myofibers.
Dissection of EDL muscle
View detailed protocolThe EDL muscle was dissected as previously described (Blackburn et al., 2019). Briefly, the skin of the hindlimb was removed and the tibialis anterior (TA) muscle was excised with a pair of dissection scissors. The tendons of the EDL were exposed and the EDL was cut from tendon to tendon with scissors.
Isolation of a single EDL myofiber
View detailed protocolIndividual myofibers were isolated from the EDL muscle as previously described (Blackburn et al., 2019). Briefly, the intact EDL muscle was placed in a 1.5 mL eppendorf tube with 800 µL of myofiber digestion buffer containing 1000 U/mL of collagenase from Clostridium histolyticum (Sigma, C0130) in un-supplemented DMEM (Gibco, 11995–065) for 1 hr. Trypsin was added to the myofiber digestion buffer at a final concentration of 0.25% to remove the myofiber associated muscle stem cells. The EDL myofibers were then transferred into 2 mL of 1 X PBS (Wisent, 311–425 CL) in a six well-plate that had previously been coated with DMEM supplemented with 10% horse serum (HS) (Wisent, 065250). The EDL was then gently pipetted up and down with a large-bore glass pipette to disassociate the myofibers.
Selection of injured and uninjured myofibers
View detailed protocolLive myofibers in the six-well plate were stained with 2 µL of 5 mg/mL of Hoechst (Molecular Probes, H1399) in 2 mL of 1 X PBS for 5 minutes in a 37 °C with 5% CO2 incubator. The myofibers were then visualized under a microscope in DAPI channel and selected based on the myonuclei location, where myofibers with a pattern of centrally located nuclei were determined to be regenerating and picked for the injury condition. Individual myofibers were then transferred to 0.2 mL microtubes using a small-bore glass pipette coated with HS.
Lysis and permeabilization of the myofiber
View detailed protocolResidual media was removed with a pipette under a microscope. Individual myofibers in 0.2 mL microtubes were put in 10 µL of ddH2O for 5 min on ice. The ddH2O was removed with a pipette under a microscope, ensuring that the myofiber remained in the tube. The myofiber was then permeabilized with 20 µL of 0.5% Triton X-100 (Sigma, T9284) in PBS for 15 min at room temperature (RT). The permeabilization buffer was removed with a pipette under a microscope and the myofiber was washed twice with 200 µL of 1 X PBS.
Tagmentation of the myofiber by Tn5 transposase
View detailed protocolTransposition and ATAC-seq library preparation for a single myofiber was adapted from previously described OMNI ATAC-Seq protocol (Corces et al., 2017). The permeabilized myonuclei were tagmented with tagmentation mixture optimized for use on a myofiber (20 µL Tagment DNA Buffer (TD Buffer) (Illumina, 20034197), 13.3 µL PBS, 0.2% Tween-20 (Sigma, P1379-1L), 0.02% Digitonin (Promega, G9441), 1.39 µL Tn5 (Illumina, 20034197) and 4.61 µL water). Each single myofiber was incubated with 6 µL of the tagmentation mixture at 37 °C for 56 minutes with periodic shaking of the tubes every 5–7 minutes. Following the transposition with Tn5, DNA was purified using a QIAquick PCR Purification Kit (Qiagen, 28104) according to the manufacturer’s guidelines.
Library preparation
View detailed protocolThe purified DNA was PCR amplified for 15 cycles using Q5 High Fidelity DNA polymerase (New England Biolabs, M0491S) with the incorporation of Illumina Nextera XT adaptors (Illumina, FC-131–1001). The libraries were then size selected with AmpureXP Beads (Beckman, Cat# A63880) at a 1: 0.85 ratio (v/v). The size selected libraries were verified for quality control by bioanalyzer as well as verification of the library size via visualization on an agarose gel stained with GelGreen dye (Biotium, 41005). Libraries were then sequenced on NovaSeq6000 Sprime Paired End (PE) 150 bp.
ATAC-Seq on MuSCs
Isolation of MuSCs by fluorescence-activated cell sorting (FACS) for ATAC-Seq
View detailed protocolMuSCs were isolated by Fluorescence Activated Cell Sorting (FACS) as previously described (Tichy et al., 2018). Briefly, hindlimb muscles from Pax7/GFP+ mice were dissected and chopped. The minced muscles were then transferred into a 15 mL Falcon tube and digested in un-supplemented F10 media (Gibco, 11550043) with 2.4 U/mL Collagenase D (Roche, 11088882001), 12 U/mL Dispase II (Roche, 39307800), and 0.5 mM CaCl2. Digestion was performed on a shaker in an incubator at 37 °C with 5% CO2 for 30 min. Following the first digestion, digested muscles were centrifuged at 600 g for 20 s and the supernatant was transferred to a 50 mL Falcon tube with 9 mL FBS (Wisent, 080450) and was kept on ice. The remaining pellet was triturated and was digested for another 15 min with additional digestion buffer added. After the final digestion, the digested muscle mixture was transferred to the 50 ml Falcon tube containing the previously digested mixture. The digested muscle mixture was then filtered through a 40-µm cell strainer (Falcon, C352340) and was centrifuged at 600 g for 18 min at 4 °C. The pelleted cells were then resuspended in 800 μL FACS buffer that is composed of 2% FBS/ PBS (v:v), 0.5 mM EDTA (Invitrogen, AM9261) and with 0.5 μL DAPI (5 mg/mL) (Invitrogen, D3671). Resuspended cells were then filtered through 40-µm cell strainer and were transferred into polypropylene round-bottom FACS compatible tubes (Falcon, 352063). MuSCs were sorted with a FACSAria Fusion cytometer (BD Biosciences) based on negative selection for DAPI and positive selection for GFP.
Lysis and transposition of MuSCs
View detailed protocolATAC-Seq on MuSCs was performed based on the previously established OMNI-ATAC-Seq protocol (Corces et al., 2017). Briefly, five thousand MuSCs were sorted by FACS into 30 μL of the ATAC lysis buffer containing 10 mM Tris-HCl (pH 7.5), 10 mM NaCl (Bioshop, 7647-14-5), 3 mM MgCl2 (Sigma, 7786-30-3), 0.1% Tween-20 (Sigma, P1379-1L), 0.1% NP-40 (Sigma, 74385), and 0.01% Digitonin (Promega, G9441) in a 0.2 mL microtube. Cells were incubated in the lysis buffer for 5 min on ice and then 3 min at room temperature (RT). Cells were then washed with 100 µL of wash buffer composed of 10 mM Tris-HCl (pH 7.5), 10 Mm NaCl, 3 mM MgCl2 and 0.1% Tween-20, and were centrifuged at 800 g for 10 min. The pellet was resuspended in 10 µL of transposition mixture (5 µL TD buffer, 3.2 µL PBS, 0.89 µL Tn5 (Illumina, 20034197), 0.1% Tween-20, 0.01% Digitonin and 0.75 µL nuclease free water). Transposition was performed for 20 min at 37 °C while shaking the tubes every 5–7 min. The DNA was then purified using a QIAquick PCR Purification Kit according to the manufacturer’s guidelines.
Library preparation for MuSCs ATAC-Seq
View detailed protocolThe eluted tagmented DNA was PCR amplified for 12 cycles with the incorporation of Illumina Nextera XT adapters using Q5 High Fidelity DNA polymerase. The libraries were then size selected with AmpureXP Beads at a 1: 0.85 ratio (v/v). The libraries were then verified by bioanalyzer and agarose gel visualization. Finally, the samples were sequenced on NovaSeq6000 Sprime Paired End (PE) 150 bp.
ATAC-Seq data processing
Request a detailed protocolThe sequencing data was processed using the GenPipes pipeline v.3.1.5 (Bourgey et al., 2019). The raw reads were trimmed using Trimmomatic v.0.36 (Bolger et al., 2014) and aligned to the mm10 genome assembly using the Burrows-Wheeler Aligner v.0.7.12 (Li and Durbin, 2009). Reads were filtered to keep only high quality alignments (MAPQ score >20) and duplicates were removed using SAMtools v.1.3.1 (Li et al., 2009). Peak calling was performed with MACS2 v.2.1.1 (Zhang et al., 2008) using piling up of paired-end fragment mode (--format BAMPE). The peak files (bed) were filtered by removing the ENCODE black listed regions (https://www.encodeproject.org/files/ENCFF547MET) using BEDTools v2.29.1 (Quinlan, 2014). Mitochondrial reads were also removed before the analysis.
Correlation analysis between the biological replicates and clustering
Request a detailed protocolIn order to perform a quantitative comparison of the read counts within accessible regions, the overlapping peaks of all replicates were merged using BEDTools v2.29.1 (Quinlan, 2014). This set of merged peaks and the BAM alignment files were used as input for the featureCounts function of Rsubread v.2.2.6 (Liao et al., 2014) to generate a raw-count matrix. The raw counts were normalized by rlog transformation using DESeq2 (Love et al., 2014) with respect to library size. Pearson correlation coefficients were calculated based on the normalized counts for each pairwise comparison. Principal component analysis (PCA) and hierarchical clustering were also performed to evaluate the similarity between the replicates.
Peak annotation analysis
Request a detailed protocolFor each condition, the BAM alignment files of the replicates were merged and peak calling was performed with MACS2 v.2.1.1 (Zhang et al., 2008). Peak sets for each condition were annotated using the ChIPseaker v.1.24.0 (Yu et al., 2015) annotatePeak function, and the UCSC Genome Browser knownGene (mm10) table.
Obtaining coverage tracks
Request a detailed protocolThe BAM alignment files were converted to bigWig format and normalized by scaling factor (--scaleFactor) with the deepTools v.2.5.0.1 (Ramírez et al., 2014) bamCoverage function.
Enrichment of genomic signal around TSS
Request a detailed protocolThe bigWig files and the TSS coordinates obtained from the UCSC Genome Browser knownGene (mm10) table were used as input for the computeMatrix function of deepTools v.2.5.0.1 (Ramírez et al., 2014). This matrix was used for plotHeatmap function to generate the heatmap.
Identification of overlapping/unique accessible regions
Request a detailed protocolFor each comparison between the conditions, overlapping and unique accessible regions were identified with DiffBind v.2.16.2 (Stark, 2011) based on the measure of confidence in the peak call by MACS2 v.2.1.1 (Zhang et al., 2008).
Analysis of differentially accessible regions
View detailed protocolThe identification of differentially accessible regions (DARs) between the conditions was done using DiffBind v.2.16.2 (Stark, 2011) and edgeR v.3.30.1 (Robinson et al., 2010). Log fold changes were calculated, and their associated p-values were corrected for multiple hypothesis testing via the Benjamini–Hochberg procedure to obtain adjusted p-values. The DARs were annotated by their nearest gene using the annotatePeaks.pl function of Homer v.4.11 (Heinz et al., 2010).
Gene set enrichment analysis
View detailed protocolGenes nearby the DARs were ranked based on the log-fold change calculated with edgeR v.3.30.1 (Robinson et al., 2010). This ranked list of genes was used as input to perform gene set enrichment analysis with the fgseaMultilevel function of the R package fgsea v.1.14.0 (Korotkevich et al., 2021). The FGSEA-multilevel method is based on an adaptive multi-level split Monte Carlo scheme, which allows the estimation of very low p-values. The Hallmark gene sets collection from the Molecular Signatures Database (MSigDB) (Korotkevich et al., 2021) was used as a reference to identify the biological processes that were significantly enriched.
Motif enrichment analysis
Request a detailed protocolThe identification of known TF motifs found in peaks overlapping the promoter region (±5 kb of TSS) was done using the findMotifsGenome.pl function from HOMER v.4.9.1 (Heinz et al., 2010). The -size parameter was set to given to use the exact peak region as target sequence. Following the screening of HOMER’s reliable motifs library against the target sequences, the motifs enriched with a p-value less than 0.05 are returned.
Animal care
Request a detailed protocolAll procedures that were performed on animals were approved by the McGill University Animal Care Committee (UACC), protocol #7512.
Data availability
The data discussed in this study have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession numbers GSE173676 and GSE171534.
-
NCBI Gene Expression OmnibusID GSE173676. ATAC-Seq of single myofibers.
-
NCBI Gene Expression OmnibusID GSE171534. ATAC-Seq of young and aged satellite cells.
References
-
Myogenic regulatory factors: The orchestrators of myogenesis after 30 years of discoveryExperimental Biology and Medicine (Maywood, N.J.) 243:118–128.https://doi.org/10.1177/1535370217749494
-
Chromatin Regulation in Complex Brain DisordersCurrent Opinion in Behavioral Sciences 25:57–65.https://doi.org/10.1016/j.cobeha.2018.07.004
-
Notch signaling is necessary to maintain quiescence in adult muscle stem cellsStem Cells (Dayton, Ohio) 30:232–242.https://doi.org/10.1002/stem.773
-
High-resolution genome-wide expression analysis of single myofibers using SMART-SeqThe Journal of Biological Chemistry 294:20097–20108.https://doi.org/10.1074/jbc.RA119.011506
-
Trimmomatic: a flexible trimmer for Illumina sequence dataBioinformatics (Oxford, England) 30:2114–2120.https://doi.org/10.1093/bioinformatics/btu170
-
The formation of skeletal muscle: from somite to limbJournal of Anatomy 202:59–68.https://doi.org/10.1046/j.1469-7580.2003.00139.x
-
ATAC-seq: A Method for Assaying Chromatin Accessibility Genome-WideCurrent Protocols in Molecular Biology 109:21.https://doi.org/10.1002/0471142727.mb2129s109
-
The chromatin accessibility landscape of primary human cancersScience (New York, N.Y.) 362:362.https://doi.org/10.1126/science.aav1898
-
Effects of aging on muscle fibre type and sizeSports Medicine (Auckland, N.Z.) 34:809–824.https://doi.org/10.2165/00007256-200434120-00002
-
The effect of weight-lifting exercise related to muscle fiber composition and muscle cross-sectional area in humansEuropean Journal of Applied Physiology and Occupational Physiology 40:95–106.https://doi.org/10.1007/BF00421155
-
Satellite Cells and Skeletal Muscle RegenerationComprehensive Physiology 5:1027–1059.https://doi.org/10.1002/cphy.c140068
-
Glutathione peroxidase 3, a new retinoid target gene, is crucial for human skeletal muscle precursor cell survivalJournal of Cell Science 125:6147–6156.https://doi.org/10.1242/jcs.115220
-
Nuclear positioning in muscle development and diseaseFrontiers in Physiology 4:363.https://doi.org/10.3389/fphys.2013.00363
-
The Dystrophin Complex: Structure, Function, and Implications for TherapyComprehensive Physiology 5:1223–1239.https://doi.org/10.1002/cphy.c140048
-
TGFβ signaling curbs cell fusion and muscle regenerationNature Communications 12:750.https://doi.org/10.1038/s41467-020-20289-8
-
Tibial muscular dystrophy is a titinopathy caused by mutations in TTN, the gene encoding the giant skeletal-muscle protein titinAmerican Journal of Human Genetics 71:492–500.https://doi.org/10.1086/342380
-
The myogenic regulatory factors, determinants of muscle development, cell identity and regenerationSeminars in Cell & Developmental Biology 72:10–18.https://doi.org/10.1016/j.semcdb.2017.11.010
-
Skeletal muscle mass and distribution in 468 men and women aged 18-88 yrJournal of Applied Physiology (Bethesda, Md 89:81–88.https://doi.org/10.1152/jappl.2000.89.1.81
-
Comprehensive analysis of tropomyosin isoforms in skeletal muscles by top-down proteomicsJournal of Muscle Research and Cell Motility 37:41–52.https://doi.org/10.1007/s10974-016-9443-7
-
Muscle injury activates resident fibro/adipogenic progenitors that facilitate myogenesisNature Cell Biology 12:153–163.https://doi.org/10.1038/ncb2015
-
Detection of circulating endothelial cells and endothelial progenitor cells by flow cytometryCytometry. Part B, Clinical Cytometry 64:1–8.https://doi.org/10.1002/cyto.b.20040
-
The Sequence Alignment/Map format and SAMtoolsBioinformatics 25:2078–2079.https://doi.org/10.1093/bioinformatics/btp352
-
featureCounts: an efficient general purpose program for assigning sequence reads to genomic featuresBioinformatics (Oxford, England) 30:923–930.https://doi.org/10.1093/bioinformatics/btt656
-
Functional muscle analysis of the Tcap knockout mouseHuman Molecular Genetics 19:2268–2283.https://doi.org/10.1093/hmg/ddq105
-
Animal models of Duchenne muscular dystrophy: from basic mechanisms to gene therapyDisease Models & Mechanisms 8:195–213.https://doi.org/10.1242/dmm.018424
-
GREAT improves functional interpretation of cis-regulatory regionsNature Biotechnology 28:495–501.https://doi.org/10.1038/nbt.1630
-
Mutations in the slow skeletal muscle fiber myosin heavy chain gene (MYH7) cause laing early-onset distal myopathy (MPD1)American Journal of Human Genetics 75:703–708.https://doi.org/10.1086/424760
-
Developmental patterns in the expression of Myf5, MyoD, myogenin, and MRF4 during myogenesisThe New Biologist 3:592–600.
-
Distinct contextual roles for Notch signalling in skeletal muscle stem cellsBMC Developmental Biology 14:2.https://doi.org/10.1186/1471-213X-14-2
-
Isolation of muscle stem cells by fluorescence activated cell sorting cytometryMethods in Molecular Biology (Clifton, N.J.) 798:53–64.https://doi.org/10.1007/978-1-61779-343-1_3
-
BEDTools: The Swiss‐Army Tool for Genome Feature AnalysisCurrent Protocols in Bioinformatics 47:11.https://doi.org/10.1002/0471250953.bi1112s47
-
deepTools: a flexible platform for exploring deep-sequencing dataNucleic Acids Research 42:W187–W191.https://doi.org/10.1093/nar/gku365
-
Perspectives on skeletal muscle stem cellsNature Communications 12:692.https://doi.org/10.1038/s41467-020-20760-6
-
Cancer cachexia decreases specific force and accelerates fatigue in limb muscleBiochemical and Biophysical Research Communications 435:488–492.https://doi.org/10.1016/j.bbrc.2013.05.018
-
edgeR: a Bioconductor package for differential expression analysis of digital gene expression dataBioinformatics (Oxford, England) 26:139–140.https://doi.org/10.1093/bioinformatics/btp616
-
Nuclear positioning in skeletal muscleSeminars in Cell & Developmental Biology 82:51–56.https://doi.org/10.1016/j.semcdb.2017.11.005
-
Stem cells for skeletal muscle repairPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 366:2297–2306.https://doi.org/10.1098/rstb.2011.0027
-
Myosin content of individual human muscle fibers isolated by laser capture microdissectionAmerican Journal of Physiology. Cell Physiology 310:C381–C389.https://doi.org/10.1152/ajpcell.00317.2015
-
Skeletal muscle adaptations with age, inactivity, and therapeutic exerciseThe Journal of Orthopaedic and Sports Physical Therapy 32:44–57.https://doi.org/10.2519/jospt.2002.32.2.44
-
Mechanisms of muscle injury, repair, and regenerationComprehensive Physiology 1:2029–2062.https://doi.org/10.1002/cphy.c100092
-
Chromatin accessibility dynamics in a model of human forebrain developmentScience (New York, N.Y.) 367:e367.https://doi.org/10.1126/science.aay1645
-
Characterization of the chromatin accessibility in an Alzheimer’s disease (AD) mouse modelAlzheimer’s Research & Therapy 12:e29.https://doi.org/10.1186/s13195-020-00598-2
-
Post-Transcriptional Regulation in Skeletal Muscle Development, Repair, and DiseaseTrends in Molecular Medicine 27:469–481.https://doi.org/10.1016/j.molmed.2020.12.002
-
The effects of endurance, strength, and power training on muscle fiber type shiftingJournal of Strength and Conditioning Research 26:1724–1729.https://doi.org/10.1519/JSC.0b013e318234eb6f
-
ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualizationBioinformatics (Oxford, England) 31:2382–2383.https://doi.org/10.1093/bioinformatics/btv145
-
Model-based Analysis of ChIP-Seq (MACS)Genome Biology 9:e137.https://doi.org/10.1186/gb-2008-9-9-r137
-
The interplay of histone modifications - writers that readEMBO Reports 16:1467–1481.https://doi.org/10.15252/embr.201540945
Article and author information
Author details
Funding
Natural Sciences and Engineering Research Council of Canada
- Vahab D Soleimani
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Christian Young at the Lady Davis Institute for Medical Research—Jewish General Hospital- core facility for his help with fluorescence‐activated cell sorting (FACS) of muscle stem cells. We thank Dr. Michael Witcher at McGill University Department of Oncology for his critical comments and careful review of an early draft of this manuscript.
Ethics
All procedures that were performed on animals were approved by the McGill University Animal Care Committee (UACC), protocol #7512.
Copyright
© 2022, Sahinyan et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 6,222
- views
-
- 537
- downloads
-
- 14
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Genetics and Genomics
The genetic basis of severe COVID-19 has been thoroughly studied, and many genetic risk factors shared between populations have been identified. However, reduced sample sizes from non-European groups have limited the discovery of population-specific common risk loci. In this second study nested in the SCOURGE consortium, we conducted a genome-wide association study (GWAS) for COVID-19 hospitalization in admixed Americans, comprising a total of 4702 hospitalized cases recruited by SCOURGE and seven other participating studies in the COVID-19 Host Genetic Initiative. We identified four genome-wide significant associations, two of which constitute novel loci and were first discovered in Latin American populations (BAZ2B and DDIAS). A trans-ethnic meta-analysis revealed another novel cross-population risk locus in CREBBP. Finally, we assessed the performance of a cross-ancestry polygenic risk score in the SCOURGE admixed American cohort. This study constitutes the largest GWAS for COVID-19 hospitalization in admixed Latin Americans conducted to date. This allowed to reveal novel risk loci and emphasize the need of considering the diversity of populations in genomic research.
-
- Genetics and Genomics
Endocrine disrupting chemicals (EDCs) such as bisphenol S (BPS) are xenobiotic compounds that can disrupt endocrine signaling due to steric similarities to endogenous hormones. EDCs have been shown to induce disruptions in normal epigenetic programming (epimutations) and differentially expressed genes (DEGs) that predispose disease states. Most interestingly, the prevalence of epimutations following exposure to many EDCs persists over multiple generations. Many studies have described direct and prolonged effects of EDC exposure in animal models, but many questions remain about molecular mechanisms by which EDC-induced epimutations are introduced or subsequently propagated, whether there are cell type-specific susceptibilities to the same EDC, and whether this correlates with differential expression of relevant hormone receptors. We exposed cultured pluripotent (iPS), somatic (Sertoli and granulosa), and primordial germ cell-like (PGCLC) cells to BPS and found that differential incidences of BPS-induced epimutations and DEGs correlated with differential expression of relevant hormone receptors inducing epimutations near relevant hormone response elements in somatic and pluripotent, but not germ cell types. Most interestingly, we found that when iPS cells were exposed to BPS and then induced to differentiate into PGCLCs, the prevalence of epimutations and DEGs was largely retained, however, >90% of the specific epimutations and DEGs were replaced by novel epimutations and DEGs. These results suggest a unique mechanism by which an EDC-induced epimutated state may be propagated transgenerationally.