Collagen VI myopathies are genetic disorders caused by mutations in collagen 6 A1, A2 and A3 genes, ranging from the severe Ullrich congenital muscular dystrophy to the milder Bethlem myopathy, which is recapitulated by collagen-VI-null (Col6a1−/−) mice. Abnormalities in mitochondria and autophagic pathway have been proposed as pathogenic causes of collagen VI myopathies, but the link between collagen VI defects and these metabolic circuits remains unknown. To unravel the expression profiling perturbation in muscles with collagen VI myopathies, we performed a deep RNA profiling in both Col6a1−/− mice and patients with collagen VI pathology. The interactome map identified common pathways suggesting a previously undetected connection between circadian genes and collagen VI pathology. Intriguingly, Bmal1−/− (also known as Arntl) mice, a well-characterized model displaying arrhythmic circadian rhythms, showed profound deregulation of the collagen VI pathway and of autophagy-related genes. The involvement of circadian rhythms in collagen VI myopathies is new and links autophagy and mitochondrial abnormalities. It also opens new avenues for therapies of hereditary myopathies to modulate the molecular clock or potential gene–environment interactions that might modify muscle damage pathogenesis.
Collagen VI (ColVI) is a major component of the extracellular matrix (ECM) crucial for muscle cell adhesion, stability and regeneration. ColVI forms a microfibrillar network involved in maintaining tissue integrity through a structural link between different components of connective-tissue basement membranes and cells (Allamand et al., 2011). ColVI is a heterotrimeric molecule composed of three α-chains encoded by distinct genes (COL6A1, COL6A2 and COL6A3). The three α-chains first assemble in equimolar association into a triple helical monomer, and successively form dimers and tetramers stabilized by disulfide bonds. The secreted tetramers give rise to a network of microfilaments that bridge the sarcolemma with the extracellular connective tissue (Bonaldo et al., 1990). Mutations occurring in the three ColVI genes cause two major disease phenotypes: a mild form, Bethlem myopathy, and a severe form, Ullrich congenital muscular dystrophy (UCMD) (Gualandi et al., 2012). Three additional ColVI genes have recently been described, COL6A4, COL6A5 and COL6A6, and code for the α4(VI), α5(VI) and α6(VI) protein chains, respectively. However, only α6(VI) is highly expressed in skeletal muscle, although not associated with basement membranes (Tagliavini et al., 2014a) and none of these additional ColVI genes are known as disease genes in humans.
ColVI myopathies are well-defined conditions, characterized by a broad spectrum of clinical features, ranging from prenatal onset with decreased fetal movements; hypotonia or torticollis, early-childhood onset, characterized by delayed motor milestones, muscle weakness and contractures; and adult onset (4th–6th decade), characterized by proximal weakness and Achilles tendon or long finger flexor contractures (De Visser et al., 2004). Owing to the progressive nature of the condition, 75% of patients over age 50 years need assistance for outdoor mobility (Jöbsis et al., 1999).
The severe phenotype with early onset, UCMD, is characterized by severe muscle weakness, striking hyperlaxity of distal joints, and proximal joint contractures. In contrast to Bethlem myopathy, UCMD is marked by an early loss of ambulation and a rapid disease progression leading to early death caused by respiratory failure (Maraldi et al., 2009). In addition, two other conditions have been associated with ColVI myopathy: autosomal-dominant limb-girdle muscular dystrophy (LGMD), and autosomal-recessive myosclerosis myopathy (Merlini et al., 2008; Scacheri et al., 2002). Recently, the COL12A1 gene has been shown to be mutated in patients with a Bethlem-myopathy-like phenotype (Hicks et al., 2014).
Despite this broad clinical spectrum, both muscle function and respiratory performance evaluation can address the diagnosis. Nevertheless, UCMD and Bethlem myopathy can be inherited accordingly to both dominant and recessive models, and generally, neither the type of mutation (missense, nonsense, splicing insertion or deletion), nor the effect of the mutation on the protein structure or function allows precise discrimination between two phenotypes (Allamand et al., 2011; Maraldi et al., 2009; Bönnemann, 2011). This suggests that the genetic background of patients might influence the expression of the mutated genes, therefore impacting on the resulting phenotype.
The Col6a1−/− mouse is a well studied animal model for ColVI myopathies. These mice carry a targeted inactivation of the α1(VI) chain and exhibit an early onset myopathic phenotype that strongly resembles Bethlem myopathy (Bonaldo et al., 1998). Muscle fibers of Col6a1−/− mice show loss of contractile strength and a number of structural defects in the hind limb muscle and especially in the diaphragm. Although the pathology associated with collagen VI myopathy is nonspecific and variable, it has been widely documented (Grumati et al., 2010; Urciuolo et al., 2013). Both muscles from mice and patients have been found to show mitochondrial dysfunction accompanied by ultrastructural alterations, increased spontaneous apoptosis and abnormal propensity of the permeability transition pore (PTP) to open (Irwin et al., 2003; Bernardi and Bonaldo, 2008; Tagliavini et al., 2014b).
Col6a1−/− mice and UCMD or Bethlem myopathy patients also show defective autophagy (Grumati et al., 2011). This suggests a mechanism of disease pathogenesis that might be in part due to defective regulation of autophagy, which leads to accumulation of abnormal mitochondria and sarcoplasmic reticulum. Moreover, in Col6a1−/− mice, changes in metabolic proteins resulting in decreased glycolysis and affecting the tricarboxylic acid (TCA) cycle fluxes lead to a different fate of α-ketoglutarate triggering lipotoxicity. The metabolic changes are associated with changes of proteins involved in mechanotransduction at the myotendineous junction, costameric or sarcomeric level, and in energy metabolism (De Palma et al., 2013). Recently, a strong depletion of the α6(VI) chain in skeletal muscle has been reported as potential biomarker of ColVI myopathies (Tagliavini et al., 2014a). In spite of these recent achievements, a full understanding of ColVI disease pathogenic mechanisms remains unknown.
Gene expression profiling and RNA sequencing are powerful tools to explore the transcriptional alterations in disease tissues (Sánchez-Pla et al., 2012). It brought great advances in cancer research and disclosed deregulated transcripts that are now being considered as prognostic or therapeutic biomarkers of neoplastic diseases (Sevov et al., 2012; Campo, 2013). Considering the availability of novel therapeutic possibilities for rare muscle disorders (Aartsma-Rus and Muntoni, 2013), finding prognostic biomarkers for disease severity or drug response will be of great benefit for patient care and treatment (Ferlini et al., 2013; Scotton et al., 2014). In this work, we explored RNA profiling of Col6a1−/− mice and Bethlem myopathy and UCMD patients to identify transcriptional patterns correlating with disease state. Such patterns can identify specific pathway alterations and pathophysiological biomarkers. This unbiased approach revealed marked changes in the expression of genes involved in the circadian rhythm pathway. We further supported these findings by protein studies performed both in mice and in patients. In addition, RNA and protein studies in the Bmal1−/− (also known as Arntl) mice, a known animal model of circadian rhythm disruption, showed profound deregulation of Col6a6 as well as of the autophagy-related genes such as p62/SQSTM1 and AMPK. Although it has been previously reported that disruption of circadian rhythms are linked to skeletal muscle remodeling, function, performance and aging, the involvement of the molecular clock genes in hereditary myopathy as ColVI disease is a new finding. These findings highlight new potential target genes involved in muscle damage possibly addressing new therapies.
Overall gene expression profile in Col6a1−/− muscles reveals changes in genes associated with muscle function
Gene expression array profiles of wild-type (WT) and Col6a1−/− mice were compared in three muscle types: diaphragm, tibialis anterior and gastrocnemius. To determine genes with statistically significant differential expression (P<0.05 in limma package) we analyzed each muscle type (either separately and in pools) and determined average expression change over all muscle types. The degree of transcriptional deregulation was different among Col6a1−/− muscles, with 11,170 differentially expressed genes in tibialis anterior, 6790 in diaphragm and 3759 genes in gastrocnemius. To further support the finding of different tissue susceptibility to ColVI deficiency, we performed sample clustering of expression profiles from each muscle type for Col6a1−/− and in WT mice. The resulting dendrogram is shown in Fig. S1A. We found that the expression profiles of Col6a1−/− and WT gastrocnemius were clustered together and therefore similar. Indeed this supports the existing knowledge about muscle involvement in the Col6a1−/− mouse model, where gastrocnemius is very mildly affected, whereas tibialis and even more diaphragm are severely affected (De Palma et al., 2013). In contrast, diaphragm and tibialis anterior Col6a1−/− muscles clustered together and were clearly separate from their WT counterparts, suggesting that ColVI deficiency severely and similarly affects these muscle types (Fig. S1).
Expression of genes associated with muscle development and function was altered in Col6a1−/− mice (Table S1). The identified changes involved several different muscle biology circuits, including structure, metabolism, transcription and regulation. There was variability between tested muscle types. A total of 60% (26 out of 44) of muscle-specific genes changed their expression in an opposite direction in tibialis anterior versus diaphragm, 23% in tibialis anterior versus gastrocnemius (7 out of 30) and 67% in gastrocnemius versus diaphragm (12 of 18). Tibialis anterior showed differential expression of myosin, troponin and calsequestrin isoforms (predominantly cardiac and slow), and upregulation of Pax3. In diaphragm and tibialis anterior, the differentially expressed genes included almost all members of the dystroglycan complex (dystrophin, dystrobrevin, dystroglycan 1, syntrophin, α-, β-, γ-sarcoglycans) and several transcription factors important for muscle-related processes, such as MyoD1, MyoG, MYF6, NFATc2 and Pax7. Embryonic isoforms of myosin were upregulated in the diaphragm. Pax7, a marker of satellite cell number, which can be regarded as a sign of muscle regeneration, was also upregulated. This muscle type variability in the transcription profile is not surprising given that these three muscles have different metabolic properties and different fiber type composition (Bonaldo et al., 1998).
Gene ontology groups and key expression regulators analyses show deregulation of the circadian rhythmic process in Col6a1−/− mice
Among the differentially expressed genes between Col6a1−/− and WT mice, 479 changed expression in the same direction in all tissues with P<0.05. We used Gene Set Enrichment analysis (GSEA) to identify gene ontology groups enriched with such genes. The top ten gene ontology groups are listed in Fig. S2. Identified gene ontology groups were broad categories and, therefore, do not specifically address any definite biological effects occurring in Col6a1−/− mice. We also performed gene ontology analysis of 334 genes differentially expressed in the ‘pooled’ data from all three muscle groups compared with ‘pooled’ WT muscles. This GSEA identified the gene ontology groups shown in Fig. 1. One of the top ten gene ontology groups was ‘rhythmic process’, which contains genes regulating circadian rhythms, as well as known genes that are regulated downstream of the molecular clock, such as DBP, and all were downregulated in Col6a1−/− muscles.
To identify key regulators that change their activity in Col6a1−/− muscles, we used sub-network enrichment analysis (SNEA) in Pathway Studio. Our first SNEA was performed using the list of genes differentially expressed in one direction in all muscle types as an input. Regulators identified by this approach are shown in Fig. 2. Genes changing their expression in one direction in all muscle types are primarily regulated by inflammatory cytokines and their receptors. In good agreement with the gene ontology analysis, SNEA of the ‘pooled’ dataset found BMAL1 (ARNTL), CLOCK, PER1 and PER3 regulators that belong to cellular circadian rhythm pathway (Fig. 2; Table S2). We then explored the intersection between the two SNEA analyses and found consistent results (Fig. 2; Table S2). Indeed, we found that CLOCK, BMAL1, GFI1B, IL2, CD28, PIM3 and inflammatory cytokines were differentially expressed genes that were also identified as being key regulators in all three comparisons.
To validate the 46 key regulators selected by SNEA, we studied the differential expression of these transcripts using TaqMan low-density arrays (TLDA). By comparing the microarray experiment data to those of the TLDA cards, we observed the same trend of fold change in 37 out of 46 genes (80%) of the diaphragm and tibialis anterior, whereas in the gastrocnemius we found 32 out of 46 genes (70%) with these features (Fig. S3A,B).
Both in the diaphragm and tibialis anterior and, to a lesser extent in gastrocnemius, the TLDA and microarray data for regulators of the rhythmic process (BMAL1, ATF5, CLOCK, DBP, EGR1, FKBP5, PER1, PER2 and PER3) and transcription factors important for muscle-related processes (MyoD1, MyoG, MYF6, PAX7) were in good agreement. This confirmed the deregulation of these two pathways in the Col6a1−/− mice (Fig. S3B).
Western blotting revealed that CLOCK levels were significantly (Student's t-test, n=3, P<0.05) increased in diaphragm and tibialis anterior muscles of Col6a1−/− mice compared with their respective WT controls, whereas they remained unchanged in gastrocnemius. MAT2A protein levels did not change significantly in the diaphragm, whereas in tibialis anterior and gastrocnemius muscles they were, respectively, increased and decreased (Fig. 3A).
Abnormal collagen VI expression in Bmal1−/− mice, a model for circadian rhythm deregulation
There are multiple genetic mouse models in which molecular clock factors have been knocked out, but of these models only the Bmal1−/− mice exhibit complete loss of behavioral rhythm under normal housing conditions (Bunger et al., 2000). We adopted this last model in order to have a homogeneous phenotype, with complete failure of circadian control. Western blot analysis of gastrocnemius muscles of Bmal1−/− mice revealed significant changes in the levels of ColVI chains, with decreased α6(VI) and increased α3(VI) levels. Interestingly, downregulation of COL6A6 was also recently observed by us in UCMD and Bethlem myopathy muscles (Tagliavini et al., 2014a). In addition to ColVI, components of the autophagy pathway were also affected in the muscle of Bmal1−/− mice with increased LC3A (also known as MAP1LC3A) and decreased p62 (also known as SQSTM1) levels, whereas Beclin 1 levels were unchanged. We also found an increase in the activation status of AMPK, as defined by a higher phosphorylated AMPK to total AMPK ratio; however, this was driven in part by a significant decrease in the levels of total AMPK (Fig. 3B).
Circadian genes show deregulation in muscle from UCMD and Bethlem myopathy patients
Gene expression profiles in three patients (one with UCMD and two with Bethlem myopathy) were compared to healthy controls by RNAseq, and we identified several altered transcripts, which we grouped into gene ontology categories (Fig. 4A). Genes that were deregulated in muscles from UCMD and Bethlem myopathy patients are reported in Table 1. In all three patients, we identified seven consistently upregulated genes belonging to transcription factor, metabolism process and cellular stress response categories. The 24 consistently downregulated genes were involved in signaling, metabolism, immunity and cell adhesion, axon guidance and circadian rhythm, including DNA repair gene ontology categories (Fig. 4B). Consistent with a link to circadian rhythms, we found that eight of the genes listed in Table 1 are known to be circadian in mouse skeletal muscle and/or changed in expression in the Bmal1−/− mouse.
In order to find gene ontology groups deregulated in different samples, we performed GSEA analysis (P-value 0.05 cut-off) for each patient separately. Bi-clustering was performed by using the Jaccard-based distance between gene ontology groups (based on the number of common genes in groups). This method allowed finding of common deregulated ‘patterns’ in different patients, even though the GSEA method shows different gene ontology groups. We found genes coding proteins belonging to postsynaptic membrane, cell junction and synapse, regulation of ion transmembrane transport, protein homo-oligomerization or protein-specific binding, cell-surface-receptor-linked signaling pathways and DNA repair or response to DNA damage, including circadian rhythm genes (Fig. 4A,B).
Considering the evidence from both mice and human studies for the deregulation of circadian genes in muscles with ColVI defects, we performed an real-time PCR expression analysis of CLOCK and MAT2A genes in another 11 UCMD and 12 Bethlem myopathy patients. We found that the CLOCK gene was significantly downregulated (P<0.0001) in all UCMD muscles. The MAT2A transcript showed a milder downregulation, although not statistically significant (P=0.0753), but only in UCMD patients (Fig. 3C). Bethlem myopathy patients did not show any significant deregulation in the analyzed circadian genes. Western blotting in four UCMD and three Bethlem myopathy patient muscles showed significantly increased CLOCK protein levels (P<0.05, ANOVA plus Tukey, n=3) in UCMD patients only, whereas MAT2A protein levels were increased in Bethlem myopathy patients only (Fig. 3D).
ColVI myopathies have been both clinically and genetically defined, and alterations in the apoptosis and autophagy circuits have been reported. Nevertheless, the link between mutations of COL6A1–COL6A3 genes and the muscle dysfunction seen in ColVI myopathy still remains elusive. We therefore adopted RNA profiling to study muscle from both Col6a1−/− mice and Bethlem myopathy or UCMD patient biopsies, searching for unique transcriptional patterns. The idea was to detect specific signatures that might enhance our understanding of the pathogenesis of ColVI myopathies. In order to strengthen our findings, we also evaluated protein expression of the transcripts that we found to be consistently deregulated in both mouse and patient muscle tissue.
Col6a1−/− and Bmal1−/− mice shows connection between ColVI myopathy and circadian genes
To our knowledge, this is the first attempt to study Col6a1−/− mice by RNA profiling. Proteomic studies were performed by Irwin et al. (2003) in Col6a1−/− mice, revealing the diaphragm as the most severely affected muscle, having a deeply subverted ultrastructure and evident loss of contractile strength, accompanied by changes in proteins involved in mechanotransduction. Our pooled data identified 334 differentially expressed genes in Col6a1−/− versus WT, showing consistent deregulation of the genes involved in metabolic processes and regulation of transcription, known to be affected in myopathies, as well as in those of the circadian CLOCK gene. This is an entirely new observation.
We show major transcriptional changes in the diaphragm and tibialis anterior, accordingly to the proteomic findings (De Palma et al., 2013). All muscles shared similar alterations in genes involved in the inflammatory response, cell proliferation, transcriptional regulation and metabolism, but differed substantially in the direction of expression change in a number of genes related to muscle functions, including transcription, metabolism, structure and regulation – a not entirely unexpected finding considering that these three muscles have different morphogenetic origins and functional requirements (Eberhard and Jockusch, 2004; Merrell and Kardon, 2013).
The CLOCK gene deregulation, we also confirmed at the protein level in Col6a1−/− mice, suggests a previously unreported relationship between circadian rhythms, the molecular clock and myopathy. The Col6a1−/− model shows both muscle and osteotendinous junction damage (Izu et al., 2011,, 2012), reflecting defective costamerogenesis and altered mechanotransduction (De Palma et al., 2013). Specific activation of both myocardin-related transcription factors (MRTFs) and ternary complex factors (TCFs), whose targets are regulators of cytoskeleton dynamics and mechanosensing, have potential to reset the molecular clock through regulation of PER1, PER2, BMAL1 and other CLOCK genes (Esnault et al., 2014). This first demonstration of a link between the ECM and the profound alterations we detected in the circadian rhythm or molecular clock pathway (BMAL1, ATF5, CLOCK, DBP, EGR1, FKBP5, PER1, PER2 and PER3 genes) supports the known data demonstrating that physical activity can regulate circadian gene expression in skeletal muscle (Wolff and Esser, 2012). In addition, providing a further bridge between protein and RNA data, intersection between the individual muscle groups and the pooled transcription data, to identify key regulators, showed that inflammatory cytokines and their receptors, known to be regulated by PER1–PER3 genes, are significantly underexpressed. The connection between ColVI and circadian genes is also strongly supported by our protein findings in Bmal1−/− mice, which exhibit significant changes in the levels of ColVI chains as well as alterations of the autophagy pathway (Fig. 3B). Given that our expression analyses were only performed at one time point, our results can only allow us to conclude that molecular clock pathway is disrupted. Further studies are needed to determine whether the clock factors are phase shifted, damped or arrhythmic.
RNA profiling of ColVI patients shows that there is a profound deregulation in ECM-mediated structural muscle remodeling and circadian genes
In the three ColVI patients (Bethlem myopathy, severe Bethlem myopathy and UCMD) studied, the highly efficient RNAseq analysis revealed many changes in transcription profile, as expected for this progressive, chronic and severe myopathy. The major changes observed, shared by all patients, were underexpression or overexpression in proteins related to membrane, cell junction or synapse, ion transmembrane transport, protein homo-oligomerization or protein-specific binding, the cell-surface-receptor-linked signaling pathway and circadian gene pathways (Fig. 4B). The overexpressed genes we detected (ZNF384, ZNF526, SCP2, SUOX, KIF1B, PCSK6 and SESN1) are all involved in the control of gene expression and metabolic pathways. As reported in UCMD muscles (Paco et al., 2013), increased expression of these transcripts is associated with muscle regeneration defects and inflammatory signatures. Our RNAseq data confirm that this myopathy is associated with structural alterations coupled with either signaling or regeneration failure (Bönnemann, 2011). The 24 downregulated genes we identified mainly involved transcripts associated with cell–cell communication or adhesion, and signaling. Deregulation of transcription, homeostasis and the cell response was also evident. This RNAseq data combines with that obtained by Paco et al. (2013) using array technology to indicate that ColVI deficiency drives both structural remodeling of muscle tissue, mediated by the ECM, adhesion, cell signaling and metabolism genes, and a massive immune response. The central role of metabolism in muscle remodeling has also been confirmed by De Palma et al. (2014), who described changes in the unfolded protein response coupled with metabolic dysregulation, autophagic impairment and mechanotransduction signaling alterations.
Most intriguingly, and in line with the mouse data, our RNAseq analysis also revealed that circadian-related CRY2 and MAT2A genes were profoundly deregulated. CRY2 is a well-established negative regulator of the molecular clock mechanism that underlies circadian rhythms, and MAT2A has a link with circadian rhythms through regulation of the enzyme required for melatonin synthesis, and mediates photoneural-circadian regulation of the pineal gland (Kim et al., 2005). Using array technology, Paco et al. (2013) also identified significant deregulation of the CLOCK pathway gene EGR1, which also supports circadian rhythm involvement in the pathogenesis of ColVI in both mice and patients. This finding prompted us to profile both MAT2A and CLOCK in larger patient cohorts, and the CLOCK gene was significantly downregulated (P<0.0001) in all UCMD patients, the more severe form of ColVI myopathy. All these data support that circadian genes, and CLOCK in particular, have a role in ColVI myopathy pathogenesis.
Circadian genes and hereditary myopathies – a new perspective
Circadian rhythms are approximate 24-h oscillations in physiology and behavior driven by the molecular clock, a transcription–translation feedback mechanism. A link between circadian genes and skeletal muscle has been identified in both animal models and in humans (Schroder and Esser, 2013). In particular, expression profiling has identified genes expressed in a circadian manner in the gastrocnemius, soleus and tibialis anterior muscles of mice (McCarthy et al., 2007; Dyar et al., 2013; Hodge et al., 2015). Analysis of human muscle samples taken 12 h apart have also shown changes in clock gene expression, consistent with the findings in mice (Zambon et al., 2003). Indeed, Zhang et al. (2012) found that the molecular clock factors BMAL1 and CLOCK transcriptionally regulate expression of MyoD1 in adult muscle, thereby suggesting a link between the daily clock cycle and periods of muscle repair and maintenance. We add to this growing body of evidence by demonstrating, by RNA and protein analysis, a new link between ColVI myopathy and deregulation of CLOCK circadian gene expression in both mouse and human tissues.
In mice, the protein data fully support the transcription results. Altered expression of CLOCK seems to be a marker of ColVI disease severity, being upregulated in the most profoundly affected muscles (tibialis anterior and diaphragm). In this way, CLOCK involvement appears to provide the missing link between autophagic pathway alterations and the physiopathogenesis of ColVI myopathies (Fig. 5). Supporting this link, many mitochondrial genes show a circadian oscillation in expression and are significantly downregulated in the muscle of the CLOCK mutant mice (McCarthy et al., 2007); in addition, the circadian animal models (both ClockΔ19 and Bmal1−/−) have a dramatic reduction (40% less) in the mitochondrial volume. In these mutant mice, the remaining mitochondria display a pathological morphology characterized by swelling and disruption of cristae, associated with an increased uncoupling of respiration (Andrews et al., 2010). These findings suggest that circadian genes can modulate the respiratory defects. Unsurprisingly, therefore, protein analysis in the Bmal1−/− mice identified a reduced expression of the α6(VI) chain. Interestingly, a severe depletion of the α6(VI) chain has recently been identified both in skeletal muscle and muscle cell cultures of ColVI patients, and this protein change has been proposed as a diagnostic marker of ColVI diseases (Tagliavini et al., 2014a). Indeed, timing plays a crucial role in muscle remodeling. For instance, Dadgar et al. (2014) demonstrated that asynchrony drives fibrosis and regeneration failure in the mdx model, finding a link between profound deregulation of global muscle transcription and the muscle damage seen in various dystrophies, including dystrophinopathies and several LGMD types. In addition, BMAL1 has recently been found to be involved in myogenic response to muscle damage, orchestrating the regeneration processes in adult skeletal muscle (Chatterjee et al., 2015).
Although these studies suggest that the circadian circuit plays a crucial role in muscle function and regeneration, the mechanism mediating these effects in skeletal-muscle-wasting diseases remains to be elucidated. In summary, our results demonstrate, for the first time, a deregulation of circadian genes in skeletal muscle of ColVI myopathy that might lead to characterizing the involvement of circadian clock circuit in many other muscle-wasting diseases, opening new therapeutical perspectives.
Putative models linking collagen VI myopathies with circadian rhythm deregulation
The new data suggest a pathogenic connection between circadian rhythm gene deregulation and collagen-VI-related myopathies. It appears that CLOCK gene deregulation could provide a link between autophagy, mitochondrial abnormalities and ColVI diseases. In light of the asynchrony recently described as the mechanism driving regeneration failure (Dadgar et al., 2014), this relationship could be expressed as one of the two following scenarios. First, that circadian rhythm deregulation is a downstream effect of ColVI deficiency and contributes to muscle damage pathogenesis (Fig. 6A). In this hypothesis, the crucial player interconnecting ColVI deficiency with circadian rhythm could be the inappropriate persistence of AKT1 activation. Owing to the reduced ColVI expression, this would lead not only to the downregulation of autophagy (Grumati and Bonaldo, 2012), but also to the activation of the HIF1A transcription factor and the consequent downregulation of GSK3-β kinase. Indeed, muscle regeneration is a finely regulated process which drives the transcription of myogenic factors through the orchestration of cell contact, growth factors and various hormones. Murine models lacking or constitutively expressing AKT genes show severe muscle deficiency or extensive hypertrophy, respectively (Wilson and Rotwein, 2007; Peng et al., 2003; Lai et al., 2004), implicating AKT genes (AKT1 and AKT2) in both the development and regeneration of muscle. Serra et al. (2007) hypothesize that this might arise through convergence of the AKT pathway and the p38 MAPK circuit at the chromatin level controlling the assembly of myogenic transcriptome and driving the regeneration process.
The second scenario is that circadian rhythm deregulation in ColVI myopathies is an independent event acting as a secondary modifier of disease severity (Fig. 6B). Our observation that the collagen VI α6 chain is downregulated in a mouse circadian model could suggest that circadian genes are genetic modifiers able to influence the primary ColVI defect. It is well recognized that the molecular clock in skeletal muscle can be modified by environmental factors including time of feeding or time of physical activity (Wolff and Esser, 2012; Hatori and Panda, 2015). Merging this concept with the known direct control that circadian genes have on MyoD1 (muscle development), PGC1 (mitochondria) and ATG14 (autophagy circuit) proteins provides the model through which circadian rhythm disruption could function independently to modify disease severity. It might therefore be enlightening to explore circadian genes in the many ColVI patients whose causative mutations are still unknown (recently mutations in COL12A1 gene have been linked to the UCMD and Bethlem myopathy phenotype) (Hicks et al., 2014).
In conclusion, our demonstration of circadian gene involvement in the pathogenesis of ColVI myopathy in both an animal model and human patients opens new avenues for both pathogenesis and treatment of hereditary muscle diseases. Timing and synchronization might have different effects on gene expression, depending on the specific gene mutation (COL6 or DMD). Different key regulators, such as AKT1 for COL6 myopathies and TGFB1 for dystrophin mutations (Dadgar et al., 2014), might be implicated, as these two muscle diseases have different pathogenic signatures, represented by defective autophagy and massive muscle necrosis, respectively. Finally, given that CLOCK is downregulated only in the severe UCMD phenotype, it should be further explored as a severity biomarker in patients with myopathies characterized by a severe muscle remodeling, with a view to applications in both disease stratification in patients and monitoring in clinical trials.
MATERIALS AND METHODS
Analysis of Col6a1−/− mice
For the focus of this work, we used standardized timing of muscle tissue collection from the Col6a1−/− mice and the Bmal1−/− mice. Mice in all colonies were maintained in normal 12-h-light–12-h-dark schedules with lights on at 7 am (Zeitgeber time – ZT0). Muscle tissues were collected between 10 am and 12 pm (ZT3–5) and frozen immediately in liquid N2. By following these strict time of collection procedures this allows comparisons across tissues and between animals. Both WT and Col6a1−/− mice belong to C57BL/6 strain regarding the fiber type composition of different muscles as previously described (Augusto et al., 2004).
Total RNA was isolated from three muscles (diaphragm, gastrocnemius and tibialis) of male mice. Mice were wild type (four) and knock out (four), at an age of 6 months. TRIzol extraction (Invitrogen), as specified by the manufacturer, was used. RNA from each muscle was quantified with Nanodrop (Thermo Scientific) spectrophotometer. Its integrity was evaluated with a 2100 Bioanalyzer (Agilent) prior to pooling equal amounts of RNA from muscles of the same kind.
Analysis of whole genome expression microarray data
We used the whole mouse genome gene expression kit including four arrays of 44,000 probes (Agilent Technologies, G4122F) to measure the transcriptional profile of muscle tissue samples of four mice. Sample labeling and hybridization were performed in triplicate according to protocols provided by Agilent (One-Color Microarray-Based Gene Expression Analysis version 5.0.1). The hybridized array was processed by Agilent scanner and Feature Extraction software (v9.5). The uniformness of array hybridization was verified through Agilent Quality Controls (Spike-in), consisting of a mixture of ten in-vitro-synthesized, polyadenylated transcripts derived from the Adenovirus E1A gene, premixed at concentrations spanning six logs and differing by one-log or half-log increments.
Microarray data analysis
Agilent's Feature Extraction Software (v9.5) was used for array image processing. Intensity data with local background adjustment by spatial detrending were log-transformed and quantile-normalized. Synonymous probes from the same gene were averaged. The software selected 28,052 probes for further analysis. To perform multi-way comparisons we utilized the limma package (Smyth et al., 2005). This package fits a separate linear model for each gene and uses moderated t-statistics to estimate differences between arrays. Empirical Bayes analysis is used to improve statistical power in small sample sizes. P-values were corrected for multiple testing using the Benjamini and Hochberg correction (Benjamini and Hochberg, 1995).
To ascribe biological changes in knockout versus WT muscles we used Gene Ontology term enrichment analysis. Gene ontology terms enriched in the differentially expressed genes were identified using Fisher exact test. We also performed Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) using Pathway Studio 9 from Elsevier.
Potential transcription regulators of the observed expression changes were identified using sub-network enrichment analysis (Sivachenko et al., 2007) with default parameters implemented in Pathway Studio. This algorithm identifies sub-networks consisting of a single expression regulator (for example, transcription factor) and its downstream genes showing concordant changes in expression between two conditions. Information about gene–gene regulations is extracted from the literature using natural language processing and stored in the Pathway Studio database. At the time of analysis, the database contained over 250,000 expression relationships between proteins. Relationships were extracted by natural language processing of more than 22 million PubMed abstracts and 880,000 full-text articles.
Microarray data analysis validation by fluidic cards
A total of 46 genes were selected for fluidic card analysis based on both differential expression between wild-type and knockout mice and literature analysis in Pathway Studio to find genes associated with the circadian rhythm pathway, muscle regeneration, autophagy and apoptosis. Selected genes were used to set up a custom TaqMan low density array (TLDA or fluidic card) according to manufacturer specifications in order to validate results obtained using the microarray platform. The major intrinsic protein of lens fiber (MIP) gene was selected as negative control for its absence of expression in the microarray data. Single fluidic card can analyze 48 genes (47 test and the 18S control housekeeping gene) simultaneously in eight samples and was run five times for each sample, it also included major autophagic effector protein Beclin 1 (Grumati et al., 2010).
150 ng for each RNA pool was retrotranscribed using High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems), according to the manufacturer's instructions and loaded in a single port of the TLDA run on ABI 7900HT System (Applied Biosystems) for 2 min at 50°C, 10 min at 95°C, followed by 40 cycles for 15 s at 97°C and 1 min at 60°C. Data analysis was performed following the ΔΔCT Method (Applied Biosystems User Bulletin #2) and using the 18S gene as calibrator and the WT samples as controls with respect to the knockout samples. Calculated mean 2−ΔΔCT values were log transformed and compared with normalized fold change microarray data.
Western blotting and RNA analysis in Bmal1−/− mice
Immunoblotting was used to evaluate the protein levels of LC3A (rabbit polyclonal, Thermo Scientific, PA1-16930; 1:1000), BECLIN1 (D40C5 rabbit mAb, Cell Signaling Technology, 3495; 1:1000), p62/SQSTM1 (Guinea pig polyclonal, Progen GP62-C; 1:1000), AMPK (rabbit polyclonal, Cell Signaling Technology, 2532; 1:1000), phosphorylated AMPK (rabbit polyclonal, Cell Signaling Technology, 2531; 1:1000), COL6A1 (COL6A1 H-200, Santa Cruz Biotechnology, sc-20649, 1:1000), COL6A3 and COL6A6 [rabbit polyclonal antisera against mouse a3(VI) and a6(VI) were kindly supplied by R. Wagener, Cologne, Germany]. Gastrocnemius muscles from 30–35-week-old male wild-type and Bmal1−/− (Andrews et al., 2010; Bunger et al., 2000) mice (both n=4) were homogenized using extraction buffer (2% SDS, 10% glycerol, 2% 2-mercaptoethanol, 50 mM Tris-HCl, pH 8.8). The homogenized samples were immediately incubated at 80°C for 5 min and protein concentration was measured using a Bio-Rad DC protein assay kit. The protein samples (50 µg/lane) were then subjected to 7.5% (for ColVI) or 10% SDS-PAGE. After electrophoresis, proteins were transferred to PVDF membrane (Millipore). The membrane was blocked in Odyssey blocking buffer (LI-COR) for 1 h and then incubated with primary antibody in 50% Odyssey blocking buffer (LI-COR) diluted with 1× PBS containing 0.1% Tween 20 for 1 h. The membrane was then washed and incubated with horseradish peroxidase (HRP)-conjugated secondary Ab (for all collagen VI proteins and p62) or Alexa-Fluor-680-conjugated secondary antibody diluted in in 50% Odyssey blocking buffer (LI-COR). The protein signal was detected either by incubation with ECL substrate and developing signals by exposing the membrane to CL-Xposure film (for HRP-conjugated secondary antibody), or by using the Odyssey infrared imager (LI-COR, for Alexa Fluor 680 antibody). The protein bands were analyzed by using ImageJ. The plots in the figures are average of four samples for both wild type and Bmal1−/−. For the Bmal1−/− studies, the animal procedures were conducted in accordance with institutional guidelines for the care and use of laboratory animals as approved by the AAALAC accredited, University of Kentucky Institutional Animal Care and Use Committees with the protocol number: 00728M2004.
Analysis of Bethlem myopathy and UCMD patients
The three ColVI patients have the following mutations: UCMD patient, heterozygous c.819_833del p.Pro274_Gly278del in COL6A1 gene; severe Bethlem myopathy patient: homozygous c.1393C>T p.Arg465* in COL6A3 gene; Bethlem myopathy patient, heterozygous c.428+1G>A p. Tyr122_Gly143del in COL6A1 gene. We obtained muscle biopsies from ColVI patients after informed consent and approval of the Ethics Committee of the University of Ferrara (no. 02/2009, 26 February 2009). The research involving human subjects was performed according to the Helsinki Declaration.
Muscle biopsy procedures were taken at a time according to the local clinic rules, ranging from 8–10 am (central Europe time) or 10–12 am (Greenwich time), depending on the sites. Fresh biopsies were frozen following routine procedures. Because of the similarity in time of collection, comparison of molecular clock components across patient samples can be made. The muscle fiber composition in these three muscles shows a predominance of fiber type I (Fig. S4). Generally, hereditary myopathies, including the collagen VI myopathies, are characterized by the ﬁber type disproportion with type I ﬁber hypotrophy and predominance (Schessl et al., 2008).
Total RNA was isolated from skeletal muscle biopsies of three patients with different phenotypes: UCMD, dominant and recessive Bethlem myopathy and from one healthy control biopsy using the RNeasy-kit (Qiagen, Chatsworth, CA) according to the manufacturer's instructions. We verified high-quality of isolated RNA and RNA integrity using an Agilent Technologies 2100 Bioanalyzer. The isolated RNA had an RNA Integrity Number (RIN) value of greater than 8.
1 μg of total RNA sample was used for poly(A)+ mRNA selection using streptavidin-coated magnetic beads, followed by thermal mRNA fragmentation. Fragmented mRNA was retro-transcribed using reverse transcriptase (Super-Script II) and random primers. The cDNA was converted into double-stranded cDNA, its ends were repaired using Klenow fragment, T4 polynucleotide kinase and T4 polymerase, and then cDNA was ligated to Illumina paired end adaptors. Size selection was performed using a 2% agarose gel, generating cDNA libraries ranging in size from 200 to 250 bp. Finally, the libraries were enriched using 15 cycles of PCR (30 s at 98°C; 15 cycles of 10 s at 98°C, 30 s at 65°C, and 30 s at 72°C; then 5 min at 72°C) and purified by the QIAquick PCR purification kit (Qiagen). The enriched libraries were diluted with elution buffer to a final concentration of 10 nM. Each library was run at a concentration of 7 pM on one Genome Analyzer (GAIIe) lane using paired-end 100-bp sequencing.
We performed mapping of reads to human genome using TopHat software. TopHat, the splice junction mapper for RNASeq reads, was run with default parameters. Supplied annotations from hg19 Refseq refGene were downloaded from the UCSC Genome Browser. The mapping results from TopHat were then passed to Cuffdiff, a tool which finds significant changes in transcript expression. In Cuffdiff quartile normalization was turned on, bias correction was set on, the minimum alignment count was set to 10, and all other parameters were set to default. Both TopHat and Cuffdiff were run on the free pubic Galaxy server (http://galaxyproject.org/). To identify consistently deregulated genes, differential probes were identified by the OK test in Cufflinks (Trapnell et al., 2010). Log ratios were calculated as the difference between patient and control.
Quantification of CLOCK and MAT2A transcripts by real-time PCR
To validate the deregulated expression of CLOCK and MAT2A genes, we enlarged our patient cohort. We analyzed 12 UCMD and nine Bethlem myopathy, three patients with intermediate phenotypes and a pool made of muscle controls. All ColVI patients have been clinically characterized with regard to ‘circadian involvement’, date and time of biopsy, if the patients were fasting or taking drugs, if the patients were ambulatory and the presence of ventilation at the time of biopsy (Table S3). Total RNA was isolated from muscle sections using the RNeasy Kit (QIAGEN, Chatsworth, CA) and reverse transcribed by using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems). In order to quantify the steady state level of CLOCK and MAT2A transcripts, commercially available TaqMan expression assays (Applied Biosystems) were used (CLOCK: Hs00231857_m1 NM_001267843 exon boundary 19–20; MAT2A: Hs01553622_g1 NM_005911 exon boundary 5–6) as for β-actin as housekeeping reference gene (ACTB Endogenous Control). Real-time PCR was performed in triplicate on the Applied Biosystems Prism 7900HT system, using 10 ng of cDNA and default parameters. Evaluation of transcripts levels was performed by the comparative CT method (ΔΔCT method). Statistical analysis was performed with a Student's t-test.
Western blotting analysis in Col6a1−/− mice and ColVI patients
For immunoblotting analysis, frozen muscle tissues from both patients and mice were ground in a dry-ice-cooled mortar and stored at −80°C. An aliquot of each frozen muscle was suspended in lysis buffer (urea 7 M, thiourea 2 M, CHAPS 4%, Tris-HCl 30 mM, and PMSF 1 mM) and solubilized by sonication on ice. Samples were pooled and protein concentrations were determined using PlusOne 2D-Quant kit (GE Healthcare, Milan, Italy). Total protein extracts (50 µg) from pooled Col6a1−/− (5 mice) and WT (5 mice) animal models, and from pooled UCMD (4 patients), Bethlem myopathy (8 patients), 30 unrelated patients with Duchenne (DMD) or Becker (BMD) muscular dystrophies, and healthy control (5 subjects) muscles, were loaded in triplicate and resolved on 14–8% gradient polyacrylamide gels. Blots were incubated with rabbit polyclonal primary antibodies as follows: anti-MAT2A (Novus Biologicals, NB110-94158; 1:1000), anti-CLOCK (Santa Cruz Biotechnology, sc-25361; 1:500). After washing, the membranes were incubated with anti-rabbit-IgG (GE Healthcare) secondary antibody conjugated to horseradish peroxidase (1:5000). The signals were visualized by chemiluminescence using the ECL Prime (GE Healthcare) detection kit and the Image Quant LAS 4000 (GE Healthcare) analysis system. Statistical analysis was performed applying the ANOVA and Tukey tests (n=3, P<0.05).
Manipulations of Col6a1−/− animals, required for the experiments described, conform to the rules of the Animal Experimental Commission of the ‘Complesso Interdipartimentale A. Vallisneri, University of Padua (Italy)’, where the work has been carried out. These regulations represent the application of the D.L. 27/1/92 n. 116 of the Italian Government, which meets Directive 86/609/EEC of the EU Council concerning the protection of animals used for experiments and scientific purposes. The presented procedures have been approved by the Animal Experimental Commission with the certificate no. 22675.
Special thanks are due to Mikhail Pyatnitskiy from Ariadne-Dx for his data analysis. We are grateful to Raimund Wagener (University of Cologne, Germany) for a3(VI) and a6(VI) for kindly providing collagen antibodies. The Eurobiobank is also acknowledged.
A. Ferlini designed the rationale of the research and supervised the work; C. Scotton, M. Bovolenta, M. S. Falzarano, E. Martoni, C. Passarelli, A. Armaroli, H. Osman and R. Selvatici performed the experimental procedures related to expression profile analysis and data validation; T. Castrignanò and G. Pesole performed the bioinformatics analysis of RNAseq. E. Schwartz, E. Kotelnikova, A. Yuryev and M. Lebowitz performed the interactome analysis and statistical analysis; C. Scotton, E. Schwartz, K. Esser and A. Ferlini performed the data interpretation; C. Rodolico, S. Messina, E. Pegoraro, A. D'Amico, E. Bertini, F. Gualandi, M. Neri, , P. Boffi, M. A. Maioli, H. Lochmüller, V. Straub, K. Bushby, P. Sabatelli, L. Merlini, R. Foley, S. Cirak, I. Zaharieva and F. Muntoni performed the clinical assessment and provided patient samples; P. Braghetta, P. Bonaldo and P. Bernardi provided and characterized the Col6a1−/− mice; D. Capitanio and C. Gelfi performed the proteomic analysis on Col6a1−/− mice and patients. X. Zhang, B. Hodge and K. A. Esser performed the proteomic analysis on Bmal1−/− mice; E. Schwartz, P. Bonaldo, P. Bernardi, F. Muntoni, H. Lochmüller and K. A. Esser critically revised the paper; C. Scotton and A. Ferlini wrote the paper. The authors M. Bovolenta and E. Schwartz contributed equally to this work.
This study was funded by the BIO-NMD European Union Seventh Framework Programme [project 241665 to A.F., E.S., F.M., H.L., C.G. and P. Bonaldo); the Fondazione Telethon [grant numbers GGP08017 to P. Bernardi, P. Bonaldo., C.G. and A.F.; GUP11007 to L.M. and P.S.]; the RD-Connect EU project 305444 [to H.L.]; and the National Institutes of Health AR066082 to K.A.E.]. The financial support of the Muscular Dystrophy Campaign Centre Grant and of the MRC Neuromuscular Centre at UCL; of the Biomedical Research Centre at Great Ormond Street Hospital and of the Great Ormond Street Hospital Children's Charity (in support of F.M.) are also gratefully acknowledged. Deposited in PMC for release after 12 months.
The authors declare no competing or financial interests.