Transcriptional and Chromatin Dynamics of Muscle Regeneration after Severe Trauma

To gain insights into the healing process after administration of trauma, we performed expression profiling by RNA sequencing (RNA-seq) and small RNA-seq (miRNA-seq) (Aguilar et al., 2015) for both the injured and contralateral tissues for multiple time points (Figure 1), and lumped the datasets into three key stages (early, 3–24 hr; middle, 48–168 hr; late, 336–672 hr). Hierarchical clustering of the RNA-seq data through time revealed clusters up- and downregulated that were associated with different stages of muscle repair and regeneration (Figure 1B). For example, chemokine ligands 2 and 7 (Ccl2 and Ccl7), which are important for the recruitment of various immune cells to the injured muscle, peaked in expression in the early period (at 24 hr after injury). Annexins 1 and 2 (Anxa1 and Anxa2), cellular membrane binding proteins that are important for promoting migration of SCs, demonstrated a different expression profile and peaked in the middle period (48–72 hr). Laminin subunits β-1 and α-4 (Lamb1 and Lama4), which are cellular adhesion molecules present in the basement membrane, exhibited a peak in differential expression in the late period (336 hr). Pathway analysis of the time-clustered RNA-seq data showed an initial burst of proinflammatory and immune-response transcripts in the early period, followed by activation, proliferation, and differentiation of myogenic precursors and extracellular matrix (ECM) remodeling in the middle and late time periods, which was consistent with previous studies of muscle tissue injury and in line with a productive healing process (Tidball, 2005, Warren et al., 2007).

Regulation of muscle regeneration has also recently been shown to occur post-transcriptionally, where microRNAs (miRNAs) bind to the 3′ UTR of mRNAs and inhibit translation (Chen et al., 2005, Kim et al., 2006). Small RNA-seq was performed at multiple time points from each stage and 841 miRNAs were detected. Of these, 143 miRNAs showed dynamic behavior and reinforced many of the results observed from the RNA-seq datasets, where early upregulated miRNAs were associated with inflammation and immune system programs and middle- and late-stage upregulated miRNAs were associated with muscle repair and regeneration (Figure S1). For example, in the early period the inflammatory miR-223 was upregulated whereas its expression subsided in the middle and late periods. In the middle period, miR126b (associated with angiogenic and chemokine signaling, Zhang et al., 2013), the miR-29 family (regulator of fibrosis and collagen expression), miR-21 (regulator of phosphatidylinositol 3-kinase [PI3K]/Akt signaling), and miR-31 and miR-206 (differentiating myoblasts) were upregulated, which was consistent with pathways observed to promote myogenic differentiation during the middle period (Cacchiarelli et al., 2010, Cacchiarelli et al., 2011). To gain insights into the consequences of differential expression of these dynamic miRNAs, we assessed pairwise correlations in expression patterns with their 3′ UTR mRNA targets. We observed multiple miRNA-mRNA targets that were reciprocally regulated during the time course, and hierarchical clustering of the correlations revealed three clusters (Figure 2A ). The first cluster contained positively correlated miRNA-mRNA targets that were significantly upregulated in the early period (miR-709, miR-6538, miR-378 family with genes Atf3, Stat5b, Fas, and Btg2) and could be ascribed to functional classes related to immune regulation and DNA-damage cellular response. The second cluster was upregulated in the middle and late periods and was associated with myogenic growth and development. For example, we detected upregulation of miR-206, miR-205, miR-203, miR-1a, miR182, and miR-31 along with several target genes (Igf2, Igfbp5) that were previously identified as critical regulators of SC differentiation (Liu et al., 2012). The third cluster possessed the largest number of positively correlated miRNA-mRNA targets and was also upregulated in the middle period. This cluster was linked to transforming growth factor β (TGF-β) signaling and ECM and cytoskeletal remodeling, and contained miR-18a, miR-21, miR-335, miR34c, and miR-142 along with target genes Tgfbi, Anxa1, Sdcbp, Cav1, and Timp2.

