Cellular metabolic processes are highly enriched in the circadian transcriptome of skeletal muscle

To identify circadian gene expression in skeletal muscle, we used a publicly available, high-resolution, circadian time-course microarray dataset from gastrocnemius muscles of male C57BL/6 mice [45,46]. These mice were housed in constant darkness, and food was provided ad libitum to eliminate the influence of external environmental cues. We chose this dataset because it has double the sampling frequency of previously published circadian muscle transcriptomes, and this allows for greater precision for circadian analysis [46,54]. Using the JTK_CYCLE statistical algorithm [47] for the reliable detection of oscillating transcripts with a 24-h periodicity, we identified 1,628 circadian mRNAs (adjusted P < 0.05). An unbiased Gene Ontology enrichment analysis of these circadian genes revealed a significant overrepresentation of cellular metabolic processes, with approximately 1,004 (62%) genes directly involved in skeletal muscle metabolic processes as well as the regulation of metabolism (Figure 1).

Figure 1 Gene ontology analysis of the skeletal muscle circadian transcriptome. Top 15 enriched GO processes listed from left to right in order of significance. Full size image

An additional benefit of using the JTK_CYCLE algorithm is its ability to determine the acrophase, or time of peak expression, of each circadian probeset. Identifying the acrophase of genes that have common ontologies may help to predict the potential timing of cellular and physiological processes. Herein, we report the acrophase according to their respective circadian times (CT), which is standardized to the free-running period of the mice under constant conditions. For the array studies, the mice were in DD for 30 h so CT 0 denotes the start of the inactive period, while CT 12 denotes the start of the active period. To identify the timing of gene expression and its relationship to metabolic processes in skeletal muscle, we annotated a subset of circadian genes by their known functions, timing of peak expression, and involvement in key metabolic pathways. We focused our analysis on metabolic functions that involve substrate (carbohydrate and lipid) utilization as well as storage and biosynthetic processes.

Lipid metabolism: genes involved in fatty-acid uptake and β-oxidation peak in the mid-inactive/light phase

Skeletal muscle expresses specialized membrane transporters to facilitate the transport of lipids into the cell [55-57]. Two lipid transport genes that encode for fatty-acid binding proteins, Fabp4 (CT 24.0) and Fabp3 (heart/muscle isoform, CT 6.0), are expressed in a circadian manner with the highest mRNA expression in the early- and mid-inactive periods, respectively. Acrophase of circadian genes involved in lipid metabolism are illustrated in Figure 2. Normalized expression traces for each gene are located in Additional files 1, 2, and 3. Previous studies have demonstrated oscillations in plasma fatty acid concentrations in mice with peak levels occurring during the inactive/light period [58-60]. Further functional analysis is required to validate the predition that the rate of fatty-acid uptake in skeletal muscle peaks during the mid-late inactive period. Upon uptake into the cell, fatty acids can be stored as triglycerides or be converted to acetyl-CoA through β-oxidation [61]. Slc25a20 encodes for an acyl-carnitine translocase that transfers fatty acids into the inner-mitochondrial matrix and reaches peak expression in the middle of the inactive period (CT 7.5) [62]. We identified multiple genes that encode for β-oxidation enzymes to be circadian and also reach peak expression around the mid-inactive phase. These include the enoyl CoA hydratase Ech1 (CT 7.0), the tri-functional enzyme subunits Hadha (CT 8.0) and Hadhb (CT 8.0), and the acetyl-CoA acyltransferase Acaa2 (CT 9.0). Malonyl-CoA, an intermediate formed during de novo fatty acid synthesis, is a potent inhibitor of β-oxidation. The striated muscle enriched gene Mlycd (CT 7.5) encodes for the malonyl-CoA decarboxylase that promotes β-oxidation by reducing cytosolic concentrations of malonyl-CoA and reaches peak expression during the mid-inactive period similar to that of the circadian β-oxidation genes. These observations suggest that rates of β-oxidation are modulated over time of day and potentially through the endogenous molecular clock in skeletal muscle [10,63,64].

