Molecular adaptations in response to exercise training are associated with tissue-specific transcriptomic and epigenomic signatures

Abstract

Regular exercise has many physical and brain health benefits, yet the molecular mechanisms mediating exercise effects across tissues remain poorly understood. Here we analyzed 400 high-quality DNA methylation, ATAC-seq, and RNA-seq datasets from eight tissues from control and endurance exercise-trained (EET) rats. Integration of baseline datasets mapped the gene location dependence of epigenetic control features and identified differing regulatory landscapes in each tissue. The transcriptional responses to 8 weeks of EET showed little overlap across tissues and predominantly comprised tissue-type enriched genes. We identified sex differences in the transcriptomic and epigenomic changes induced by EET. However, the sex-biased gene responses were linked to shared signaling pathways. We found that many G protein-coupled receptor-encoding genes are regulated by EET, suggesting a role for these receptors in mediating the molecular adaptations to training across tissues. Our findings provide new insights into the mechanisms underlying EET-induced health benefits across organs.

Keywords: ATAC-seq; DNA methylation; GPCR; RNA-seq; RRBS; chromatin accessibility; endurance training; sex differences; tissue specificity; transcriptome.

Graphical abstract

Figure 1

Gene expression profiling in eight tissues from control rats (A) PCA plot based on the normalized FPKM values of all expressed genes per rat tissue sample (eight tissues, five male and five female samples/tissue). Tissue type and sex are color-coded and indicated in the key. Numbers in brackets signify the proportion of variance explained by each principal component. (B) Hierarchical cluster analysis of the RNA-seq data from the eight rat tissues using 17,406 genes. Normalized gene expression values were log2 transformed and are depicted as colors ranging from red to blue as shown in the key. (C) Barplot representing the total number of expressed genes per tissue. (D) Pairwise comparison of gene expression between all eight tissues, as shown by Pearson correlation coefficients. (E and F) Scatterplots showing gene expression correlation between SKM-GN and heart (E) and between liver and hippocampus (F).

Figure 2

DNAme landscapes of rat SKM-GN and liver tissues (A) Rank order distribution of aggregated promoter DNAme levels (M-values, y axis) in SKM-GN (red line) and the corresponding values in liver (blue dots). There are 20,093 gene features per tissue (n = 10); each dot represents a gene feature. (B) Rank order distribution of aggregated GB methylation levels in SKM-GN (red line) and the corresponding values in liver (red dots). (C and D) Pearson correlation analysis between SKM-GN and liver promoter methylation levels from (A) and between SKM-GN and liver GB methylation levels from (B). (E and F) Correlation between either promoter (E) or GB (F) methylation levels and the number of CpG sites in the corresponding region in SKM-GN. (G) Hierarchical clustering heatmap of promoter and GB methylation levels for all gene features across all eight tissues. Clustering uses complete linkage and Pearson correlation. Methylation levels are displayed as colors ranging from red to blue as shown in the key. (H and J) Rank order distribution of promoter methylation levels (black line) and corresponding GB methylation levels (red dots) for either SKM-GN-expressed genes (H) or SKM-GN non-expressed genes (J). (I and K) Pearson correlation analysis between promoter and GB methylation for either SKM-GN-expressed genes (I) or SKM-GN-non-expressed genes (K). (L and M) Violin plots of promoter and GB methylation levels for either SKM-GN-expressed genes (L) or SKM-GN-non-expressed genes (M).

Figure 3

TSGs show distinctive promoter and GB methylation profiles compared to other tissues (A and B) Scatterplot of log2 FC of either promoter methylation (A) or GB methylation (B) vs. log2 FC of gene expression in liver vs. SKM-GN. The two most densely populated quadrants represent hypomethylated, highly expressed genes in either liver (NW quadrant) or SKM-GN (SE quadrant). Pearson correlation coefficients are shown. (C and D) Simple linear regression analysis of the RNA levels of liver (C) or SKM-GN (D) TSGs vs. the corresponding levels in other tissues. (E and F) Correlation between either promoter or GB methylation levels of liver TSGs and the corresponding levels in other tissues. (G and H) Same as in (E) and (F) for SKM-GN TSGs. (I and J) Simple linear regression analysis of promoter or GB methylation levels of liver TSGs vs. the corresponding levels in other tissues. (K and L) Same as in (I) and (J) for SKM-GN TSGs. (M and N) Comparison of promoter or GB methylation levels of liver TSGs vs. the corresponding levels in other tissues. (O and P) Same as in (M) and (N) for SKM-GN TSGs. (M–P) Mean ± SEM. ∗Adjusted p value < 0.0001, except in (O) for heart (<0.002) and lung (<0.001), calculated using two-way ANOVA with Tukey’s multiple comparisons test.

Figure 4