Long noncoding RNAs (lncRNAs) are a new class of transcripts that have recently been discovered and play diverse regulatory roles in muscle differentiation and disease (Neguembor et al., 2014, Cesana et al., 2011, Gong et al., 2015). Previously identified lncRNAs were intersected with the RNA-seq datasets and 2,444 lncRNAs were detected, with 124 lncRNAs exhibiting dynamic behavior. Figure 2B shows the differentially expressed lncRNAs clustered into three groups, and similar to the mRNA and miRNA datasets, we view upregulated lncRNAs involved in inflammatory and immune-related processes (Csf1r, Il1bos, Ly9, Treml1, CD68, and Ncf1) in the early time points followed by increases in expression of lncRNAs that were previously shown to be critical for myogenic differentiation (H19 and linc-MD1) in the middle and late periods (Dey et al., 2014). Uniquely, in the early period we also detected upregulated lncRNAs that were found in other injured muscle tissues (myocardial infarction-associated transcripts Mirt1 and Mirt2), suggesting a common regulatory scheme (Zangrado et al., 2014).

SCs from injured tissues were isolated using flow cytometry (Sca-1⁻, CD45⁻, Mac-1⁻, Ter-119⁻, β1-integrin⁺, Cxcr4⁺; Figure 3A ) at time points associated with activation, proliferation, and differentiation (3, 24, 48, 72, and 168 hr), and RNA-seq and miRNA-seq was performed. A total of 17,636 mRNAs were detected and 9,953 genes were differentially expressed at one or more time points (Figure 3B). Of the time points sampled, 168 hr possessed the highest number of genes that were differentially expressed (584) compared with 434 genes for 3 hr, 108 for 24 hr, 405 for 48 hr, and 241 for 72 hr. These datasets were compared with previously published microarray datasets (Liu et al., 2013, Pallafacchina et al., 2010, Farina et al., 2012) of isolated SCs and excellent agreement was observed. Hierarchical clustering of the datasets through time revealed five clusters, and gene ontology (GO) analysis of the differentially expressed genes within each cluster illustrated significant pathway perturbations (Figure 3B). A set of genes that remained highly expressed over the time course were in the first cluster and were associated with ATP synthesis, ion channel activity, and regulation of glucose and insulin. The second and third clusters were low in expression from 3 to 48 hr and upregulated from 72 to 168 hr (Figure 3C) and contained multiple genes associated with myogenic differentiation (MyoD1, MyoG, Des, Ash2l, Hes6), cell-cycle regulation (Ccna2, Ccnb1, Dnmt1, Birc5), cytoskeletal proteins (Cdh15, Tmem8c/Myomaker), and mitochondrial metabolism. The fourth and fifth cluster displayed a relatively high level of expression from 3 to 48 hr that decreased from 72 to 168 hr and were associated with development, and contained markers of quiescent SCs (Pax7, Foxo3; Figure 3C).

miRNA-seq was performed on the fluorescence-activated cell sorting (FACS)-sorted SCs, and 341 miRNAs were detected for at least one time point (Figure 3D). Of these, 107 miRNAs were differentially expressed (Figure S1) and many of the dynamic miRNAs mirrored expression patterns observed from whole tissues. For example, in the early period miR-22 was downregulated in both SCs and at the tissue level. miR-22 has recently been shown (Lu et al., 2015) to inhibit Hdac4 expression, which is a well-known negative regulator of myogenesis. During this time period, the miR-181 family (modulators of PI3K signaling and metabolic adaption needed for subsequent proliferation) and miR-191 (promotes cell migration via induction of Tgf-β signaling) also showed decreases in expression, which is consistent with transcriptional changes in metabolism, Tgf-β-signaling, and PI3K signaling observed in the tissue at that time. In the middle time periods (≥72 hr), changes in expression of let-7c (SMAD and Tgf-β signaling inhibitor), miR-16 (regulator of proliferation and cell-cycle genes, Liu et al., 2008), miR-148a (modulates fate commitment), miR-182 (repressor of negative regulators of cell-cycle genes), and several other members of the let-7 family (associated with cell cycling and self-renewal) were detected compared with low expression in the earlier time points (3–48 hr). miR-486a and miR-486b were also upregulated at these time points, and were previously shown to potentiate proliferation by targeting p85α, Igf1r, and Igf1, which is consistent with a decrease in expression of these genes at the same time points (72 and 168 hr). Multiple other miRNAs exhibited similar expression profiles (low expression from 3 to 48 hr and high expression from >72 hr) such as 21a, 125a, 127, 199a, 206, 411, and 541, consistent with previous findings for isolated SCs (Arnold et al., 2011).