Figure 2 Schematic acrophase diagram of circadian genes involved in lipid metabolic processes. The relative location of the circadian genes (italicized) in respect to the x-axis indicates acrophase or time of peak expression calculated by the JTK_CYCLE algorithm. Location of substrates and pathways does not represent peak substrate concentrations and/or rates of individual pathways as these were not measured in our analysis. White/grey shading is representative of the inactive and active phases, respectively. Full size image

Nuclear receptors are known to be potent transcriptional regulators of metabolism as they sense changes in environmental conditions and induce appropriate changes in the expression of metabolic genes [65-69]. The nuclear receptor Estrogen-related receptor alpha (Esrra, CT 7.5) and the nuclear co-activator PPARγ coactivartor-1 beta (Ppargc1b, CT 7.0) are both circadian genes in skeletal muscle with peak expression occurring at the mid-inactive phase. These factors have been shown to promote mitochondrial biogenesis, fatty-acid uptake (targets Fabp3), and β-oxidation [70,71]. The nuclear co-repressor Nrip1, also known as Rip140, is a potent negative regulator of skeletal muscle oxidative metabolism and has been shown to suppress expression of the fatty-acid transporter, Fabp3, in skeletal muscle [72-74]. NRIP1 suppresses gene expression by binding nuclear receptors (including PPARs and estrogen-related receptors) and recruiting histone deacetylases [75]. Interestingly, peak expression of Nrip1 occurs during the beginning of the active period (CT 13.0) and may therefore act as a molecular brake to oxidative metabolism as the muscle transitions from lipid to carbohydrate utilization during the early active phase.

Lipid metabolism: lipogenic genes reach peak expression at the end of the active/dark phase

The lipogenic genes Acly (CT23.0), Acaca (CT 23.0), and Fasn (CT 22.5) involved in de novo fatty-acid synthesis, or the conversion of excess carbohydrates into fatty acids, reach peak expression at the end of the active phase (Figure 2) [61,76]. Scd1 (CT 24.0) encodes the enzyme that catalyzes the rate-limiting reaction of monounsaturated fatty-acid formation to promote lipid bilayer fluidity and lipogenesis [77,78]. The genes Srebf1 (CT 24.5), Srebf2 (CT 24.0), and Mlxip (CT 23.5) encode transcription factors that target carbohydrate response elements within lipogenic gene promoter regions (Acly, Acaca, and Fasn) and are also circadian with peak expression at the end of the active phase [79,80]. Consistent with our results, Srebf1 oscillations have been reported in the liver and genome-wide binding studies have shown a circadian recruitment pattern of SREBF1 to the promoters of lipogenic genes with maximal binding during the active (fed) stage [81-84].

The gene Pnpla3 (CT 21.0), also known as adiponutrin, promotes lipogenesis by converting LPA to phosphatidic acid (PA) [85]. The gene Lpin1 (CT 24.0) which encodes for the lipin-1 enzyme is responsible for converting phosphatidic acid (PA) to diacylglycerol (DAG), the upstream metabolite required in phospholipid biosynthesis [86,87]. The highly regulated, committing step in triacylglycerol (TG) synthesis, addition of a fatty-acyl-CoA to DAG, is performed by the product encoded by Dgat1 (CT 24.5), which is also expressed in a circadian manner [88]. Once a TG molecule is formed, it can be elongated by enzymes encoded by Acsl5 (CT 23.0) or Elovl5 (CT 22.5) [89,90]. The observation that circadian lipogenic genes reach peak expression levels around the end of the active phase suggests that skeletal muscle promotes storage of excess energy at the end of the active/absorptive phase.

Carbohydrate metabolism: genes involved in carbohydrate catabolism peak in the early active/dark phase

