Skeletal muscle BMAL1 is necessary for transcriptional adaptation of local and peripheral tissues in response to endurance exercise training

Abstract

Objective: In this investigation, we addressed the contribution of the core circadian clock factor, BMAL1, in skeletal muscle to both acute transcriptional responses to exercise and transcriptional remodeling in response to exercise training. Additionally, we adopted a systems biology approach to investigate how loss of skeletal muscle BMAL1 altered peripheral tissue homeostasis as well as exercise training adaptations in iWAT, liver, heart, and lung of male mice.

Methods: Combining inducible skeletal muscle specific BMAL1 knockout mice, physiological testing and standardized exercise protocols, we performed a multi-omic analysis (transcriptomics, chromatin accessibility and metabolomics) to explore loss of muscle BMAL1 on muscle and peripheral tissue responses to exercise.

Results: Muscle-specific BMAL1 knockout mice demonstrated a blunted transcriptional response to acute exercise, characterized by the lack of upregulation of well-established exercise responsive transcription factors including Nr4a3 and Ppargc1a. Six weeks of exercise training in muscle-specific BMAL1 knockout mice induced significantly greater and divergent transcriptomic and metabolomic changes in muscle. Surprisingly, liver, lung, inguinal white adipose and heart showed divergent exercise training transcriptomes with less than 5% of 'exercise-training' responsive genes shared for each tissue between genotypes.

Conclusions: Our investigation has uncovered the critical role that BMAL1 plays in skeletal muscle as a key regulator of gene expression programs for both acute exercise and training adaptations. In addition, our work has uncovered the significant impact that altered exercise response in muscle and its likely impact on the system plays in the peripheral tissue adaptations to exercise training. Our work also demonstrates that if the muscle adaptations diverge to a more maladaptive state this is linked to increased gene expression signatures of inflammation across many tissues. Understanding the molecular targets and pathways contributing to health vs. maladaptive exercise adaptations will be critical for the next stage of therapeutic design for exercise mimetics.

Keywords: Circadian biology; Exercise; Inflammation; Metabolism; Signal transduction; Transcription.

Figure 1

Loss of skeletal muscle Bmal1 alters the transcriptional response to acute exercise. (A) Experiment 1. Vehicle and iMSBmal1KO mice completed a 30-minute, work-matched 15° inclined treadmill bout and gastrocnemius muscles were collected 1-hour post. (B) Exercise bout work was matched between genotypes. (C) RNA-sequencing was used to assess the number of DEGs 1-hour post exercise in vehicle and iMSBmal1KO mice. (D) Venn overlap of vehicle and iMSBmal1KO mice exercise responses. (E) Pathway analysis of shared exercise response genes. (F) Z-scores of exercise responsive clock factors and common exercise response genes shown between genotypes. Associated pathway databases are indicated by (R) Reactome, (K) KEGG, and (BP) Biological Processes. Pathway analysis of vehicle exercise responsive genes (G) and iMSBmal1KO exercise responsive genes (H).

Figure 2

Transcriptome and chromatin accessibility analysis of iMSBmal1KO mice. Vehicle treated and iMSBmal1KO mice underwent baseline physiological phenotyping including RNA and ATAC-sequencing analysis at ZT13 on gastrocnemius tissue. (B) Number of DEGs genes with loss of skeletal muscle Bmal1 (q = 0.05). Numbers in bars indicate up or downregulation. (C) Pathway analysis of differentially expressed RNAs split by direction of change. Numbers in bars indicate the number of genes present within the pathway. (D) Number of differentially accessible regions as assessed by ATAC-seq. Numbers in bars indicate the number of DNA regions/peaks becoming more or less accessible with loss of Bmal1. (E) Pathway analysis of differentially accessible chromatin regions split by direction of change. Numbers in bars indicate the number of genes present within the pathway. (F) Upset plot demonstrating shared features shared across acute exercise responses (DEGs), genes/chromatin regions affected by loss of muscle BMAL1 and (BMAL1:CLOCK ChIP-seq peaks, Gabriel et al. 2021 Sci Adv). (G) Pathway analysis of ‘BMAL1 regulated’ exercise responsive features. (H) Example BMAL1 ChIP-seq plots illustrating BMAL1 binding in promoter regions of Nr4a3 and Irs2.