To further probe the regulation of the coding and noncoding transcriptional programs, we used chromatin immunoprecipitation followed by sequencing (ChIP-seq) to globally map the chromatin state of various cis-regulatory elements at nine time points. Validated antibodies for histone H3 lysine 4 trimethylation (H3K4me3), a modification associated with promoters, and H3K4me1 and H3K27ac, associated with poised and active enhancer regions (Ernst et al., 2011), respectively, were used to enrich chromatin. Each tissue from each time point was immunoprecipitated using the three antibodies, and ChIP-seq maps from the same antibody and time point were merged, resulting in 54 chromatin state maps covering >1.5 billion reads.

In the early time period 93,149 sites were enriched for H3K4me3 and 104,890 sites for the middle period, followed by 141,948 sites for the late period (Figure 4). Of these, 1,606 were differentially enriched in the early time points relative to the controls, 6,733 in the middle time points, and 6,121 in the late time points. To integrate sites that gained or lost the H3K4me3 modification during the time course with transcriptional activity (fragments per kilobase of transcript per million mapped reads [FPKM] > 1), we merged RNA-seq results with H3K4me3 enrichments. Approximately 58% of sites were found to associate with transcriptional activity (FPKM > 1) and 4,241 sites exhibited dynamics (acquisition or loss of the histone modification) through at least one stage (early to middle, middle to late, early to late).

In contrast to H3K4me3, which has previously been shown to be largely static during chromatin remodeling events (Garber et al., 2012), H3K4me1 and H3K27ac, which demarcate enhancer elements, were observed to be highly dynamic (Rada-Iglesias et al., 2011). H3K4me1 resides on both poised and active enhancers, whereas H3K27ac marks active sites (Creyghton et al., 2010). In the early time period 93,938 sites enriched for H3K4me1 were identified, followed by 106,353 sites for the middle period and 28,663 sites for the late period. H3K27ac showed similar behavior with 32,285 sites enriched in the early time period, 43,780 sites in the middle period, and 13,748 sites in the late period. The highest number of sites that acquired H3K27ac was found in the 72-hr period (28,178 sites) compared with 11,535 for the 48-hr stage and 4,067 for the 168-hr stage (Figure 4). In contrast to H3K4me3 enrichments, which were primarily found at or near transcriptional start sites, the enhancer marks H3K4me1 and H3K27ac were more broadly distributed across inter- and intragenic loci (Figure 4). Since H3K4me1 demarcates both poised and active enhancers and H3K27ac marks only active enhancers (Ostuni et al., 2013), the total number of enhancer elements (identified by H3K4me1⁺, H3K27ac⁺, and H3K4me3⁻) was determined, and the ratio of the enriched enhancers was quantified as active (H3K4me1⁺ and H3K27ac⁺) or poised (H3K4me1⁺ and H3K27ac⁻). A wide spectrum of binding patterns was observed for the two enhancer categories across the nine time points, with the majority of differential active enhancers occurring at intragenic and intergenic loci for all time periods. The highest number of differential active enhancers (injured⁺ and control⁻) was found for the middle period, with the largest number of active enhancers occurring in the 72-hr period.

GO term analysis of the enriched differential enhancer peaks for the early period demonstrated over-represented terms, similar to the transcriptional groups (immune response, chemokine and cytokine signaling, inflammation, and metabolism of lipids and lipoproteins; Figure 5A ). The result of these chromatin and transcriptional enrichments in the early period is consistent with upregulation of cytokines and chemokines (Il-6 and TNFα), and increases in expression of pleiotropic transcription factors such as components of the activator protein 1 (AP-1) complex (c-Fos, Atf3, JunB), STATs, NF-κβ, EGR, and PU.1 (Figures 5B and S2), which stimulate a permissive chromatin state and have previously been shown to induce SC activation genes such as Myod1, Myf6, and c-Myc (Toth et al., 2011).

In the middle time period (48–168 hr after injury), multiple genes associated with Tgf-β signaling were differentially expressed, such as several Smads (Figure S3), which have been shown to interact with chromatin remodeling complexes such as histone acetyltransferases p300 and CBP (CREB-binding protein) to induce H3K27 acetylation (Mullen et al., 2011). Consistent with this view, we observed the highest number of sites that acquired H3K27ac in the middle time period (Figure 4). GO annotation of the differential enhancer peaks for the middle period revealed multiple enrichments for anti-inflammatory cytokines, G-protein-coupled receptors, and the Rho family of small GTPases (Figures 5 and S4). This result is in agreement with the peak in expression of SRF target genes (Kuwahara et al., 2005) and anti-inflammatory cytokines interleukin-4, -10, and -13 (Il-4, Il-10, Il-13), which induce M2 macrophage polarization and are essential components for resolution of inflammation and tissue repair (Figure 5). Increases in expression of these factors were also mirrored by upregulation of the macrophage-derived matrix metalloproteinase 12 (Mmp12), which cleaves and inactivates CXC chemokines (Cxcl1, -2, -3, -5, and -8) and monocyte chemotactic proteins (Ccl2, -7, -8, and -13), inhibiting leukocyte flux to the injured site and abrogating the amount of proinflammatory molecules present in the injured tissue (Figure S5).