Glycolysis, the breakdown of glucose to form pyruvate, is primarily regulated at two enzymatic reactions catalyzed by the hexokinase and phosphofructokinase enzymes [91]. We observe that the hexokinase-2 (Hk2) gene is circadian with peak expression occurring at the beginning of the active phase (CT 12.0). Acrophase of circadian genes involved in carbohydrate metabolism are illustrated in Figure 3. Normalized expression traces for each gene are located in Additional files 1, 2, and 3. Hk2 is responsible for the first step in glycolysis by phosphorylating glucose to make glucose-6-phosphate, thereby trapping glucose within the cell [92]. The rate-limiting step of glycolysis involves the catalysis of fructose-6-phosphate to the highly unstable fructose-1,6-bisphosphate by the enzyme phosphofructokinase-1 (PFKM) [93,94]. A potent allosteric activator of PFKM is fructose-2,6-bisphosphate, which is the product of the other phosphofructokinase isozyme phosphofructokinase-2 (PFK2) [95]. Three genes (Pfkfb-1,3,4) that encode phosphofructokinase-2 subunits are circadian with peak expression occurring during the mid- and late-inactive phases (CT 10.0, CT 4.5, and CT 12.0, respectively).

Figure 3 Schematic acrophase diagram of circadian genes involved in carbohydrate metabolic processes. The relative location of the circadian genes (italicized) in respect to the x-axis indicates acrophase or time of peak expression calculated by the JTK_CYCLE algorithm. Location of substrates and pathways does not represent peak substrate concentrations and/or rates of individual pathways as these were not measured in our analysis. White/grey shading is representative of the inactive and active phases, respectively. Full size image

Glycolytic flux through the Kreb’s cycle is controlled by pyruvate dehydrogenase complex (PDH) [96,97]. PDH decarboxylates pyruvate to form acetyl-CoA, which is a substrate for the Kreb’s cycle. The activity of PDH is regulated at the posttranslational level. Phosphorylation by kinases (PDKs) inhibits PDH activity, while dephosphorylation by phosphatases (PDPs) activates the complex [98,99]. The Pdk4 gene, which encodes for a PDH kinase that inhibits PDH, reaches maximal expression at the mid-inactive phase (CT 6.0). This expression pattern is similar to that of the β-oxidation genes and suggests that skeletal muscle substrate preference is pushed toward utilization of lipids over carbohydrates during the mid- to late-inactive phase. Conversely, the PDH phosphatase gene, Pdp1, peaks at the beginning of the active phase (CT 10.0) in a similar temporal fashion compared to the glycolytic enzymes described above. This temporal regulation of Pdp1 may therefore help increase glycolytic flux during the active phase. Dyar et al. observed similar expression patterns of Pdk4 and Pdp1 in skeletal muscle and were first to report a shift to carbohydrate utilization at the beginning of the active phase [43].

Adrb2 encodes for the β2-adrenergic receptor (β 2 AR) involved in the fight-or-flight response in peripheral tissues [100,101]. Agonist (that is, catecholamine) binding is well established to evoke a cell-signaling cascade that promotes glucose uptake, glycogenolysis, and lipolysis to provide a readily available source of energy for skeletal muscle [102-104]. Adrb2 is expressed in a similar pattern to that of the glycolytic activating genes as it peaks at the beginning of the active phase. Interestingly, the expression of Adrb2 coincides with that of oscillating epinephrine concentrations in mammals, which has previously been identified as peaking at the beginning of the active phase in mouse models [105]. The G-protein receptor kinase, encoded by Adrbk1, phosphorylates the β 2 AR, thereby rendering it susceptible to receptor-mediated endocytosis via β-arrestin proteins encoded by Arrdc3 and Arrb1 [106-108]. Adrbk1, Arrdc3, and Arrb1 are all expressed in a circadian manner and antiphasic to the expression of Adrb2. These observations suggest there is a time of day difference in adrenergic signaling and that sensitivity to epinephrine may be highest in skeletal muscle during the active period while being desensitized prior to the inactive period.

Carbohydrate metabolism: genes involved in carbohydrate storage peak at the mid-active/dark phase

Excess carbohydrates are stored as glycogen in skeletal muscle which accounts for approximately 70 to 80% of whole body stores [109]. Unlike the liver, skeletal muscle glycogen content is not responsible for maintaining blood glucose concentrations but serves as a rapidly accessible energy depot for active contractions [110]. Glycogenesis is regulated by both glucose-6P concentrations and the enzymatic activity of glycogen synthase [111,112]. The gene Ppp1r3c (CT 20.0) reaches peak expression around the mid-inactive phase and encodes a regulatory subunit of the protein phosphatase-1 (PP-1) responsible for activating glycogen synthase while also inhibiting glycogen breakdown (Figure 3) [113]. Enzymatic activity of PP-1, and subsequent activation of glycogen synthase, is regulated downstream of the insulin signaling pathway [114].