Figure 3

Graded exercise performance and body composition analysis of vehicle and iMSBmal1KO mice. The experimental design is illustrated in (A) which shows the timecourse and parameters of exercise training, physiological testing and tissue collection. At ∼4 months of age, iMSBmal1^fl/fl were treated with either a vehicle, or with tamoxifen for 5 days (2 mg/day) to induce Cre-recombination and loss of Bmal1 specifically in skeletal muscle. Both vehicle and iMSBmal1KO mice were randomly assigned to 6 weeks of exercise training performed at ZT13, or control groups. Graded exercise tests were performed prior to, after 2, 4 and 6-weeks of training with the speed adjusted accordingly to make training progressive. MRIs were performed before, after 3 and 6 weeks of training. Mice from all groups were euthanized, and tissue was collected at ZT13, 47 h after their final graded exercise test. (B) Illustrates the weekly work done (AU) during exercise training. (C) Measurements of body mass were taken before, after 3 and 6 weeks of training and the subsequent percentage change calculated for each individual. (D) Total fat mass (g) measured before, after 3 and 6 weeks of training. (E) Percent changes in fat mass (g) were calculated from the beginning of the training period to the middle and end of the training period. (F) Total lean mass (g) measured before, after 3 and 6 weeks of training. (G) Percent changes in lean mass (g) were calculated from the beginning of the training period to the middle and end of the training period.

Figure 4

iMSBmal1KO gastrocnemius muscles present a larger and divergent transcriptional adaptation to 6-weeks of exercise training. (A) Volcano plot illustrating DEGs between control and iMSBmal1 WT exercise trained mice. A total of 142 (75 downregulated, 67 upregulated) genes were differentially expressed using an FDR cut off of <0.05. (B) HOMER Motif analysis performed on up and downregulated DEGs identifies a number of transcription factors regulating the transcriptional response to normal exercise adaptation. (C) KEGG Pathway analysis of genes associated with the vehicle treated transcriptional response to exercise training. (D) Volcano plot illustrating DEGs between iMSBmal1KO control and iMSBmal1KO exercise trained mice. A total of 780 (283 downregulated, 497 upregulated) genes were differentially expressed using an FDR cut off of <0.05. (E) HOMER Motif analysis performed on up and downregulated DEGs identifies a number of transcription factors regulating the transcriptional response to exercise adaptation in the absence of skeletal muscle Bmal1. (F) KEGG Pathway analysis of genes associated with the iMSBmal1KO transcriptional remodeling to exercise training. (G) Venn diagram overlapping of DEG responses to vehicle and iMSBmal1KO exercise training at P < 0.01. Minimal overlap in the transcriptional exercise responses is present between vehicle treated and iMSBmal1KO mice using a lenient significance cut off (63 out of 2107, 2.99% overlap). (H) Pie chart of the 63 common exercise training responsive genes in both genotypes reveals common and oppositely regulated gene reprogramming. (I) Analysis of DEGs affected by loss of BMAL1 and exercise training responsive DEGs in the iMSBma1KO mice, reveals 94 common genes. Pearson's correlation of the 94 common DEGs (Log² fold changes) reveals that genes affected with Bmal1KO are not rescued with exercise training. (J) Heatmap presenting Log² Z-scores for core clock factors/clock output genes and genes associated with glucose handling and glycogen utilization. Exercise training does not rescue the downregulation of these genes in iMSBmal1KO mice. (K) Heatmap presenting Log² Z-scores for Esrra and Ppargc1a in response to iMSBmal1KO and exercise training.

Figure 5

Upregulation of genes associated with pyruvate oxidation, TCA cycle, β-oxidation, mitochondrial complexes and mitoribosome subunits after exercise training in the iMSBmal1KO mouse. (A) Schematic of biological processes involved in providing acetyl-CoA from carbohydrates and lipids for the TCA cycle. Subsequently, the TCA cycle provides NADH and FADH₂ for the electron transport chain. Components are numbered/labelled for clarity of the associated genes presented in panels B and C. (B) Heatmaps presenting genes associated with pyruvate oxidation, TCA cycle and β-oxidation, (Log² Z-scores). (C) Heatmap presenting genes associated with mitochondrial complexes I–V and the mitochondrial ribosome (Log² Z-scores).

