We sought to determine the mechanisms underlying failure of muscle regeneration that is observed in dystrophic muscle through hypothesis generation using muscle profiling data (human dystrophy and murine regeneration). We found that transforming growth factor β–centered networks strongly associated with pathological fibrosis and failed regeneration were also induced during normal regeneration but at distinct time points. We hypothesized that asynchronously regenerating microenvironments are an underlying driver of fibrosis and failed regeneration. We validated this hypothesis using an experimental model of focal asynchronous bouts of muscle regeneration in wild-type (WT) mice. A chronic inflammatory state and reduced mitochondrial oxidative capacity are observed in bouts separated by 4 d, whereas a chronic profibrotic state was seen in bouts separated by 10 d. Treatment of asynchronously remodeling WT muscle with either prednisone or VBP15 mitigated the molecular phenotype. Our asynchronous regeneration model for pathological fibrosis and muscle wasting in the muscular dystrophies is likely generalizable to tissue failure in chronic inflammatory states in other regenerative tissues.
Duchenne muscular dystrophy (DMD) is a relatively common inborn error, caused by lack of dystrophin in the membrane cytoskeleton of muscle fibers (Hoffman et al., 1987). Myofibers show increased plasma membrane instability and can undergo cell death (necrosis) either individually or in groups. Necrotic myofibers are able to regenerate within the preexisting basal lamina through satellite cell activation, myoblast proliferation, fusion, and maturation—a process that takes ∼2 wk. Robust regeneration is seen in young DMD patient muscle as well as dystrophin-deficient dogs, mice, cats, and pigs (Chen et al., 2005; Kornegay et al., 2012a; Klymiuk et al., 2013). Successful regeneration is associated with retained muscle function at a young age. However, as a DMD patient ages, muscle regeneration gradually fails, with endomysial fibrosis and fatty replacement of the muscle. The fibrofatty infiltration is typically commensurate with muscle wasting and loss of muscle function. Some muscles in both dystrophin-deficient DMD patients and dogs, such as the extraocular muscles and sartorius, seem able to maintain regeneration and are functionally spared (Kornegay et al., 2012b; Nghiem et al., 2013). Most dystrophin-deficient muscles in cats and mice are able to successfully regenerate throughout the normal life span (Gaschen et al., 1992). The failed regeneration seen in the DMD muscle has features that are shared with aging muscle, including chronic inflammatory state and loss of mitochondrial function (Baron et al., 2011). However, the age-, muscle-, and species-specific progressions from successful to failed regeneration are poorly understood, and there has not been a unifying model able to explain the observed transition to muscle wasting and the loss of muscle function (Serrano et al., 2011).
The gradual loss of muscle tissue and disorganization of muscle structure in most muscles in a DMD patient has been attributed, in part, to a progressive exhaustion of satellite cell function in response to the chronic drive to repair degenerating muscle fibers (Blau et al., 1983). However, there are conflicting studies regarding the number and nature of muscle stem cells, and loss of satellite cells with age is questioned as there appear to be higher numbers of satellite cells in DMD muscle than there are in normal muscle (Watkins and Cullen, 1986, 1988; Maier and Bornemann, 1999; Kottlors and Kirschner, 2010). Reports of diminished DMD myoblast proliferation have been associated with low minimal telomere length in DMD whole muscle (Aguennouz et al., 2011) and myoblasts (Decary et al., 2000). However, this remains controversial as another investigation detected no overall reduction in mean telomere length in DMD muscle cells (Oexle et al., 1997).
An alternative hypothesis promoted for failed regeneration is the excessive connective tissue proliferation (fibrosis), a feature in common with many chronic inflammatory states (Serrano et al., 2011). Fibrosis could impair vascularization of muscle tissue, create an insufficient environment for effective satellite cell proliferation and differentiation, and alter cell–cell communication during myofiber regeneration. Histologically, fibrosis is characterized by the exaggerated expression and deposition of extracellular matrix components, particularly collagen I. TGF-β, a key cytokine promoting fibrosis, appears to have a critical role in fibrosis in DMD (Chen et al., 2005; MacDonald and Cohn, 2012) as well as most fibrotic states in most tissues (Akhurst and Hata, 2012; Bowen et al., 2013). Strong induction of TGF-β pathways is seen in both successful and unsuccessful remodeling (Karkampouna et al., 2012). Indeed, TGF-β is important for normal muscle regeneration, and the molecular mechanisms underlying a change in normal versus pathological roles of TGF-β pathways have remained elusive for muscular dystrophy and other chronic inflammatory states (Bhattacharyya et al., 2013). The likely importance of TGF-β pathways in DMD is underscored by the recent finding of two genetic modifiers of disease severity, both of which map to TGF-β pathways (SPP1/osteopontin [Pegoraro et al., 2011; Bello et al., 2012]; LTPB4 [Flanigan et al., 2013]).
Here, we used a data-driven approach involving three large mRNA profiling experiments (two human muscle biopsy and one murine regeneration series) to generate a novel hypothesis regarding the cellular underpinnings of failed regeneration in muscle disease (asynchronous regeneration). The hypothetical model was then tested using an experimental approach of asynchronous bouts of staged degeneration/regeneration in localized microenvironments of wild-type (WT) muscle.
Development of bioinformatics methods for increased sensitivity of biochemical network selection
It is well documented that different analytical methods applied to the same microarray dataset can provide very different interpretations of the data, often with only 30% concordance of lists of differentially regulated genes (Bakay et al., 2002; Seo et al., 2004, 2006; Irizarry et al., 2005). We reasoned that increased sensitivity for selection of differentially regulated genes could be obtained by using multiple probe set algorithms, multiple microarrays (e.g., U133A and U133B arrays on the same samples), and multiple queries of existing literature for biochemical pathways (Fig. S1). Such an iterative approach would be expected to greatly increase sensitivity of gene selection but also increase false positives. However, using an integrative approach with a composite score of gene selection and associated networks, we could decrease false positives by focusing on the top-ranked composite networks (Fig. S1). In effect, as different sources of data are assembled to derive a probability determining the specificity of a network, the most highly scored networks have considerable statistical power and significance.
To our knowledge, this two-tier composite approach has not been previously described. This approach would be expected to be most effective on larger and more complex datasets, and for this reason, we used a human 12-disease group, 117-biopsy dataset, and 234-microarray dataset as a dataset for model generation, with multiple validation and model extension datasets.
Visualization of domain knowledge identifies a potential disease grouping by the histological phenotype of fibrosis and failed regeneration
Seven heat maps with clear patterns of association of disease groups were generated and visually inspected for possible shared pathogenesis between the clustered disease groups. Many of the disease groupings were expected and were not studied further. For example, juvenile dermatomyositis (JDM) is classified as an autoimmune disorder and shows extensive inflammation (immune cells) infiltrating the muscle (Tezak et al., 2002; Chen et al., 2008). Indeed, a strong clustering of JDM distinct from all other disorders was seen with both A and B chips and all probe set algorithms, and the members of this cluster were transcripts involved in immune cell infiltration (unpublished data). Likewise, most muscular dystrophies show some degeneration and regeneration of myofibers, and these gene clusters were evident as a recapitulation of embryonic developmental programs (unpublished data).
One particular grouping of the seven primary heat maps was believed to be particularly interesting based upon domain knowledge and was studied in some depth (Fig. 1 A). This grouping showed a strong up-regulation of transcripts in three disorders: DMD, JDM, and ED-D (Emery–Dreifuss muscular dystrophy autosomal dominant; LMNA [lamin A/C] mutations). The same cluster was up-regulated, but to a lesser extent, in Becker muscular dystrophy (BMD; dystrophin mutation), LGMD2A (limb girdle muscular dystrophy type 2A; CALP3 [calpain 3] mutation), and LGMD2B (limb girdle muscular dystrophy type 2B; DYSF [dysferlin] mutation). The clustering was seen for the A microarray and PLIER (Probe Logarithmic Intensity Error) probe set algorithm (Fig. 1, heat map) and also all other applied probe set algorithms (not depicted). This grouping could be potentially significant, as these three disorders often show particularly severe histological patterns of fibrosis and failed regeneration early in the disease process.
As shown in Fig. S1, this grouping was expanded to both A and B microarrays and both PLIER and dCHIP (DNA chip analyzer) probe set algorithms. Applying IPA (Ingenuity Pathway Analysis) to this clustering, we found the top-ranked networks to share many component members, and these were merged into a 56-member network (Fig. 1 B). Fold changes were displayed using both node size and intensity of color, with the numeric value shown under each node. Statistical support for differential regulation of component members is shown for the PLIER algorithm (Table 1).
This network was centered on TGF-β, a well-studied cytokine associated with pathological fibrosis in many tissues (Akhurst and Hata, 2012; Bowen et al., 2013) and included many of the known components of fibrotic tissue (collagens and lysine oxygenases). Another notable component of the network was HIF1A (hypoxia-inducible factor 1, α subunit), a known transcription factor strongly induced by tissue ischemia (Koh and Powis, 2012). The nature of the network suggested that it was simply reflective of a more severe pathology (greater fibrosis, greater failed remodeling, and regeneration), and DMD, JDM, and LMNA often (but not always) show more severe pathology at the initial diagnosis compared with the other disorders in the dataset.
To test whether the network was more reflective of tissue pathology than clinical and molecular diagnosis, a second validation dataset was generated that was then analyzed by unsupervised clustering. This new dataset contained 49 patient muscle biopsies analyzed on HG-U133 plus 2.0 microarrays. Samples studied were 6 normal volunteers, 17 DMD, 11 BMD, 7 LGMD2I (limb girdle muscular dystrophy type 2I; FKRP [fukutin-related protein deficiency]), and 8 LGMD2B.
Unsupervised clustering using the 56 network members from Fig. 1 resulted in clustering of the 49 biopsies into two distinct groups indicated by the dendrogram at the top (Fig. 1, purple left group with uniform down-regulation of transcripts and blue right group with uniform up-regulation; and Fig. 2). The two branches did not correspond well to clinical/molecular diagnosis, with Becker, Duchenne, FKRP-deficient and dysferlin-deficient patients seen in both groups. To test whether the segregation of groups was caused by histopathology of the muscle studied, the archived interpretation of the severity of the dystrophy was accessed (all read by the same author, E.P. Hoffman; Fig. 2, normal, mild, moderate, or severe dystrophy groups). Mild was defined by evidence of degeneration and regeneration but little or no fibrosis or inflammation. Moderate showed increased inflammation and fiber size variation but minimal or moderate fibrosis. Severe showed fibrofatty replacement of muscle and evidence of failed regeneration. Normal or mild pathology was seen in 20/22 patients in the left branch but only 3/27 in the right branch. In contrast, moderate or severe pathology was seen in 2/22 patients in the left branch and 24/27 in the right branch. Thus, we concluded that this cluster and associated networks were more reflective of severity of histological pathology than clinical/molecular diagnosis.
Query of the putative failed regeneration network in normal staged tissue regeneration identifies time-resolved subnetworks
In many respects, the identification of the TGF-β/fibrosis network was quite expected, as fibrosis is associated with end-stage tissue pathology and with failed regeneration in many tissues, including muscle. However, the majority of patient muscle biopsies were studied at initial diagnosis, before end-stage clinical disease. Moreover, there is increasing evidence that TGF-β is important for normal staged muscle regeneration (Philippou et al., 2008; Kollias and McDermott, 2008; MacDonald and Cohn, 2012). To query the relationship of this cluster with normal regeneration, we overlaid these network members over our previously published 27–time point murine in vivo muscle degeneration/regeneration time series (Fig. 3; Zhao et al., 2003, 2006; Zhao and Hoffman, 2004).
All network members for which we could identify cross-species transcript units and probe sets did indeed show significant regulation of these transcripts in normal regeneration both by heat map (Fig. 3 A) and by statistical comparisons (Table 2). However, we noted from the heat map (Fig. 3 B) that different members of this pathology-associated network were expressed at different time points during normal regeneration. To visualize this temporal-specific pattern, we grouped the 27–time point murine series into four time frames corresponding to relatively distinct periods of muscle regeneration (myoblast proliferation—days 1–3; myotube formation—days 3.5–4.5; myofiber maturation—days 4.5–14; late-stage adaptation—days 20–40; Fig. 3 B). This parsed the TGF-β network in Fig. 1 B into a series of time-resolved subnetworks shown in Fig. 3 B, where the relative expression levels were again coded by both symbol size and color intensity within each subnetwork. The extracellular matrix components commonly associated with pathological fibrosis were seen in the myofiber maturation stage but not other stages.
This analysis evoked a novel hypothesis regarding the distinction between successful regeneration and unsuccessful remodeling (fibrosis). We reasoned that normal regeneration is synchronous, with specific pathways and pathway members having appropriate temporal expression patterns during the 2–3-wk muscle regeneration process. In contrast, dystrophic muscle was an asynchronous remodeling process, with different regions of the small muscle biopsy studied in different stages of regeneration.
This analysis established the hypothetical model in which asynchronous remodeling may be the driver for failed regeneration, muscle wasting, and early death in the severe dystrophies. By this model, the complete TGF-β–centered network shown in Fig. 1 is only possible in dystrophic muscle. In normal muscle, it is parsed into subnetworks resolved by time. We then sought to test this model experimentally.
Establishing a model of asynchronous muscle regeneration
We reasoned that a potential experimental model to test asynchronous regeneration would be to induce localized regions of muscle degeneration/regeneration in WT mice, with repetitive regional injuries in close proximity within the muscle at defined time points. The experimental design is shown in Fig. S2. In brief, low doses of notexin (10 µl notexin at 5 mg/µl) were injected intramuscularly along with a blue tattoo dye. At a defined time point later, a second injection was performed with the same low dose of notexin, using a red tattoo dye. Mice were sacrificed 13 d after the second injection, a time point when regeneration should be largely complete for both areas. We then focused on specific areas of the muscle, where both injection sites were separated by a relatively small area (50–300 myofibers)—the hypothesis being that myofibers in the microenvironments immediately at the injection sites would be synchronously regenerating, whereas the myofibers in between the two injection sites would receive regeneration cues from both sites (cross talk microenvironment).
In an initial experiment, injections were spaced by 0 (two injections at the same time), 1, 2, 4, 5, and 10 d (four muscles at each time point). Cryosections from the belly of each muscle were analyzed for fibrotic content (Fig. S3 A). This analysis showed that increasing time between injections resulted in increasing content of fibrosis of the muscle (Fig. S3 A). This suggested that the timing of injections influenced the amount of fibrotic tissue deposited in the muscle as a whole tissue.
We hypothesized that asynchronous remodeling may result from interference in intercellular signaling when two nearby groups of myofibers are at different stages of regeneration. To test this, we used laser capture microscopy (LCM) to collect groups of myofibers at the sites of injury, as identified by red and blue tattoo dyes, and the myofibers in between the sites of injury. We focused on the mouse groups of 0, 4, and 10 d between injections.
Five regions were excised using LCM from each muscle: (1) blue dye area (injection 1; 17 d after injection for 4-d series and 23 d after injection for 10-d series), (2) red dye area (injection 2; 13 d after injection for both series), (3) muscle fibers in between the two injections sites (called in-between), (4) normal nonregenerating (intact) myofibers (identified by peripheral nuclei), and (5) single regeneration site (single dye area with no evidence of the second dye). These areas were identified by centralized nuclei myofibers located in the external area of one injury site (either blue or red).
LCM-excised regions were used for mRNA isolation and mRNA profiling using Illumina bead arrays. Unsupervised chip clustering showed excellent signal/noise balance, as indicated by the clustering of microarray replicates of LCM from adjacent muscle sections (Fig. S2 C). For data analysis, we first compared region 4 (uninjured) versus region 5 (single regeneration). There were very few differentially regulated transcripts, suggesting that the muscle in region 5 had indeed successfully regenerated to completion. Likewise, there were few statistically significant differences between region 1 (first injury) versus region 5 and region 2 (second injury) versus region 5 (single site), consistent with successful remodeling (unpublished data).
In contrast, comparison of the tissue in-between injuries (area 3) showed mRNA expression profiles that were quite distinct from neighboring regions and also remarkably different for the 4-d series (Fig. 4) versus the 10-d multiinjury series (Fig. 5). Specifically, the 4-d in-between myofiber region showed a relative loss of mitochondrial-related transcripts and an increase in inflammatory transcripts (Fig. 4, left). The 10-d in-between myofiber region showed increase of extracellular matrix components (Fig. 5, left). Statistical comparison of fold changes and p-values between the 4- and 10-d multiple injection sites showed highly significant differences, with the mitochondrial changes specific to the 4-d time series (Table 3) and the extracellular matrix transcripts specific to the 10-d time series (Table 4). Alignment of the subcellular localizations of the components of the top-ranked Ingenuity networks for the 4 and 10 d confirmed the preponderance of intracellular protein changes in the 4 d, in contrast to the extracellular matrix components in the 10 d (Fig. S3, B and C).
We then overlaid the gene networks in Fig. 4 over normal staged muscle regeneration, separating the up-regulated and down-regulated transcripts (Fig. 4, right). This showed that the down-regulated mitochondrial transcripts were highly expressed in uninjured and fully regenerated muscle but were strongly down-regulated both in the 4-d in-between area as well as the 4-d time point during normal staged regeneration. This was consistent with the myofibers in between the two injury sites having suspended regeneration at the 4-d time point. Similarly, the up-regulated inflammatory transcripts in the in-between area separated by 4 d were highly expressed at the 4-d time point of normal staged regeneration but not at the later fully regenerated time point (Fig. 4). Thus, both the down-regulated mitochondrial transcripts and the up-regulated inflammatory transcripts in the in-between area appeared arrested at the time point separating the asynchronous injury times.
The 10-d staged injury series was also consistent with the suspended development model. Connective tissue remodeling transcripts were highly up-regulated in the in-between area separated by 10 d, and this pattern was consistent with the 10-d time point in normal regeneration (Fig. 5, right). Statistical analysis showed that this pattern was largely specific to the 10-d staged remodeling and not seen in the 4-d staged remodeling (Table 4).
The microarray data from the 10-d interval injury model showed that the myofiber region in between asynchronous injuries showed inappropriate persistent expression of connective tissue-remodeling transcripts, and this would be consistent with the pathological findings of localized tissue fibrosis. Although we had previously shown an overall increase in fibrosis that was highest in the day 10 staged series muscle (Fig. 3 A), we refined this analysis to the localized area in between the two injection sites (Fig. 6). Histological visualization of the area between the two injection sites showed dramatic fibrotic replacement of muscle in the 10-d series, which was not seen with the 0- and 4-d series (Fig. 6).
Targeted experiments confirm aberrant expression of fibrotic transcripts (10 d apart) and persistent inflammation and loss of mitochondrial capacity (4 d apart)
To validate the mRNA microarray findings, we selected specific transcripts and proteins to evaluate using an independent series of mice that had been subjected to the same 4- and 10-d asynchronous muscle injuries. Tattoo dyes were used to visualize the first and second injection sites. Tissue for quantitative RT-PCR (qRT-PCR) was obtained by LCM. Immunohistochemistry and enzyme histochemistry were also performed.
The 4-d interval injury model showed an increase in inflammatory transcripts and a relative loss of mitochondrial-associated transcripts (Fig. 4). qRT-PCR confirmed the mRNA profiling data, showing preferential loss of mitochondrial transcripts in the region in between the two injections sites, relative to other regions of the same muscle tissue and the area between two notexin injections made on the same day (Table 5). To confirm the relative loss of mitochondrial oxidative capacity predicted by expression microarray and qRT-PCR microarray data for the 4-d in-between areas, we used standard enzyme histochemistry (NADH, cytochrome c oxidase [COX], succinate dehydrogenase [SDH], and a combination of NADH and COX enzymatic staining). The region between the 4-d interval injection sites showed a relative loss of mitochondrial enzymatic activity stains compared with other injection paradigms, confirming the microarray data (Fig. S4).
To validate the up-regulation of proinflammatory transcripts, we studied three proteins by immunofluorescence localizations in three double injury mice (injuries spaced 0, 4, and 10 d apart; CD11b, CD74, and CD163; Fig. 7). These proteins were selected based on their role in proinflammatory pathways as well as DMD pathogenesis. CD11b is a leukocyte marker that induces innate immune response resulting in phagocytosis and cell-mediated cytotoxicity (Solovjov et al., 2005). CD11b immunostains showed the highest expression in the in-between region at 4 d but not 10 d or the simultaneous injection control (Fig. 7 A).
Human leukocyte antigen class II histocompatibility antigen γ chain (CD74) was selected because of its role in regulating adaptive immune response upon injury and its induction of NFK-β pathways (one of the earliest proinflammatory markers of DMD; Cresswell, 1994; Chen et al., 2005). CD74 showed mRNA up-regulation in the 4-d between area and not the 10-d area (Fig. 4 and Table 3). Immunostaining of CD74 showed expression in the region between the 4-d injection sites but not the region in between the 10-d injection sites (Fig. 7 B). CD74 was also not seen in the injection sites themselves.
CD163 is a marker of specific macrophage subtypes (M2c) in dystrophin-deficient muscles, is associated with increased tissue repair (Villalta et al., 2009), and is a known acute phase-regulated receptor (Madsen et al., 2001). The immunostaining patterns for CD163, cell surface markers of monocytes and macrophages class II (M2c), showed similar patterns (Fig. 7 C).
For the 10-d interval injuries, we had shown increases in connective tissue by histological methods (Fig. 6 and Fig. S3 A). To further validate this finding by immunostaining of specific proteins, CD206 (MRC1 [mannose receptor, C type 1]; macrophage marker) and MMP9 (metalloproteinase implicated in dystrophy pathogenesis; Dahiya et al., 2011a,b; Miyazaki et al., 2011) were studied (Fig. 8). Both showed strong expression of the myofiber in-between region of injuries spaced 10 d apart, with lack of immunostaining in 4-d injuries or 0-d control.
LCM of human DMD patient muscle shows patterns consistent with asynchronous remodeling
Muscle biopsies from DMD patients have long been recognized to show grouped degeneration and regeneration and focal fibrosis. This histological pattern is consistent with the asynchronous remodeling shown in our murine experimental model. To provide evidence for asynchronous remodeling in human DMD muscle at the molecular level, we studied a frozen diagnostic muscle biopsy from a young child with DMD (3 yr). Frozen sections were stained by hematoxylin and eosin to identify neighboring regions of myofibers both within fascicles and between fascicles (Fig. 9 A), and different subregions (1–6) were identified and then isolated by LCM and subjected to Illumina mRNA profiling. Histochemical staining for mitochondrial enzyme function showed the variable loss of mitochondrial oxidative phosphorylation activity typically seen in DMD biopsies (Fig. 9 B). There were variable degrees of inflammation and fibrosis, as detected by Masson’s trichrome and Gomori’s trichrome (Fig. 9, C and D). Unsupervised clustering of the microarrays from adjacent sections showed clear and consistent differences in expression profiles from the six regions, with replicates of profiles clustering together, yet each region showing a distinct pattern from other regions.
As we were unable to specifically stage the pathology of the human biopsy, the specific timing of the degeneration and regeneration cycles of the six regions could not be determined. Nonetheless, the profiles were consistent with asynchronous bouts of degeneration/regeneration leading to chronic inflammation and mitochondrial function loss associated with the murine 4-d spaced multiple injuries and the fibrosis associated with the 10-d multiple injuries.
Pharmacogenomics of glucocorticoids on the proposed intercellular signaling networks
Daily glucocorticoids are considered standard of care in DMD (Bushby et al., 2010). The mechanisms of action of glucocorticoids in DMD are thought to be through antiinflammatory mechanisms as in other chronic inflammatory states (Hoffman et al., 2012). However, extensive side effect profiles, such as muscle wasting typically associated with glucocorticoids, and differences in response between human and mouse models have complicated dissection of the molecular mechanism of action (Sali et al., 2012). Recently, we developed a novel dissociative steroid, VBP15, that is similarly antiinflammatory yet lacks most side effects associated with glucocorticoid response element–mediated transactivation subactivities (Heier et al., 2013). We hypothesized that both prednisone and VBP15 might serve to suppress inappropriate intercellular signaling in the regions between asynchronously injured and remodeling areas of muscle, particularly those associated with the proinflammatory 4-d asynchronous cross talk regions.
We tested the effects of drug treatment on the TGF-β–based asynchronous remodeling networks both in the mdx mouse model as well as the WT asynchronous remodeling experimental system. For the mdx model, we obtained frozen muscle samples from a previously described preclinical trial of prednisone and VBP15 (Heier et al., 2013). Drug treatments were started at 10 wk of age and continued to 26 wk of age (4 mo of treatment). Skeletal muscles (gastrocnemius) of mice treated with 5 mg/kg/d prednisone and 15 mg/kg/d VBP15 were tested by mRNA profiling (n = 5 per group) using Illumina bead arrays. As a baseline study, the human TGF-β–centered network, including HIF1A and S100A11 nodes, was queried in untreated mdx versus WT mice, and this showed strong up-regulation of the majority of network members (22/29; Fig. 10 A). A marked exception was HIF1A, a key ischemia-responsive gene; however, the rapid expression of HIF1A upon necropsy may lead to postmortem expression in WT animals. Comparing prednisone- and VBP15-treated mdx mice with untreated mdx mice showed both drugs to significantly down-regulate the majority of the pathway members (15/29 prednisone and 12/29 VBP15), including TGF-β ligand and HIF1A. These data suggest that a possible mechanism of action of both prednisone and VBP15 is through suppression of the TGF-β–related networks, leading to more successful regeneration.
We then used the same data analysis methods for the WT asynchronous remodeling experimental system, similarly determining the effect of drug treatment on expression of the TGF-β network. As a baseline, expression profiles of LCM samples (cross talk area vs. single injury sites) was performed for both 4- and 10-d asynchronous remodeling samples, with both showing up-regulation of the TGF-β network in the cross talk microenvironment (10/29 for 4 d and 16/29 for 10 d; Fig. 10). The 10-d cross talk microenvironment showed more extensive up-regulation of the TGF-β networks compared with 4 d, consistent with the findings presented earlier (Fig. 5). Drug treatment of the WT 4-d asynchronous mouse model was performed 1 d before first injury and then continuously during the 13 d after recovery from the second injection (17 d total), using either 5 mg/kg/d prednisone or 15 mg/kg/d VBP15. LCM and mRNA profiling of the cross talk microenvironment was compared with single injury microenvironment in the same muscle. This showed significant down-regulation by prednisone of TGF-β network members in the cross talk region (Fig. 10).
To validate these findings at the immunohistochemical (protein) level, the WT 4-d asynchronous model drug-treated muscles were immunostained for CD11b, CD74, and CD163 and compared with untreated mice. This analysis showed the expected high expression of each protein marker in the 4-d multiple injury sites (Fig. S5). However, daily administration of either prednisolone or VBP15 markedly reduced the immunopositivity for these antigens (Fig. S5). These data are consistent with both glucocorticoids and VBP15 serving to inhibit inappropriate cross talk in asynchronously remodeling muscle.
We used a large-scale data integration approach to query the molecular underpinnings of progressive failure of regeneration in the muscular dystrophies (Fig. S1). Using an initial mRNA profiling dataset of 12 disease groups (117 biopsies and 234 microarrays), we developed an iterative composite approach to identify clusters of transcripts that were associated with the more severe muscle disease classes (DMD, JMD, and LMNA; Fig. 1). Functional clustering of the diagnostic transcripts showed a 56-member TGF-β–centered network consistent with tissue fibrosis. Unsupervised clustering of an independent mRNA profiling dataset of 49 biopsies from four muscle disease groups (DMD, BMD, LGMD2I, and LGMD2B) and controls showed a strong correlation of the TGF-β–centered network with pathological severity rather than diagnostic criteria (Fig. 2). However, all members of this same 56-gene TGF-β–centered network were expressed during specific stages of normal muscle regeneration, using a 27–time point staged murine regeneration dataset (Fig. 3 A). These data analyses established the hypothesis that normal regeneration and progressive dystrophy involved the same proteins; the only distinction was one of timing and degree of expression. The TGF-β network was parsed into stage-specific subnetworks in normal muscle regeneration yet occurred within the same muscle at the same time in severely dystrophic muscle (Fig. 3 B). We hypothesized that asynchronous bouts of degeneration and regeneration in neighboring regions of muscle tissue could result in inappropriate cross talk between the regions.
To test this hypothesis, we developed an experimental model for asynchronous regeneration of muscle in neighboring microenvironments. An initial injury was induced in normal mouse muscle using notexin with the injection site labeled with tattoo dye. A second adjacent injection, labeled with a second tattoo dye, was subsequently made, at defined times (0–10 d). The mice were sacrificed 13 d after the second injection when regeneration would be expected to be largely complete (Fig. S2). Then, regions of the muscle containing both injection sites and regenerating myofibers in between the sites (cross talk microenvironment) were studied by LCM and in situ immunostaining methods. We found that increasing the time between injections led to a greater degree of abnormal collagen deposition (fibrosis; Fig. S3), indicating that the time period between injections was an important determinant of the response of the muscle to the damage. To better define the response of specific subregions of asynchronously remodeling muscle, we used LCM and mRNA profiling as well as immunostaining, histology, and enzyme chemistry experiments. Myofibers in between the two injury sites showed an expression pattern that seemed suspended within the developmental stage that separated the two injections, with the 4- versus 10-d between staged injections sharing little similarity with each other (Fig. 4 and Fig. S3, B and C). The myofibers between the sites of the 10-d repeat injections showed strong expression of transcripts associated with a profibrotic state, with high expression of collagens, and other connective tissue remodeling networks (Fig. 5). This showed considerable overlap of network members, such as TGF-β, COL1A1, COL4A1, COL4A2, and COL6A1, with the human muscle biopsy TGF-β network associated with fibrosis (Fig. 1 B) and was similar to the myotube maturation subgroup of normal staged (synchronous) muscle regeneration (Fig. 3). To assign statistical support to these visualizations of LCM data networks, we took specific network members (mRNAs) and defined expression changes (p-values and fold changes) in specific time points of regeneration, 4-d in-between LCM data and 10-d in-between LCM data (Table 4). The expression patterns were consistent with inappropriately prolonged and extensive connective tissue remodeling, in which signaling led to suspension of progression through the 10-d time window of normal regeneration. In other words, the tissue region between the injection sites staged 10 d apart remained fixed in the 9–10-d regeneration program and appeared unable to successfully mature past this developmental time point.
Our data showed that normal regeneration and failed regeneration involve the same networks, differences largely being caused by the synchrony (episodic) nature of normal regeneration and the asynchronous (repeated injury) nature of failed regeneration. With failed regeneration, regenerating myofibers are likely receiving incorrect temporal cues from their neighbors, and this inappropriate cross talk leads to inability to progress through the normal temporal stages of regeneration.
Glucocorticoids are derivatives of cortisol, the diurnal hormone that sets day/night cycles in most animals and plants. Pharmacological doses of glucocorticoids, typically given on a daily basis, are considered antiinflammatory and remain among the most effective drugs for chronic inflammatory states since their initial use in the 1940s. There is increasing recognition of the role of cortisol and clock genes in establishing cell cycle regulation throughout the body and that disturbances or blunting of diurnal cortisol fluctuations are associated with the proinflammatory state (Niimi et al., 2007; Holderfield and Hughes, 2008). We recently showed that primary bronchial epithelial cells from asthmatic patients are intrinsically dyssynchronous, are proinflammatory, and poorly repair wounds, yet this can be normalized with pulsed glucocorticoids (Freishtat et al., 2011; Alcala et al., 2014).
We reasoned that the effectiveness of glucocorticoids in DMD may be caused by resynchronization of remodeling, thereby rescuing the inappropriate cross talk described in this paper. To test this, we treated WT mice with prednisolone, VBP15 (a dissociative steroid), and vehicle/placebo and then tested for the drug regimens to rescue the inappropriate cross talk in our asynchronous regeneration model. Drug treatments of both mdx mice and the WT asynchronous model were tested (Fig. 10 and Fig. S5). As hypothesized, both prednisolone and VBP15 reduced evidence of cross talk in both mdx and the WT asynchronous models. These data suggest that part of the mechanism of action of both prednisone and VBP15 is through suppression of TGF-β networks induced by asynchronous remodeling.
There is an extensive literature supporting many of the findings in our study, although we believe that this study is the first to synthesize these observations into a single model. Specifically, transient TGF-β activity is essential for maintaining normal tissue repair and regeneration, whereas persistent TGF-β activity results in excessive fibrosis and, ultimately, scarring in skin and internal organs. In dystrophic muscle, TGF-β may affect tissue fibrosis by elevating collagen secretion by fibroblasts and converting myoblasts to fibroblasts. Furthermore, TGF-β can inhibit terminal myogenic differentiation and thus impair muscle repair process (Alexakis et al., 2007).
Our model may provide an explanation for the enigmatic pathology of the mdx mouse. Dystrophin-deficient mouse muscle shows normal histology in early weeks postpartum (0–3 wk) and then a stage of widespread necrosis followed by relatively successful regeneration. The mdx muscle in the range of 2–12 mo shows a mild dystrophy relative to human DMD muscle. Our model suggests that the widespread necrotic phase (2–3 wk) leads to well-synchronized regeneration; the severity of the initial necrotic phase may be beneficial in that it may avoid the poorly synchronized regeneration and pathology-associated signaling cross talk that occurs in human DMD muscle.
The regenerative asynchrony model may be applicable to chronic pathologies in nonmuscle tissues. Asthma involves chronic inflammation of the airway, where pulsed glucocorticoids are standard of care. In other studies, we have shown that normal epithelial wound repair is synchronous, that remodeling of asthmatic cells is asynchronous and proinflammatory, and that pulsed glucocorticoids restores regenerative synchrony and the nonproinflammatory state (Freishtat et al., 2011; Alcala et al., 2014). One could extend these models to chronic liver pathology, poor skin wound healing, or possibly all chronic inflammatory tissue states leading to fibrosis.
Materials and methods
Human patient muscle samples and mRNA profiles
Flash-frozen neuromuscular disease patient muscle biopsies were taken for diagnostic purposes. Normal control biopsies were obtained from baseline samples from exercise experiments in young adult volunteers.
Two datasets were analyzed. The first (test set) has been previously published (Bakay et al., 2006). This set contained Affymetrix mRNA profiles from 117 patient muscle biopsies using both HG-U133A and HG-U133B microarrays (n = 234 microarrays total; GEO accession no. GSE3307). The muscle disease groups were DMD (n = 10), amyotrophic lateral sclerosis (ALS; n = 9), acute quadriplegic myopathy (n = 5), BMD (n = 5), CAL (calpain III gene mutations; LGMD2A; n = 10), DYSF (LGMD2B; n = 10), ED-R (Emery–Dreifuss muscular dystrophy X linked; n = 4), ED-D (n = 4), FKRP (LGMD2I; n = 7), facioscapulohumeral dystrophy (FSH; n = 14), JDM (n = 21), and normal human skeletal muscle (NHM; n = 18). Biopsies were taken generally at the time of diagnosis and were typically from the vastus lateralis. The original published dataset analyzed 13 groups, but we removed the hereditary spastic paraplegia group from analyses reported here because of a high proportion of outlier transcripts and also removed four JDM samples that were not made public. To generate the mRNA profiles, frozen muscle biopsies were solubilized, total RNA was isolated, RNA was converted into biotinylated cRNA probes, hybridized to Affymetrix microarrays, and detected using fluorescent streptavidin, and microarrays were scanned. Quality control metrics used were as previously reported (Tumor Analysis Best Practices Working Group, 2004).
A second set of 49 human patient mRNA profiles was generated using HG-U133 Plus 2.0 microarrays, and this dataset is new to this study. These datasets contained profiles from 6 normal controls, 17 DMD (absence of dystrophin), 11 BMD (present but abnormal dystrophin), 7 LGMD2I (FKRP deficiency, a glycosylation defect), and 8 LGMD2B (DYSF). Patients had a broad range of ages, clinical severity of their disease, and histopathological findings, although all neuromuscular disease patients showed evidence of a dystrophic process (degeneration/regeneration of muscle fibers). Hematoxylin and eosin stains were performed on frozen sections, and the amount of fibrotic replacement (fibrosis) visually approximated for by the same evaluator (E.P. Hoffman) into normal, mild, moderate, or severe fibrosis categories.
Murine muscle regeneration time series
Staged degeneration/regeneration of adult WT mice was induced by intramuscular injection of cardiotoxin into the entire gastrocnemius muscle using a 10-needle manifold (Zhao and Hoffman, 2004). Mice were sacrificed at 27 time points after the injection, and muscles showing homogenous regeneration were selected for mRNA profiling (n = 2 per time point). Frozen muscle samples were solubilized, total RNA was isolated, RNA was converted to biotinylated cRNA, cRNAs were hybridized to MG-U74Av2 microarrays, and fluorescent images were scanned (GEO accession no. GSE469). Analyses and validations of this degeneration/regeneration time series have been previously published (Zhao et al., 2003, 2006; Zhao and Hoffman, 2004).
Biostatistical clustering methods
We developed a novel method of defining gene clusters within a large multigroup dataset, using the 12-group human muscle biopsy dataset (Fig. S1). Mapped image files from microarrays are interpreted using probe set algorithms that provide probe set normalization, chip normalization, and project normalization so that quantitative signals can be assigned to each probe set (transcript) represented on the microarray. We and others have shown that different probe set algorithms often lead to very different interpretations of the same microarray image files (Seo et al., 2004, 2006; Seo and Hoffman, 2006). To avoid loss of sensitivity for gene selection by using a single probe set algorithm, we simultaneously interpreted the muscle biopsy data using three distinct probe set algorithms (dCHIP, MAS 5.0, and PLIER; Fig. S1). Both A and B microarrays were used, leading to two images with distinct probe sets representing each sample.
Two-way clustering on the rows (genes) and columns (muscular disorders) of the data matrix was performed using the VISDA (Visual Statistical Data Analyzer) software to cluster genes (Fig. S1; Wang et al., 2000, 2007; Zhu et al., 2008). Within each gene cluster, the clustering of the muscular disorders was performed using hierarchical clustering with correlation as the similarity measure after a preprocessing step in which all samples within one muscle group were taken as an average.
This led to a series of probe set algorithm-specific gene clusters that grouped the muscular dystrophies into larger groups based upon the specific microarray used (U133A chip or U133B chip). Thus, PA cluster i refers to PLIER algorithm U133A chip dataset, MA cluster i refers to MAS 5.0 algorithm U133A chip dataset, and dA cluster i refers to dCHIP algorithm U133A chip dataset (Fig. S1 and formulas in this section). Each gene cluster was then tested for gene–gene and gene–protein interactions between cluster members using IPA. The top-ranked networks were then retained for each cluster i (PA network i, MA network i, and dA network i; Fig. S1). We then merged networks and assigned cumulative scores. To accomplish this, we implemented programs to query overlaps among networks and then derived a total score for each network based on its overlaps with all other networks of each gene-disorder pattern. The probability for two networks, S1 and S2 with n1 and n2 genes, respectively, having m overlapped genes can be calculated as
given that S1, S2, and the overlapped genes belong to the global network with N genes (N = 26,400, the total number of mammalian genes in Ingenuity Pathway Knowledge base), and C represents a mathematical operation of combination. C of N and m gives the number of subsets of m distinct elements of a set of N elements. Because this probability indicates the likelihood of the overlap of the two networks occurring randomly, −log10(p(m)) may be used as a score for specifying the reliability of the overlap. In addition, the IPA dynamically computes relevant networks by focusing on the genes of interest and delivers a set of scored networks for each gene cluster. A composite score for each network was then calculated by summing up the IPA score and the total score based on network overlaps; the networks were ranked based on their composite scores.
Cross-species query of time series data
Once a gene cluster of interest was selected, probe sets were tested for recognition of canonical transcripts (e.g., 3′ end of well-documented transcript units), and then, these human transcript units were mapped to corresponding rat and mouse transcript units using Affymetrix NetAffx Analysis Center. The Affymetrix probe set IDs corresponding to the transcript units were then queried against the mouse time series data. All mouse image files were analyzed using the PLIER probe set algorithm. Candidate gene clusters were visualized in mouse datasets using heat maps (hierarchical clustering).
All differentially regulated transcripts visualized were tested for statistical significance using t test, without correction for multiple testing. For time series data, adjacent time points were grouped and statistically compared with a relevant control region of the time series.
Asynchronous regeneration protocol
WT BL57 mice of 4–8 wk of age were used for all asynchronous remodeling experiments. To induce asynchronous bouts of regeneration, a 20-guage needle was filled with 10 µl of 5-µg/ml notexin and then dipped in tattoo dye (blue for first injection; red for second injection). Both gastrocnemii were surgically exposed, and the initial injection was made distally to proximally along the axis of the muscle. The wound was closed for the number of days indicated, and then, the second injection was performed after the time annotated in each experiment. The initial injection site could be visualized in the muscle as a result of the tattoo dye, and the second injection was initiated adjacent the first site. Mice were euthanized 13 d after the second injury, and gastrocnemii were collected.
For drug treatments with 5 mg/kg/d prednisolone, 15 mg/kg/d VBP15, or vehicle, the drug was delivered in sugar syrup orally. We had five mice in each group. Drug treatments were started 24 h before the initial notexin injection and then continued daily until sacrifice.
Van Gieson measures of fibrotic tissue
Reinjured muscle sections were stained by the Van Gieson method (n = 4 muscles per group). Images were taken and analyzed with ImageJ (National Institutes of Health) to measure the extent of the connective tissue. We assayed the amount of connective tissue per field of tissue by measuring the number of pixels occupied by red staining and then measuring the number of pixels covered by the muscle tissue. The data are then represented as the percentage of muscle that is collagenous connective tissue.
Laser capture microdissection
Flash-frozen muscle biopsies or specimens were sectioned at 10-µm thickness on RNase-free glass slides (CellCut; Molecular Machines & Industries, Inc.). Slides were immediately placed in 99% ethanol to fix. Rapid hematoxylin and eosin stain using RNase-free solutions (Molecular Machines & Industries, Inc.) was then performed followed by laser capture microdissection (CellCut Plus instrument) to identify and cut the areas of interest. We used tubes with an adhesive cap to collect LCM-isolated regions (Molecular Machines & Industries, Inc.).
In brief, 15 consecutive sections for mice at 4 and 10 d were subjected to LCM to isolate each region on each muscle, and cRNA amplification was performed from the combined sections. 3–4 of the 15 LCM sections (per area per mouse) were successful in mRNA profiling, and these showed excellent chip-based clustering. Although the LCM was performed on only one mouse for 4 d and one mouse for 10 d (in part because of the need to identify closely spaced tattoo dye injection sites), additional mice were used for immunostaining validations.
LCM mRNA profiling
RNA isolation was performed using PicoPure kit (Applied Biosystems), which includes on-column DNase digestion. Two rounds of T7-mediated cRNA amplification were performed (RiboAmp HS kit [Applied Biosystems] and then TotalPrep-96 RNA Amplification kit [Life Technologies]). The mRNA profiling on asynchronous mouse LCM experiments was performed in triplicate on serial sections. Murine WG-6 v2.0 Expression BeadChips (Illumina) containing 45,200 transcript probes or HumanHT-12 v4 Expression BeadChips (Illumina) containing 47,000 transcript probes were used. Results of expression arrays were read using the HiScan SQ System (Illumina). Statistical analyses (P ≤ 0.05 with correction for multiple testing [False Discovery Rate]) and a fold change ≤1.5 or ≥1.5) were performed to identify dysregulated genes in the areas of interest using Partek Genomics Suite software (Partek, Inc.). IPA was used to query protein–protein and gene–protein interactions.
TaqMan gene expression assay
The RiboAmp HS kit was used to perform two rounds of amplification on total RNA isolated from LCM samples (n = 3 muscles per group) and convert amplified RNA to cDNA. Quantitative gene expression analysis was performed by TaqMan-based qRT-PCR on a real-time PCR system (7900 HT Fast; Applied Biosystems). We used TaqMan Gene Expression kits obtained from Applied Biosystems to run qRT-PCR: Cox6b1 (Mm00824357_m1), Ndufs6 (Mm02529639_u1), Ndufb5 (Mm00452592_m1), Ndufb6 (Mm01208591_g1), and HPRT (hypoxanthine guanine phosphoribosyltransferase; Mm01318747_g1). The expression of HPRT was analyzed as endogenous control simultaneously, and the target gene expression in each sample normalized to HPRT. Thermal cycling conditions were 95°C for 10 min, 50 cycles at 95°C for 15 s, and 60°C for 1 min.
Monoclonal rat anti–mouse CD11b (eBioscience), mouse anti–human CD74 (Abcam), mouse anti–human macrophage surface antigen (AM-3K/CD163), mouse anti–human mannose receptor (CD206; Abcam), and mouse anti–human MMP9 (Santa Cruz Biotechnology, Inc.) monoclonal primary antibodies were used to perform immunofluorescence assays. The Mouse on Mouse kit (Vector Laboratories) was used to reduce background staining when using mouse monoclonal antibodies on mouse tissue sections (n = 3 muscles per group). Goat anti–mouse (Invitrogen) and goat anti–rat (Invitrogen) secondary antibodies conjugated with Cy3 or Alexa Fluor 488 were used. Nuclei were stained with DAPI (OCHEM). The images were acquired using a fluorescence microscope (BX61; Olympus).
NADH, SDH, COX, costaining of NADH and COX, Masson’s trichrome, and Gomori’s modified trichrome were performed on 10-µm biopsy frozen sections to evaluate histology and mitochondrial enzymatic activities (Dubowitz and Sewry, 2007). Murine gastrocnemius 10-µm sections (n = 3 muscles per group) were stained following a standard Van Gieson protocol (Dubowitz and Sewry, 2007). A minimum of nine nonoverlapping images were analyzed in ImageJ for the percentage of connective tissue per field of muscle tissue.
Drug treatment of
mdx and WT asynchronous remodeling mice
Drug-treated mdx mice were from the previously reported preclinical trial comparing 4 mo of 5 mg/kg/d prednisone treatment to 15 mg/kg/d VBP15 treatment (Heier et al., 2013). Frozen archival muscle samples were obtained for WT, mdx untreated, and the two drug-treated groups (prednisone and VBP15). RNA was isolated, and mRNA profiling was performed using Illumina bead arrays, according to the manufacturer’s protocol (n = 5 per group). Bioinformatics analyses were performed using TGF-β pathway candidates. Drug treatment of the WT asynchronous remodeling mice (n = 5 per group) was performed using the same drugs and treatment regimens (5 mg/kg/d prednisone; 15 mg/kg/d VBP15), with treatment starting the day before the first surgery until the day before euthanasia.
Microscopic image acquisition
All images were acquired using a microscope (BX61; Olympus), with U Plan S Apochromat 20×/0.75 NA dry objective, at room temperature. Photography was performed using a camera (DP11; Olympus) with DP controller 126.96.36.1996 and DP Manager 188.8.131.52 software.
Online supplemental material
Fig. S1 shows a flow chart of the data integration and biochemical network analysis on the 12-group human muscle transcriptional profiling data. Fig. S2 shows LCM methods and chip clustering of mRNA profiles of LCM regions of multiply injured muscle. Fig. S3 shows quantitative analysis of endomysial fibrosis in asynchronous remodeled muscle samples. Fig. S4 shows histochemical stains for mitochondrial activity showing a relative loss of mitochondrial oxidative phosphorylation in the area between injuries spaced 4 d apart. Fig. S5 shows treatment of asynchronous remodeling mouse muscle with prednisone or VBP15 suppresses the proinflammatory proteins in between injuries spaced by 4 d.
This work was supported by the US Department of Defense (W81XWH-05-0334), National Institutes of Health (R01 NS029525, U54 HD053177, and R24 HD050846), Children’s National Medical Center Board of Visitors, and the Muscular Dystrophy Association USA. The development of VBP15 is supported by the Congressionally Directed Medical Research Programs of the Department of Defense, Muscular Dystrophy Association Venture Philanthropy program, and the National Institutes of Health National Center for Advancing Translational Sciences, Therapeutics for Rare and Neglected Diseases program.
Competing financial interests: E.P. Hoffman and K. Nagaraju are coinventors of patent WO 2011127048 A3 (Non-hormonal steroid modulators of NF-κB for treatment of disease), and this has been assigned to ReveraGen Biopharma. E.P. Hoffman and K. Nagaraju own founder shares of ReveraGen Biopharma and are members of the Board of Directors. The authors declare no further competing financial interests.
Author contributions: S. Dadgar contributed to the experimental design, data generation and interpretation, and manuscript writing. Z. Wang contributed to the bioinformatics methods development and conduct, figures, and manuscript writing. H. Johnston-Carey contributed to the experimental design, data generation, and writing. A. Kesari contributed to the human muscle biopsy expression profiling design and conduct. K. Nagaraju contributed to the experimental design and data interpretation. Y.-W. Chen contributed to the human muscle biopsy expression profiling design, conduct, and data interpretation. D.A. Hill contributed to the immunohistochemical and histological stains and data interpretation. T. Partridge contributed to the experimental design and data interpretation. M. Giri contributed to the bioinformatics methods and data interpretation. J. Nazarian contributed to the experimental design and data interpretation. R. Freishtat to the contributed to the experimental design and data interpretation. J. Xuan contributed to the bioinformatics methods development. Y. Wang contributed to the bioinformatics methods development. E.P. Hoffman contributed to the experimental design and interpretation and manuscript writing.
amyotrophic lateral sclerosis
Becker muscular dystrophy
cytochrome c oxidase
Duchenne muscular dystrophy
laser capture microscopy
normal human skeletal muscle
S. Dadgar, Z. Wang, and H. Johnston contributed equally to this paper.