Insulin promotes an anabolic signaling cascade that works in opposition to that of adrenergic signaling to drive glycogen and lipid storage. Previous reports have identified a ‘counter-regulatory’ role of the insulin receptor to selectively inhibit β 2 AR signaling through phosphorylation and subsequent internalization of the receptor [101,115]. Interestingly, the genes that encode the insulin receptor substrate-1, Irs1 (CT 22.0), and its downstream PI3-kinase target, Pik3r1 (CT 19.0), are both circadian with peak expression occurring at the late-active phase while the genes involved in suppressing PI3-kinase, Pik3ip1 (CT 8.0), and the insulin-receptor substrate-1, Fbxo40 (CT 5.0), reach peak expression during the inactive phase [116,117]. These data suggest that the molecular clock may act to prime skeletal muscle to store excess glucose during the mid- to late-active phase. This prediction is further supported by previous studies that report skeletal muscle glycogen content as having a diurnal rhythm with the highest levels occurring during the mid-active phase [118-120]. Skeletal muscle glucose uptake is primarily controlled via the presence/absence of the glucose transporter GLUT4/Slc2a4 in the plasma membrane (sarcolemma) and transverse tubules. A t-SNARE syntaxin-4 interacting protein, encoded by Stxbp4, has previously been shown to repress GLUT4 insertion into the plasma membrane in the absence of insulin signaling [121-123]. The gene Tbc1d1 encodes for Rab-GTPase that represses GLUT4 translocation in the absence of insulin- or contraction-induced signaling cascades [124-126]. Interestingly, Tbc1d1 and Stxbp4 are both expressed in a circadian manner and reach peak expression in the middle of the active phase (CT 19.0). Previous reports have identified Tbc1d1 as a circadian gene in skeletal muscle and other tissues [43,127]. Together, these gene products may play a role in reducing glucose uptake at the end of the active phase by repressing GLUT4 translocation and/or insertion into the plasma membrane. This temporal separation of anabolic and catabolic signaling processes in skeletal muscle may be vital for maintaining a tight regulation of serum glucose levels, and disruption of which may contribute to the metabolic phenotypes often reported in clock-mutant mice models.

Generation of an inducible skeletal muscle-specific mouse model of Bmal1 inactivation

Use of the high-resolution microarray data set allowed for the identification of mRNAs expressed in a circadian pattern, but this could be due to the intrinsic molecular clock or could be a response to external behavioral (feeding/activity) or neural/humoral cues [24,128,129]. To determine the role of the intrinsic skeletal muscle molecular clock in the temporal regulation of metabolic gene expression, we generated an inducible mouse model to inactivate Bmal1 specifically in adult skeletal muscles. Upon treatment with tamoxifen in 12-week-old adult mice, we detect recombination of exon-8 (that is, DNA binding region) of the Bmal1 gene specifically in skeletal muscle (Figure 4A), confirming the tissue specificity of the mouse model. We waited until 12 weeks of age to limit possible developmental effects as BMAL1 has been shown to promote myogenesis [20,130]. As seen in Figure 4A, recombination was not detected in the skeletal muscle or non-muscle tissues of vehicle-treated mice (iMS-Bmal1 +/+). Western blot analysis confirmed the depletion of BMAL1 protein in the skeletal muscle of the iMS-Bmal1 −/− mice with no effect on the liver (Figure 4B). Tamoxifen-induced loss of Bmal1 in adult skeletal muscle resulted in significant and expected gene expression changes of genes involved in the core clock mechanism. In particular, genes directly activated by the BMAL1/CLOCK heterodimer, such as Rev-erbα and Dbp, are markedly downregulated in iMS-Bmal1 −/− but not in the iMS-Bmal1 +/+ samples (Figure 4C). Collectively, these results demonstrate the effective loss of BMAL1 protein and disruption of core-clock gene expression in the iMS-Bmal1 −/− muscle tissue.