A primary determinant of successful muscle regeneration after injury is ECM remodeling (Calve et al., 2010), which provides a scaffold for proliferating SCs to migrate, differentiate, and fuse in the correct orientation. In the middle period, increases in expression and multiple enhancer enrichments were viewed for genes associated with ECM remodeling, such as components of the basement membrane (myoferlin, laminin, collagen VI genes, annexins Figures 6 and S5), glycoproteins, and Mmps. In addition, during the middle period chromatin remodeling and upregulation of genes associated with muscle development and architecture were also detected (Figure 6).

To understand this critical step further, we conducted additional analysis of the differential active enhancer sites during the middle period by evaluating the number of enriched active enhancers that overlapped with MyoD and MyoG ChIP-seq data (Cao et al., 2010, The Mouse ENCODE Consortium, 2014, Mousavi et al., 2013). The largest increase in MyoD binding sites at active enhancers was found at 72 hr (Figure S6), which is consistent with the largest increase in expression (Figure 6E) and known association between MyoD and p300/CBP on E-box motifs of target genes during myogenic differentiation (Palacios and Puri, 2006). The greatest increase in MyoG binding sites at active enhancers was viewed at 336 hr, which agrees well with observation of late myogenic differentiation programs in that time period. Compiling the enriched enhancer sites coincident for MyoD or MyoG binding and performing GREAT (genomic regions enhancement of annotations tool) analysis revealed over-represented terms nominally associated with myogenic differentiation such as different types of growth factor signaling (insulin growth factor [Igf], fibroblast growth factor [Fgf], hedgehog), p38 mitogen-activated protein kinase (p38-MAPK), and PI3K signaling (Figure S6 and discussed further below).

IGF signaling plays a crucial role in muscle repair and regeneration after injury, most notably during the transition between myoblast proliferation and differentiation whereby Igf2 activates the Akt pathway and targets bound by MyoD. In line with this observation, an increase in expression of Igf1 was detected in the early period (Pelosi et al., 2007), which subsided in the middle period when several Igfbps and Igf2 and members of the Akt pathway were observed to be upregulated along with miR-483 (which is embedded in the second intron of Igf2 and is an SRF target) (Figure 7A ; Schiaffino and Mammucari, 2011). One component of the IGF pathway, the PI3K component, is essential for skeletal muscle regeneration and is composed of three primary classes (classes I, II, III) that are structurally and functionally distinct. The three class IA PI3K p110 catalytic subunits (α, β, and δ) are expressed in skeletal muscle (Matheny and Adamo, 2010) and of these subunits, p110α and p110β have been shown to positively influence myoblast proliferation and differentiation (Matheny et al., 2012, Matheny et al., 2015). In the early and middle periods following injury, we observed differential enhancer binding and increases in expression of PI3K catalytic subunits p110α (Pik3ca), p110β (Pik3cb), and p110δ (Pik3cd) as well as the PI3K p85α (Pik3r1) regulatory subunit for the injured muscle compared with the contralateral control (Figure 6A). Increases in expression and active enhancers were also viewed for class IB p110γ (Pik3cg) and its regulatory subunit p101 (Pik3r5) during the healing process, which are expressed abundantly in immune cells. TF binding analysis of the enriched enhancer regions for the regulatory subunit p85α showed motifs for signal mediators of Il-4 and the anti-inflammatory response, which is consistent with activation of pathways during this time period (Figure 7B).

The temporal increases in PI3K expression (and other positive myogenic regulatory pathways) between 24 and 72 hr post injury were also concomitant with increases in expression of negative feedback regulators such as Mg53 (disruptor of upstream IGF signaling via ubiquitin ligases, Yi et al., 2013), Xbp1 (Acosta-Alvear et al., 2007), and Usf1 (which compete with or co-occupy binding sites for MyoD and Mef2), and Id1, Id2, and Snai1 (which reduce the binding affinity of MyoD and Mef2; Soleimani et al., 2012). The negative regulators began to decrease in expression starting at 72 hr (Figure S7), further reinforcing the observation of myogenic differentiation at this time.