Promoter and GB methylation levels are inversely correlated with gene expression and CA across tissues (A–C) For each tissue (n = 10), expressed genes were ranked by expression level (log2 FPKM ≥ 0.5) and arranged into seven bins numbered 1 through 7 (in bold) on the x axis. In parentheses are the numbers of genes per bin. The range of log2 gene expression levels per bin is as follows: (1), ≥0.5 to <1.5; (2), ≥1.5 to <2.5; (3), ≥2.5 to <3.5; (4), ≥3.5 to <4.5; (5), ≥4.5 to <5.5; (6), ≥5.5 to <6.5; (7), ≥6.5. The x axis represents the log2 mean gene expression levels per bin ±SEM. Plotted is the correlation between RNA expression and either promoter accessibility (A), promoter methylation (B), or GB methylation (C) in each bin. (D) Pearson correlation analysis between (1) RNA expression (in green) and either promoter or GB methylation and (2) promoter CA (in blue) and either promoter or GB methylation. ∗Adjusted p value (<0.0001) from one-way ANOVA with Dunnett’s multiple comparisons test.

Figure 5

Accessible TSS proximal regions are associated with lower promoter methylation levels and tend to harbor tissue-shared motifs (A) CA at promoter regions in SKM-GN. (B) CA for the genes exhibiting ATAC-seq signal at either the TSS proximal (0 ± 200 bp), upstream TSS proximal (−200 to −3,000 bp), or downstream TSS proximal promoter region (+200 to +3,000 bp). Line graphs are color-coded by tissue type. (C and D) Distribution of the ATAC-seq peaks detected in SKM-GN (C) and liver (D) in the same promoter regions as in (B). On the y axis is the number of genes that mapped to each region set. (E–G) RNA expression (E), promoter methylation (F), and GB methylation (G) for the same genes as in (B). (B, E−G) Mean ± SEM (n = 10). ∗Adjusted p value < 0.0001; see statistics in Table S6. (H) Motif enrichment analysis in each tissue was performed on either TSS proximal, upstream TSS proximal, or downstream TSS proximal ATAC-seq peaks in a given tissue using HOMER (see the STAR Methods). Dot plot shows top TF binding motifs identified in each of the three promoter regions. The number of ATAC-seq peaks in each tissue is similar as in (B). Color indicates the enrichment log10-adjusted p value, and dot size indicates the fold enrichment. (I and J) Heatmap (I) and hierarchical clustering (J) of transcript levels for the top 30 TF motifs enriched in each tissue. Clustering uses complete linkage and Manhattan distance. Transcript levels are displayed as colors ranging from red to blue as shown in the key.

Figure 6

EET-regulated genes across tissues largely correspond to TEGs (A) Overlap between ERGs following 8 weeks of EET (n = 10) and TEGs in the age- and sex-matched control rats (n = 10). The percentage of either upregulated or downregulated ERGs overlapping with TEGs is indicated in parentheses. (B) UpSet plot showing how per-tissue ERGs intersect. Rows correspond to the sets of ERGs per tissue; columns represent the intersections between these sets. (C) Clusters of the top two over-represented pathways among downregulated and upregulated ERGs per tissue. Color indicates the enrichment log10-adjusted p value, and dot size indicates the fold enrichment. (D and E) Heatmap of GTPase- (D) and GPCR-encoding genes (E) that are differentially regulated by EET in the indicated tissues. Complete lists of the differentially regulated GTPase- and GPCR-encoding genes are provided in Tables S7 and S8. FC in gene expression values were log2 transformed and are displayed as colors ranging from red to blue as shown in the key.

Figure 7

EET elicits sex-biased transcriptional and epigenetic changes (A) Correlation analysis between DNAme signals and ATAC-seq signals in same-sex tissues across all training (1, 2, 4, and 8 weeks of EET) and control groups. F, female (n = 5); M, male (n = 5). Per-sample Pearson correlations were computed and grouped by tissue, and the correlation coefficient values were plotted in reverse scale. Boxes depict the median, 25th, and 75th percentile. Whiskers extend to the largest and smallest point up to median ± 1.5∗IQR. Points beyond that range are shown as individual outliers. (B and C) Number of differentially methylated OCRs (B) and promoters (C) in same-sex tissues following 8 weeks of EET. Differential methylation using limma voom with a false discovery rate (FDR) <0.05. (D) Shown per tissue is the number of ERGs that are either male specific (“M”), female specific (“F”), or common to both sexes (“M&F”). (E) Top two enriched pathways among ERGs that are either “M,” “F,” or “M&F” in the indicated tissue. Color indicates the enrichment log10-adjusted p value, and dot size indicates the fold enrichment. (F–I) Left, pie charts showing the number of ERGs (ERGs, upregulated ERGs [uERGs], or downregulated ERGs [dERGs], as specified) that are either “M,” “F,” or “M&F” in the indicated tissue. Right, presented are the subsets of ERGs that are enriched for the indicated pathway among all the ERGs from the corresponding pie chart (“All”) or among the ERGs that are either “M,” “F,” or “M&F.” Percentage of ERGs and significance of enrichment (q values) are indicated in parentheses. (J) LV analysis of DNAme signals overlapping with OCRs in the indicated male tissues, 8 weeks EET vs. control groups. n = 5 samples per group; t test with multiple testing correction; p value < 0.05. (K) Enriched pathways among the ERGs overlapping with the LV1, LV3, and LV9 gene sets for the indicated tissues. Color indicates the enrichment log10-adjusted p value, and dot size indicates the fold enrichment.