The mitochondrial multi-omic response to exercise training across rat tissues

Abstract

Mitochondria have diverse functions critical to whole-body metabolic homeostasis. Endurance training alters mitochondrial activity, but systematic characterization of these adaptations is lacking. Here, the Molecular Transducers of Physical Activity Consortium mapped the temporal, multi-omic changes in mitochondrial analytes across 19 tissues in male and female rats trained for 1, 2, 4, or 8 weeks. Training elicited substantial changes in the adrenal gland, brown adipose, colon, heart, and skeletal muscle. The colon showed non-linear response dynamics, whereas mitochondrial pathways were downregulated in brown adipose and adrenal tissues. Protein acetylation increased in the liver, with a shift in lipid metabolism, whereas oxidative proteins increased in striated muscles. Exercise-upregulated networks were downregulated in human diabetes and cirrhosis. Knockdown of the central network protein 17-beta-hydroxysteroid dehydrogenase 10 (HSD17B10) elevated oxygen consumption, indicative of metabolic stress. We provide a multi-omic, multi-tissue, temporal atlas of the mitochondrial response to exercise training and identify candidates linked to mitochondrial dysfunction.

Keywords: HSD17B10; acetylome; aerobic; exercise; metabolism; metabolomics; mitochondria; multi-omics; proteomics; transcriptomics.

Figure 1.. Training-induced changes in biomarkers of mitochondrial volume.

A) Experimental design. Nineteen tissues were collected from male and female rats that remained sedentary (SED) or completed 1, 2, 4, or 8 weeks of progressive treadmill exercise training. Collected tissues were assayed for epigenomics (8 tissues), transcriptomics (19 tissues), proteomics (7 tissues), PTMs (phosphoproteome 7 tissues, acetylome and ubiquitylome 2 tissues), and metabolomics (19 tissues). Mitochondria-associated transcripts and proteins were selected using MitoCarta 3.0 and mitochondrial metabolites from a previously published dataset (ref. 27). HIPPOC = Hippocampus, HYPOTH = Hypothalamus, SMLINT = Small Intestine, SKM-GN = Gastrocnemius Skeletal Muscle, SKM-VL = Vastus Lateralis Skeletal Muscle, WAT-SC = Subcutaneous White Adipose Tissue, BAT = Brown Adipose Tissue, VENACV = Vena Cava. Created using BioRender.com. B) Correlation between mtDNA quantification and the percentage of mitochondrial RNA-seq reads. Dashed line represents rho=0.5. C) Training response of biomarkers of mitochondrial volume after 8 weeks of training. Cell marks: X=not significant (p>0.05), ?= tissues in which the biomarker was not assessed. Color scale is proportional to the ANOVA-test z-score. D) Comparison of the number of significant training responses of the mitochondrial biomarkers (p<0.05). E-F) Visualization of biomarker data in SKM-GN (E), and liver (F). Each boxplot represents abundance per sex and time group. ANOVA statistics are provided for each tissue and sex combination. The whiskers extend from the hinge to the largest and lowest values, but max 1.5 x(the interquartile range). *indicates timepoint-specific significance at *p<0.05, **p<0.01 and ***p<0.001 (displayed if ANOVA p-value <0.05). The y-axis range differs for visualization purposes.

Figure 2.. The multi-omic mitochondrial response to training across tissues.

A) Heatmap of the number of mitochondria-associated analytes that significantly changed over the training time-course in at least one sex (5% FDR). Each cell represents results for a single tissue and data type. Numbers indicate training-differential mitochondrial analytes and colors indicate the proportion of differential to measured MitoCarta analytes. B) UpSet plot of the training-differential MitoCarta transcripts across tissues, identifying tissue-specific changes and differential features shared by many tissues. Numbers above vertical bars indicate transcripts in the connected tissues. Horizontal bars indicate the total number of differential transcripts per tissue. Only interaction sizes of 6 genes or more are shown. Enrichment results (MitoCarta pathways) are shown for two selected gene sets. C) MitoCarta pathway enrichments for the 8-week timepoint in the top 9 responding tissues. The 8-week male and female differential transcripts were identified using our graphical analysis. The plot shows the top pathway from each MitoCarta subcategory with the greatest number of enrichments, taken from the sum-of-log combined p-value per tissue and pathway. Each point represents a significant enrichment in a given node, where the direction of the triangle indicates the direction of effect (△= up, ▽ = down), the fill color indicates the corresponding sex (blue=male, pink=female), while a black triangle indicates sex-consistent enrichment. Because enrichment analyses are performed separately, significant enrichment can be found in multiple sets. For example, increased amino acid metabolism proteins in the heart for both the male-specific and sex consistent differential gene sets.