Figure 6

Loss of Bmal1 in skeletal muscle disrupts the metabolomic rewiring in response to exercise training and causes depletion of TCA cycle intermediates. (A) Table of differentially abundant metabolites (q < 0.05) with loss of skeletal muscle BMAL1, the vehicle exercise training response and the iMSBmal1KO exercise training response. (B) Venn diagram analysis of the different comparisons reveals little overlap and divergent metabolomic adaptations to exercise training. (C) Metaboanalyst KEGG metabolite pathway analysis of the iMSBmal1KO exercise training response. Numbers indicate the number of metabolites identified within the pathway. (D) Heatmap illustrating normalized metabolite abundance for amino acids. (E) Heatmap illustrating normalized metabolite abundance for metabolites relating to glycolysis and the TCA cycle. (F) Boxplots of highlighted amino acids changing with loss of BMAL1 and exercise training. (G) Boxplots of highlighted TCA cycle metabolites/phosphocreatine changing with loss of BMAL1 and exercise training. (H) Summary diagram of TCA cycle metabolites affected with exercise training in the iMSBmal1KO (Blue indicates higher abundance, red indicates lower abundance. Significance for one way ANOVA analysis is presented as ∗P ≤ 0.05, ∗∗P ≤ 0.01, ∗∗∗P ≤ 0.001, ∗∗∗∗P ≤ 0.0001.

Figure 7

Loss of skeletal muscle BMAL1 indirectly disrupts transcriptional homeostasis of peripheral tissues. (A) Upset plot of the differentially expressed genes affected by loss of skeletal muscle BMAL1 across tissues. Set size indicates the total number of differentially expressed genes regulated by loss of BMAL1 in skeletal muscle in each tissue. Numbers above bars indicate the number of differentially expressed genes. Overlapping differentially expressed genes which are common between tissues are indicated by the connected points below bars. In summary when comparing vehicle control and iMSBmal1KO controls we identified 684 (425 upregulated, 259 downregulated) differentially expressed genes in inguinal white adipose tissue (iWAT), 199 (86 upregulated, 114 downregulated) in liver, 1177 (665 upregulated, 512 downregulated) in the heart and 110 (82 upregulated, 28 downregulated) in lung tissue where identified. Pathway analysis for differentially expressed genes in inguinal white adipose tissue (B), liver tissue (C), heart tissue (D) and lung tissue (E). Red segments of bars indicate downregulation with loss of BMAL1, blue segments of bars indicate upregulation with loss of BMAL1. Associated pathway databases are indicated by (R) Reactome, (K) KEGG and (BP) Biological Processes.

Figure 8

Divergent transcriptional reprogramming of peripheral tissues in response to endurance exercise training in the absence of skeletal muscle BMAL1. (A) Upset plot of the gene sets associated with the training response of each tissue in both genotypes. Set size indicates the total number of differentially expressed genes regulated by exercise training in each tissue in both genotypes. Numbers above bars indicate the number of differentially expressed genes. Overlapping differentially expressed genes which are common between tissues are indicated by the connected points below bars. Differentially expressed genes common between genotypes within the same tissue are indicated by red connected points. (B) Summary Venn diagram of genotype specific exercise training responsive genes in vehicle (409 DEGs) and iMSBmal1KO mice (423 DEGs) with shared features (32 DEGs) in inguinal white adipose tissue. Associated inguinal white adipose tissue pathway analysis for vehicle (C) and iMSBmal1KO mice (D). (E) Summary Venn diagram of genotype specific exercise training responsive genes in vehicle (316 DEGS) and iMSBmal1KO mice (541 DEGs) with shared features (30 DEGs) in liver tissue. Associated liver pathway analysis for vehicle (F) and iMSBmal1KO mice (G). (H) Summary Venn diagram of genotype specific exercise training responsive genes in vehicle (228 DEGs) and iMSBmal1KO mice (765 DEGs) with shared features (33 DEGs) in heart tissue. Associated heart pathway analysis for vehicle (I) and iMSBmal1KO mice (J). (K) Summary Venn diagram of genotype specific exercise training responsive genes in vehicle (192 DEGs) and iMSBmal1KO mice (617 DEGs) with shared features (16 DEGs) in lung tissue. Associated lung pathway analysis for vehicle (L) and iMSBmal1KO mice (M). Red segments of bars indicate downregulation with exercise training, whereas blue segments of bars indicate upregulation with exercise training. Associated pathway databases are indicated by (R) Reactome, (K) KEGG and (BP) Biological Processes.

Figure 9

Loss of skeletal muscle BMAL1 causes rewiring of inflammatory and metabolic pathways across multiple tissues after exercise training. Deseq2 normalized counts were used to produce within tissue Z-scores for all expressed genes for all 4 experimental groups. Average Z-scores were then produced from all genes found within the given KEGG pathway for vehicle exercise trained (gray bars) and iMSBmal1KO exercise trained mice (dark blue bars) for each tissue. Tissues are indicated by green bars (Gastrocnemius Muscle), light blue bars (iWAT), yellow bars (Liver), red bars (Heart) and purple bars (Lung). Data is presented for KEGG pathways related to (A) the Immune System, (B) Endocrine System and (C) Amino Acid, Carbohydrate and Lipid Metabolism. Blue signifies relative upregulation of associated pathways and red signifies relative downregulation of the genes within the associated pathways.

Supplementary Figure 1

(A) DNA recombination gel for gastrocnemius skeletal muscle samples. Double band PCR product indicates successful recombination (loss of BMAL1). Single band indicates non-recombined Bmal1. Recombination gels for (B) heart, (C) liver, (D) inguinal white adipose and (E) lung. (F) Genome browser of exercise trained RNA-seq samples and read counts alignment across tissues with and without skeletal muscle BMAL1. Exon 8 is highlighted in red and is lost during recombination but is expressed in all other groups/tissues. iMSBmal1KO muscle does not express exon 8 leading to truncated mRNA product which is not translated. (G) Western blot for BMAL1 protein in vehicle control, iMSBmal1KO and iMSBmal1KO exercise trained mice gastrocnemius muscles and heart. The assay confirms loss of BMAL1 is specific to skeletal muscle and is not rescued by exercise training. (H) BMAL1 signal intensity (AU) relative to total protein. Full blots and total protein gels are available in Supplementary File 2.

Supplementary Figure 2

Physiological assessment of vehicle treated and iMSBmal1KO mice. (A) Differences in age (weeks) between cohorts 1 and 2 at termination. Mice in the older cohort are indicated by square data points, mice in the younger cohort are indicated by circle data points. (B) Age in weeks (termination). (C) Body weight (g). (D) MRI measured lean mass as a percentage of total body weight. (E) MRI measured fat mass as a percentage of total body weight. (F) Forelimb grip strength (N/gBW). Rotarod time to fall assessment in seconds (average of 3 trials) (G) and maximal speed (RPM) (H). Cohort 1 is delineated by circles, cohort 2 is delineated by squares. Vehicle treated mice are delineated by blue circles, iMSBmal1KO mice are delineated by red circles. Significance was determined by use of an unpaired t-test where ns indicates P > 0.05, ∗ indicates P < 0.05, ∗∗ indicates P < 0.01, ∗∗∗ indicates P < 0.001 and ∗∗∗∗ indicates P < 0.0001.

Supplementary Figure 3

Panel (A/B) Maximum work-done tests were performed prior to, after 2, 4 and 6-weeks of training with the speed adjusted accordingly to make training progressive. (C) The amount of food consumed per week was measured and is presented as a proportion of body weight. (D) Measurements of body mass were taken weekly. Habitual cage activity as measured by infrared cage monitors revealed no change in total activity levels between groups. (E) Total Activity. (F) Activity restricted to the dark/active phase. (G) Activity restricted to the light/inactive phase. Significance for two-way ANOVA analysis is presented as ∗P ≤ 0.05, ∗∗P ≤ 0.01, ∗∗∗P ≤ 0.001, ∗∗∗∗P ≤ 0.0001.

Supplementary Figure 4

Histological phenotyping of gastrocnemius muscle. (A) Gastrocnemius muscle wet weight is presented as milligrams per gram of body weight. Myovision 2.0 software was used to automatically assess the average fiber cross-sectional area (μm²) (B), the proportion of different fiber types (expressed as a percentage of all fibers) (C). (D) Gastrocnemius SDH intensity value (Log10 AU). Myovision 2.0 assessments of the total number of nuclei per fiber cross-section (E), the number of myonuclei per fiber cross-section (F). (G) illustrates the % of fibers with central nuclei and (H) illustrates the % of fibers that were embryonic myosin positive. (I) Representative dystrophin and DAPI images for each respective group. Scale bars indicate 100 μm. (J) Representative immunofluorescent images depicting whole cross-section gastrocnemius muscles from each experimental condition labeled with Type IIA/Type IIX/B. Scale bars indicate 1800 μm. Data are presented as mean ± standard deviation or min-max values for box plots. Significance for one-way ANOVA analysis is presented as ∗P ≤ 0.05, ∗∗P ≤ 0.01, ∗∗∗P ≤ 0.001.

Supplementary Figure 5

Training Study. Loss of Bmal1 in skeletal muscle causes transcriptional disruption of circadian rhythms, substrate metabolism and energy sensing pathways. (A) Volcano plot illustrating the DEGs between vehicle control and iMSBmal1KO controls. A total of 1217 (682 downregulated, 535 upregulated) genes were differentially expressed using an FDR cut off of <0.05. (B) HOMER Motif analysis performed on up and downregulated DEGs (q = 0.05), identifies a number of transcription factors regulating the transcriptional response to loss of BMAL11 in skeletal muscle. (C) KEGG Pathway analysis reveals disruption of circadian rhythms, metabolic pathways and cell cycle/energy sensing pathways with loss of BMAL1 in skeletal muscle.

Supplementary Figure 6

(A) Heatmap representing the 94 genes overlapping the vehicle control vs. iMSBmal1KO Control statistical comparison and iMSBmal1KO Control vs. iMSBmal1KO exercise trained statistical comparison.

Supplementary Figure 7

ZT13 core clock factor gene expression (Deseq2 normalized counts) in the gastrocnemius muscle in response to muscle specific Bmal1KO, exercise training and the combination of the two treatments. Note, Arntl/Bmal1 mRNA is missing exon 8 as displayed in Supplementary Figure 1 and is not translated. Data is presented as mean ± standard deviation. Significance for one way ANOVA analysis is presented as ∗P ≤ 0.05, ∗∗P ≤ 0.01, ∗∗∗P ≤ 0.001, ∗∗∗∗P ≤ 0.0001.

Supplementary Figure 8

Metabolic indices of vehicle treated and iMSBmal1KO mice with and without 6-weeks exercise training taken at ZT13. (A) Blood hemoglobin A1C percentage (Hba1c%), (B) resting blood glucose (mg/dL), (C) resting plasma insulin (ng/dL), (D) gastrocnemius muscle glycogen concentration (μg/mg), (E) gastrocnemius muscle triglyceride concentration (nmol/mg), (F) liver glycogen concentration (μg/mg), (G) liver triglyceride concentration (nmol/mg). Significance for one way ANOVA analysis is presented as ∗P ≤ 0.05, ∗∗P ≤ 0.01, ∗∗∗P ≤ 0.001, ∗∗∗∗P ≤ 0.0001.

Supplementary Figure 9

(A) HOMER transcription factor enrichment analysis of exercise responsive DEGs across genotypes and tissues presented as −Log10 P-values.

Supplementary Figure 10

Rescue correlations for all genes affected by loss of skeletal muscle BMAL1 and their exercise trained responses in iMSBmal1KO mice for (A) gastrocnemius muscle, (B) inguinal white adipose tissue, (C) liver, (D) heart and (E) lung. iWAT and heart display relatively strong relationships (R² = 0.51–0.67) suggesting exercise training in the iMSBma1KO mice brings their transcriptomes closer to vehicle control mice. Muscle, liver and lung show weaker relationships (R² = 0.04–0.22).