Figure 4 Characterization of iMS-Bmal1 −/− mice. Recombination assay (A) of genomic DNA isolated from muscle and non-muscle tissues from tamoxifen-treated (iMS-Bmal1 −/−) and vehicle-treated (iMS-Bmal1 +/+) mice at 17 to 18 weeks of age (5 weeks postinjection). Recombination of the Bmal1 gene (exon 8) yields a 572-bp PCR product. The non-recombined allele is detected at 431 bp. Western blot (B) analysis of BMAL1 expression in iMS-Bmal1 −/− and iMS-Bmal1 +/+ liver and gastrocnemius samples. Note that the original blot containing both muscle and liver samples was cut, and brightness/contrast was altered to enhance the visibility of Bmal1 in the muscle samples. (C) Real-time PCR results of time-course expression values for Bmal1, Rev-erbα, and Dbp in the iMS-Bmal1 +/+ (black) and iMS-Bmal1 −/− (red). Representative wheel running rhythms (D) of iMS-Bmal1 −/− and iMS-Bmal1 +/+ mice. White and black bars (top) indicate light and dark phases. 12 L/12D represents the 12-h light/12-h dark cycle. 12D/12D represents constant darkness conditions. Tick marks indicate wheel running activity. Representative chi-squared periodograms (E) of iMS-Bmal1 −/− and iMS-Bmal1 +/+ mice indicating approximate 24-h period lengths in both mice. Full size image

iMS-Bmal1−/− display normal circadian activity rhythms

We used voluntary wheel running to assess circadian behavior in the iMS-Bmal1 mice 22 to 29 weeks posttreatment. We did not detect any significant differences in entrainment to light under 12-h light/12-h dark conditions between iMS-Bmal1 +/+ and iMS-Bmal1 −/−, and analysis of activity rhythms under constant darkness did not reveal any changes in circadian behavior (Figure 4D,E). Clock-lab analysis indicates that both iMS-Bmal1 +/+ and iMS-Bmal1 −/− exhibit approximate 24-h period lengths (23.85 ± 0.083 and 23.77 ± 0.138 h, respectively) with no differences in amplitude, the relative strength of the rhythm. These data are consistent with other studies and confirm that inactivation of BMAL1 in skeletal muscle does not directly alter circadian activity patterns [43,131]. Therefore, gene expression changes observed in this model are more likely to be downstream of the endogenous molecular clock mechanism in skeletal muscle.

Expression of key circadian metabolic genes are significantly altered in iMS-Bmal1 −/− skeletal muscle

Gene expression analysis of iMS-Bmal1 +/+ and iMS-Bmal1 −/− muscle tissue reveals that the intrinsic molecular clock, even in constant conditions, plays a role in temporally regulating carbohydrate and lipid metabolism. We performed our transcriptome analysis at 5 weeks postrecombination to identify early gene expression changes caused by the loss of the clock mechanism in skeletal muscle. Analyzing gene expression at this time point also limits potential off-target effects of tamoxifen treatment by allowing for a sufficient wash-out period. We found that the circadian genes involved in carbohydrate metabolism were most affected by loss of Bmal1. The expression of the glycolytic enzymes, Pfkfb1, Pfkfb3, and Hk2 as well as the PDH phosphatase, Pdp1 were all significantly downregulated in the gastrocnemius (Figure 5A). In addition, expression of the adrenergic receptor, Adrb2, was also significantly decreased. These genes are convincing clock-controlled candidates in skeletal muscle as they have circadian expression patterns similar to that of known clock-controlled genes (peak expression during inactive to active phase transition), and their loss of expression following Bmal1 inactivation is indicative of direct transcriptional regulation by the clock. By targeting these genes, the molecular clock mechanism can precisely regulate the timing of carbohydrate utilization to occur during the active phase. The observation that circadian genes involved in glucose utilization are diminished in our model is in agreement with the muscle-specific Bmal1 knockout model generated by Dyar et al. in which they report significant decreases in glucose oxidation and insulin stimulated glucose uptake in their muscle tissues [43].