Figure 3.. Adrenal and brown adipose tissue demonstrate differential transcriptional dynamics in response to endurance training.

A) Graphical representation of the mitochondria-associated training-differential analytes in brown adipose tissue (BAT). Each node represents one of nine possible states (row labels, with F for females and M for males, seven states shown) at each of the four sampled training timepoints (column labels). Edges are drawn between nodes to represent the path of differential analytes over the training time-course. This graph includes the five largest paths. Both node and edge size are proportional to the number of analytes represented by the node or edge. B) Expression patterns of Ucp1 and Ucp2 in females (left) and males (right). *indicates FDR<0.05. C) Correlation between changes in Ucp2 expression and chromatin accessibility (intronic region of UCP2, chr1:165508254–165509507). Each point represents the average for n=5 animals per timepoint and sex. D) Expression changes of examples of known PPARGC1A/PGC1ɑ interactors. All are upregulated in males after 1 week and downregulated in females after 8 weeks, with exception of Jund, which is regulated in opposite directions (FDR<0.05). Data is shown as mean +/−SE for each gene per timepoint and sex. E) DREM analysis results of the MitoCarta genes in female adrenal gland predict several transcription factors to be involved in the 1-week responses.

Figure 4.. Endurance training induces largely sex-consistent increases in metabolic protein abundance in skeletal muscle.

A) The dynamics of the molecular training response visualized by constructing a summary graph in which rows represent nine combined states (row labels, with F for females and M for males, seven states shown) and columns represent the four training timepoints. Nodes correspond to a combination of time, sex, and state. An edge connects two nodes from adjacent timepoints. The differential abundance trajectory of any given training-regulated analyte is represented by drawing a path through the nodes in this graph. This graph represents the five largest trajectories for mitochondria-associated training-differential analytes in SKM-GN. Both node and edge size are proportional to the number of analytes represented by the node or edge. B) Network view of pathway enrichment results corresponding to the analytes of the week 8, sex-consistent upregulation nodes in SKM-GN and SKM-VL. Nodes indicate enriched pathways (10% FDR), and an edge represents a pair of nodes with a similarity score of at least 0.3 between the gene sets driving each enrichment. Node fill color indicates for which –ome or –omes a pathway is significant, while border color indicates if the pathway is significant in one or both tissues. Node size is proportional to the number of differential analyte sets for which the pathway is significantly enriched. Clusters of enriched pathways were defined using Louvain community detection, and are annotated with high-level biological themes. C) Fatty acid oxidation pathway enrichment for the SKM-GN proteome. Only significant genes are shown. Rows are clustered using hierarchical clustering. D) Log2 fold changes of significant differential protein phosphorylation sites in Complex I proteins. All phosphorylation changes are significant in females, whereas all except Ndufs5_T93 are significant in males at 8 weeks. Data is shown as mean +/−SE for each gene per timepoint and sex.

Figure 5.. Endurance training alters the cardiac mitochondrial acetylome.

A) Graphical representation of the mitochondria-associated training-differential analytes in cardiac muscle. Each node represents one of nine possible states (row labels, with F for females and M for males, seven states shown) at each of the four sampled training timepoints (column labels). Edges are drawn between nodes to represent the path of differential analytes over the training time-course. This graph includes the five largest paths. Both node and edge size are proportional to the number of analytes represented by the node or edge. B) Number of significantly up- and downregulated mitochondria-associated transcripts and proteins at each timepoint, with color representation based on the main MitoCarta pathway association of each analyte. C-D) Network view of pathway enrichment results corresponding to the analytes C) downregulated in both sexes after 8 weeks (the 8w_F-1_M-1 node in (A)) and D) upregulated in both sexes after 8 weeks (the 8w_F1_M1 node in (A)). Nodes indicate enriched pathways (10% FDR), and an edge represents a pair of nodes with a similarity score of at least 0.3 between the gene sets driving each pathway enrichment. Node fill color indicates for which –omes a pathway is significant, while a black border color indicates if the pathway is significant in both the down- and upregulated nodes. Node size is proportional to the number of differential analyte sets for which the pathway is enriched. Clusters of enriched pathways were defined using Louvain community detection, and are annotated with high-level biological themes. E-F) Correlation between changes in protein levels and acetylation levels in males (E) and females (F). Orange=MitoCarta proteins, gray=other proteins. G) Acetylation and phosphorylation changes (FDR<0.05) of metabolic proteins in males and females after 8 weeks (sites changing in only one sex are not illustrated). Each lollipop represents a specific acetylation (rounded top) or phosphorylation (diamond top) site, red=increase, blue=decrease. Multiple lollipops on the same protein indicates several sites changed. Hadha had 8 differentially acetylated sites, with only 6 illustrated due to space constraints. H) Site-specific acetylation changes in ACAT1 in males (top panel) and females (bottom panel), and in I) ACO2 in males (left panel) and females (right panel). All displayed sites were differentially acetylated overall (taking all timepoints and sexes into account, FDR<0.05), and sites that reach timewise significance (FDR<0.05) are highlighted with black frames.