Figure 5 Differentially expressed circadian, metabolic genes in iMS-Bmal1−/− skeletal muscle. Average expression changes of the circadian carbohydrate (A) and lipid (B) genes in iMS-Bmal1−/− gastrocnemius averaged over circadian times 18, 22, 26, 30, 34, and 38. Tibialis anterior and soleus gene expression changes (Dyar et al.) averaged over circadian times 0, 4, 8, 12, 16, and 20. The red line denotes control (iMS-Bmal1 +/+) gene expression values. *P ≤ 0.05; **P ≤ 0.01; ***P ≤ 0.001. Full size image

Lipid metabolic processes appear to be elevated as the nuclear co-repressor, Nrip1, involved in repressing β-oxidation was significantly decreased with loss of Bmal1 (approximately 21% decrease, Student’s t test P value = 0.019). Previous studies have shown that knockout of Nrip1 results in an increase in succinate dehydrogenase staining of gastrocnemius muscle consistent with a shift to slow oxidative fiber types [72]. Interestingly, the fatty-acid transporter, Fabp3, and the β-oxidation genes, Hadha and Hadhb, were significantly elevated in the iMS-Bmal1 −/− gastrocnemius tissues (Figure 5B). Two circadian genes involved in triacylglycerol elongation, Pnpla3 and Elovl5, were also increased in the iMS-Bmal1 −/−. Altogether, we report significant expression changes in circadian genes that are key regulators of metabolism in skeletal muscle. We think that the gene changes observed in iMS-Bmal1 −/− are either directly or indirectly regulated downstream of BMAL1/molecular clock in skeletal muscle and not due to changes in external cues as circadian activity patterns in iMS-Bmal1 −/− are indistinguishable from vehicle-treated controls. The observation that circadian genes involved in carbohydrate and lipid metabolism are disrupted in iMS-Bmal1 −/− highlights a fundamental importance of the intrinsic molecular clock in temporal regulation of substrate utilization and storage in skeletal muscle in the absence of external cues.

iMS-Bmal1−/− gene expression changes reveal a fast to slow fiber-type shift

Skeletal muscle comprises different fiber types that are differentiated based on contractile function as well as predominant substrate utilization [132-135]. For example, fast-type skeletal muscles (type IIX/IIB) primarily rely on ATP generated from anaerobic metabolism (glycolysis/lactic-acid fermentation) to provide quick energy sources required for short bursts of activity, while slow-type skeletal muscles and fast-type IIA muscles rely on oxidative metabolism to promote a more sustained and less fatigable bout of contractions. We analyzed changes in gene expression related to fiber type following Bmal1 ablation in adult skeletal muscle and included both circadian and non-circadian transcripts. We identified a selective increase in slow-type sarcomeric genes in the gastrocnemius muscles with a limited effect on fast-type sarcomeric genes (Figure 6A,B). We chose the list of ‘slow’ and ‘fast’ sarcomeric genes, because these have been shown to be significantly enriched in either slow-soleus or fast-EDL myofiber preparations [136]. Additionally, calcium handling genes and nuclear receptors common in slow-fiber muscles (for example, Casq2, Atp2a2, Ankrd2, Csrp3.) were significantly increased in iMS-Bmal1 −/− (Table 1). Similar to the changes observed for the circadian metabolic genes, we see that non-circadian metabolic genes involved in carbohydrate metabolism are significantly decreased, while genes involved in lipid metabolism are increased (Tables 2 and 3). This switch from a fast to a slow fiber type mRNA profile is in agreement with the observed metabolic changes as slow fiber type muscles rely more heavily on oxidative metabolism compared to fast-type skeletal muscle.

Figure 6 Increase in slow type sarcomeric genes in iMS-Bmal1 −/−. Average gene expression changes of slow (A) and fast (B) type sarcomeric genes in iMS-Bmal1 −/− compared to control values (red line). *P ≤ 0.05; **P ≤ 0.01; ***P ≤ 0.001. Full size image

Table 1 Fiber-type specific gene expression changes in iMS- Bmal1 −/− Full size table

Table 2 Metabolic genes upregulated in iMS- Bmal1 −/− Full size table