Figure 6.. Training-induced mitochondrial adaptation in the liver through protein acetylation.

A) Graphical representation of the mitochondria-associated training-differential analytes in the liver. Each node represents one of nine possible states (row labels, with F for females and M for males, seven states shown) at each of the training timepoints (column labels). Edges are drawn between nodes to represent the path of differential analytes over the training time-course. This graph includes the five largest paths. Both node and edge sizes are proportional to the number of analytes represented by the node or edge. B-C) Correlation between changes in protein abundance and acetylation in B) male and C) female liver. Pink=MitoCarta proteins, gray=other proteins. D) Acetylation and phosphorylation changes (FDR<0.05) of mitochondrial metabolic proteins in male and female liver after 8 weeks of training (sites changing in only one sex are not illustrated). Each lollipop represents a specific acetylation (rounded top) or phosphorylation (diamond top) site, red=increase, blue=decrease. Multiple lollipops on the same protein indicates several sites changed. Proteins with more significant differential sites than could be fitted into the illustration due to space were; Atp5c 9 sites, Atp5a1 9 sites, Atp5h 13 sites, and Idh2 15 sites in total. E) Protein expression changes in Sirt3 and Sirt4. Females are represented by circles and males by triangles. Data is shown as mean +/−SE for each gene per timepoint and sex. *indicates significance at FDR<0.05. F) Site-specific acetylation changes in HMGCS2 in males (left) and females (right). All displayed sites were differentially acetylated overall (taking all timepoints and sexes into account, FDR<0.05), and sites that reach timewise significance (FDR<0.05) are highlighted with black frames.

Figure 7.. Training results in opposite regulation of mitochondrial proteins compared to type II diabetes and cirrhosis.

A) Significance of the overlap between the exercise-regulated differential proteins compared to identified proteins in case-control proteomics disease cohorts. The horizontal line represents p=0.05. MI=Myocardial Infarction, HCM=Hypertrophic Cardiomyopathy, NASH=Non-alcoholic Hepatosteatosis, Cirr=Cirrhosis, T2D=Type 2 Diabetes, HF=Heart Failure. N=8 studies with both differential proteins and the background set of all quantified proteins. Each presented dataset is named by disease and first author. B) Significance of the opposite directionality (Fisher’s exact test) when comparing the fold change sign of the overlapping proteins from (A) plus one study on NAFLD=Non-alcoholic Fatty Liver Disease. C-D) GeneMANIA protein-protein networks of proteins in opposite directions. C) Sex-consistent, week-8, skeletal muscle differential proteins that had opposite direction of effect in two separate T2D cohorts. D) Liver cirrhosis network of the 8-week female differential proteins that were both significant and had opposite direction of effect in the liver cirrhosis cohort. E) Protein abundance of HSD17B10 in HepG2 cells after siRNA-based knockdown (n=3) compared to negative controls (scramble siRNA, n=3). Number represents p-value. Representative western blot shown above bars. F-G) Oxygen consumption rate (OCR) in knockdown (n=10) versus control cells (n=9) assayed using Seahorse. Data is shown for all collected timepoints in F, with the non-mitochondrial respiration shown in the yellow section. Each calculated variable from F is shown in G. Numbers represent FDR-corrected p-values. H) Protein abundance of electron transport chain components measured using Western blot (n=3 knockdown, n=3 controls).