Multi-omic identification of key transcriptional regulatory programs during endurance exercise training in rats

Domain: biology
Published: 2026-03-21
Editorial outcome: minor_revision

Abstract
Transcription factors play a key role in regulating gene expression. We conduct an integrated analysis of chromatin accessibility, DNA methylation, mRNA expression, protein abundance and phosphorylation across eight tissues in fifty rats of equally represented sexes following endurance exercise training to identify coordinated epigenomic and transcriptional changes and determine key transcription factors involved. We uncover tissue-specific endurance exercise training associated changes and transcription factor motif enrichment across differentially expressed genes, accessible regions, and methylated regions. We discover distinct routes of training-induced regulation through either epigenomic alterations providing better access for transcription factors to affect target genes, or via changes in transcription factor expression or activity enabling target gene responses. We identify transcription factor motifs enriched among correlated epigenomic and transcriptomic alterations, differentially expressed genes correlated with exercise-related phenotypic and cell type composition changes, and training-induced activity changes of transcription factors whose target genes are enriched for differentially expressed genes. This analysis elucidates the unique gene regulatory mechanisms mediating diverse transcriptional responses to training across tissues.

Introduction
Regular exercise impacts health and mitigates disease processes throughout the body. Exercise maintains muscle function, improves cardiovascular wellness and cognitive performance, and lowers the risk of cardiovascular disease and many other disorders, including cancers and dementia1. The molecular processes mediating the adaptations induced by exercise training across tissues are poorly understood.

As the main regulators of gene transcription, transcription factors (TFs) act via the recruitment of other factors, co-activators, or co-repressors to cis-regulatory elements at the promoter or distal regions of target genes. The access of TFs to cis-motifs partly depends on chromatin organization. Hence, along with changes in chromatin accessibility and other epigenetic modifications, including DNA methylation, TFs govern constitutive gene expression as well as gene responses to stimuli in tissues. TFs are critical exercise-response mediators2,3 and, in skeletal muscle, exercise training-induced transcriptomic changes have been associated with different TFs than those induced by acute exercise4.

Our companion multi-tissue analysis of the molecular response dynamics during endurance exercise training found that the majority of differentially regulated genes are tissue-specific, whereas a small proportion are shared across multiple tissues5,6. Thus, gene responses to training are likely mediated through the combinatorial action of tissue-enriched and shared transcriptional regulators. Shared exercise-induced TF regulation can elicit tissue-specific functions, as seen with PPARγ, which is implicated in PGC1α-stimulated mitochondrial biogenesis in muscle7, regulation of adipogenesis8, and hippocampal BDNF activity and its cognitive effects9. As complex regulatory patterns drive tissue-specific gene regulation10, they likely are involved in mediating the diverse effects of exercise training on tissues. Therefore, in order to unravel the complex regulatory processes driving endurance exercise training (EET) response, we need epigenomic as well as transcriptomic data covering diverse tissues that can all be affected by exercise. However, few studies have evaluated training-induced genome-wide changes in RNA expression, chromatin accessibility, and DNA methylation11, and the majority have concentrated on skeletal muscle.

We leveraged the study design of the Molecular Transducers of Physical Activity Consortium (MoTrPAC) EET study in rats5 to characterize the TFs mediating gene responses to training across multiple tissues. During 8 weeks of EET, genome-wide transcriptome, chromatin accessibility, and DNA methylation were assayed in 8 tissues from age-matched male and female rats: gastrocnemius (SKM-GN), heart, hippocampus (HIPPOC), kidney, liver, lung, brown adipose (BAT), and white adipose (WAT-SC). Protein abundance and phosphorylation data was assayed in six of these tissues (SKM-GN, heart, kidney, liver, lung, and WAT-SC) and we identified relationships between EET responses in TF activity at the proteomic level with changes in TF target EET response. Acronyms and abbreviations used in this paper are summarized in Table 1. Through integrative analysis of all omes for each tissue, we establish a map of the regulatory transcriptional responses to training across tissues.

Results
Characterization of epigenetic and transcriptional responses to endurance training
To understand the epigenetic and transcriptional response mechanisms elicited during eight weeks of EET, we analyzed ATAC-seq, RNA-seq and RRBS profiles generated in skeletal muscle (gastrocnemius; SKM-GN), heart, hippocampus (HIPPOC), kidney, liver, lung, brown adipose tissue (BAT), and subcutaneous white adipose tissue (WAT-SC) from fifty rats divided into groups of 10 animals (5 males, 5 females). Animals were subjected to either 1, 2, 4, or 8 weeks of training (Fig 1a) or to no training (untrained controls), consisting of 8 weeks of housing and stationary treadmills. Protein abundance and protein phosphorylation were also measured in six of the eight tissues: SKM-GN, heart, kidney, liver, lung, and WAT-SC. We identified differentially accessible regions (DARs; F test adjusted p-value < 0.1), differentially methylated regions (DMRs: F test adjusted p-value < 0.1), and differentially expressed genes (DEGs; F test adjusted p-value < 0.1) between EET and control groups (Fig 1b). To characterize the transcriptional and epigenomic changes induced by EET across tissues, we evaluated the tissue-specificity of DARs, DMRs, and DEGs. Although most expressed genes, open chromatin sites, and methylated CpG sites were detectable in multiple tissues, the majority of DARs (90%), DMRs (91%), and DEGs (66%) were identified in only one tissue (Fig 1c, Fig. S1). This suggested that gene regulatory responses to EET were largely specific to individual tissues, which was in line with another MoTrPAC publication6.

We then examined the distribution patterns of log2 fold change (LS2FC) in gene expression, L2FC in chromatin accessibility, and L2FC in methylation across time points and sexes (Fig. S2, S3, S4 respectively). While the ratios of up- to down-regulated analytes (i.e. DEGs, DARs, DMRs) were similar across the majority of tissues, L2FC patterns differed between tissues and omes. Overall, methylated regions showed the largest spread of L2FC values compared to expressed genes and chromatin accessibility peaks, a pattern that was even more pronounced when comparing DMRs to DEGs and DARs. DEGs and DMRs were more likely to have larger L2FC values than non-significant features in their respective omes; however, DARs did not share this trend. Multiomic power analysis of the five omes included in this study identified RNA-seq as the ome with a consistently higher power and the phosphoproteome with consistently the least power (Fig. S5).

Sex differences were the most frequent separator of samples in most tissues along top PC axes in RNA-seq (Fig. S6), ATAC-seq (Fig. S7), and RRBS (Fig. S8). We observed sex differences among DEGs in HIPPOC, BAT, and WATT-SC. Sex differences were also prevalent among DARs and DMRs in a majority of tissues. Overall, the proportion of DEGs showing concordant changes (L2FC) across time points and sexes was higher than that of DARs and DMRs (Fig. 1d). Notably, heart, SKM-GN, and kidney exhibited the most consistency in expression changes across all groups, whereas WAT-SC showed the least. In BAT, we detected pronounced variations in DEG as well as in DAR profiles between earlier and late time points in both sexes (Fig. S2g and S3g). Gene set enrichment analysis among week 1 or week 8 DEGs in males and females highlighted varying pathway enrichment patterns that were, in most cases, consistent across sexes and time points (Fig S9). SKM-GN and heart shared enrichment for oxidative phosphorylation and cardiac muscle contraction pathways, as well as markers for Parkinson’s, Huntington's, and Alzheimer's diseases. WAT-SC DEGs were enriched for the chemokine signaling pathway and immune-related diseases, including systemic lupus erythematosus, asthma, and primary immunodeficiency.

Next, we investigated whether alterations in cell type proportions contributed to the expression changes observed. Cell type deconvolution analysis (See Methods) identified large-scale changes in immune cell type proportions that were attributable to training duration in BAT (Fig. 1e and Fig. S10a) and to sex independent of training in WAT (Fig. 1f and Fig. S10b). Changes in dendritic cell proportion in heart, fibro/adipogenic progenitor (FAP) cell proportion in kidney, and monocyte and endothelial cell proportions in lung were also notable. Significant percentages of DEGs in many tissues were either correlated or anti-correlated with changes in immune cell type composition in response to EET (Fig. S11), with the most consistent correlations occurring in BAT or WAT-SC. Large percentages of DEGs were also correlated with changes in natural killer (NK) and granulocyte/monocyte progenitor (GMP) cells in liver and monocytes in lung. A full list of correlations between DEGs and cell type compositions can be found in Supplemental Data S1. Because the immune cell type composition differences observed between female and male WAT-SC were independent of exercise training, and because differential analysis was conducted in male and female samples separately, the DEGs identified in WAT-SC in each sex do not likely reflect any cell type composition changes occurring during training.

To further characterize the epigenomic changes induced by EET across tissues, we examined the genomic distribution of DARs vs. that of all open chromatin regions detected (Fig. 1g,h). Compared to all accessible regions, DARs were significantly enriched at proximal promoters across all tissues except WAT, which was excluded due to a scarcity of DARs (Fig 1h,k, Methods). Consistent with previous studies12,13, open chromatin peaks across tissues were predominantly located in intronic and distal intergenic regions. Given the importance of the proximal promoter in the regulation of gene transcription14, the enrichment of DARs in this region suggested that EET results in the transcriptional activation of target genes. DMRs were significantly enriched at downstream regions across all tissues except WAT, relative to all methylation sites in each tissue (Fig 1i,j,l). Unlike DARs and DMRs, open chromatin regions that mapped to DEGs, which we refer to as DEG-associated peaks (DEGaPs), shared a similar genomic distribution as the peaks associated with all expressed genes (Fig S12); however, ATAC-seq peaks that mapped to either expressed genes or DEGs contained a higher proportion of intronic peaks and a lower proportion of distal intergenic peaks compared to all open chromatin peaks in the dataset. This was in part explained by a subset of unexpressed genes only containing distal intergenic peaks, suggesting that some gene-body chromatin accessibility was likely to be important for gene expression.

Differential analytes in a tissue also tend to have higher baseline enrichment in that tissue (Figure 2). Upon calculating the z scores of baseline expression of DEGs, accessibility of DARs, and methylation of DMRs across all 8 tissues, we found, in the case of both RNA-seq and ATAC-seq data, DEGs and DARs in a given tissue were, on average, more highly expressed (Figure 2a,d) and more highly accessible (Figure 2b,e), respectively, at baseline in the tissue with the significant training response. Conversely, DMRs in a given tissue were more hypomethylated at baseline in the tissue with a significant training response (Figure 2c,f). This pattern of tissue specificity was mostly maintained when subdividing differential analytes into up- and down-regulated subsets (Fig S13); however, down-regulated DMRs were less likely to follow this pattern, most notably seen with BAT. While there is some overlap in significant training responses across tissues, the majority of training response analytes represent tissue-enriched analytes, even if they are not tissue-specific.

Identification of distal correlated epigenetic regulation despite few DARs and DMRs mapping to adjacent DEGs
We next sought out DAR-DEG associations by assessing the concordance between chromatin accessibility and gene expression changes. We assigned each DAR to the nearest gene and determined the fraction of DARs that were annotated to DEGs. Applying a hypergeometric test, we found that BAT, SKM-GN, and liver showed a higher proportion of overlap between DARs and DEGs (Fig 3a). The substantial overlap between DARs and DEGs in BAT may be related to the EET-induced increase in immune cell populations (see Fig. 1e). SKM-GN and liver showed the highest count of DARs among all tissues (Fig. 1b). Despite hundreds of DARs in both kidney and lung, only a few of their nearest genes were DEGs (Fig 3a). Similarly, we investigated DMR-DEG associations by assigning each DMR to the nearest gene and determining the fraction that was annotated to DEGs. Only BAT and WAT-SC have significant overlaps between DMRs and DEGs (Fig 3b). BAT and WAT-SC are the two tissues with the largest number of DMRs in the dataset by a considerable margin (Fig 1b).

The binding of TFs to distal open chromatin regions can regulate gene transcription15,16. Given the modest proportion of DARs mapped to adjacent DEGs, we extended the search window and sought relationships outside the closest gene for a given DAR. With respect to the location of DARs relative to the TSS of the nearest DEG, in all tissues, the majority of nearest DAR-DEG pairs reflected a normal distribution with a median centered approximately 1 Mb away from the nearest DEG, and a substantial left tail representing closer pairs (Fig 3c). BAT, SKM-GN, and liver contributed most of the DARs proximal to a DEG, confirming our earlier observations (Fig 3a). A similar pattern is seen when measuring the distance between a DMR and the nearest DEG (Fig 3d) and the distance between a DMR and the nearest DAR (Fig 3e). BAT and WAT-SC contain the majority of closest DMR-DEG pairs likely because of their greater DMR-DEG overlap (Fig 3b) and BAT contains the majority of closest DMR-DAR relationships. WAT-SC DMR-DAR relationships tend to be the most distant, largely because of the lack of DARs, reducing the likelihood of neighboring interactions.

DARs that were adjacent to DEGs tended to be more highly correlated with gene expression changes (Fig 3f, Fig S14), a pattern predominantly seen in SKM-GN and liver. In contrast, we found less of a connection between distance to TSS and correlation in training response for DMRs adjacent to DEGs, although there were still a large number of highly correlated adjacent pairs, especially seen in BAT and WAT-SC (Fig 3g, Fig S15). Interestingly, while instances of adjacent DMRs and DARs are more limited in this study, in tissues more populated by DAR-DMR pairs, positive correlation tended to occur when the DMR was upstream of the DAR, and correlation was decreased or negative when the DMR was downstream of the DAR (Fig 3h, Fig S16), most evidently seen in liver.

DAR-DEG pairs are associated with distinct pathways in each tissue and SP2, BMYB, and BMAL1 represent key regulatory TFs
We focused on DARs and DEGs located within 500 kb from each other and isolated those that were either positively or negatively correlated across samples at all time points and sexes (t test of Pearson correlation p-value < 0.05; Supplementary Data S2). Pathway enrichment analysis among DAR-DEG pairs identified largely tissue-specific enrichment patterns for each tissue (Fig S17). In agreement with the training-associated increase in immune cell types inferred from cell type deconvolution analysis (see Fig. 1e), DAR-DEG pairs in BAT showed enrichment for several immune pathways. Lung also showed considerable enrichment for immune-associated pathways, suggesting activating roles in antigen defense by exercise. By contrast, liver DAR-DEG pairs were primarily enriched for lipid biosynthesis and metabolic processes, while heart DAR-DEG pairs were enriched for muscle movement and filament sliding, and SKM-GN DAR-DEG pairs were enriched for myofiber synthesis and muscle contraction (Fig S17). These results suggested that the correlated epigenetic and transcriptional changes induced by training affected tissue-specific functions.

To identify key regulators of the training response in each tissue, we analyzed TF motif enrichment at the DARs of DAR-DEG pairs (Fig. 4a). Notably, we identified SP2 in lung as the most significantly enriched TF with 33.33% of DAR-DEG pairs in lung containing an SP2 binding site vs. a general enrichment of 4.35% (p-value = 0.00032). The DAR-DEG pair targets of SP2 exhibited all negative correlations in training response (Fig. 4b) between gene expression and chromatin accessibility and included the genes Mpeg1 (Fig. 4e), Gm2a (Fig. 4f), Bcl11b (Fig. 4g), Dennd1c (Fig. S18a), and Oas2 (Fig. S18b). Each gene and associated peak shared a similar pattern of training response as well, reflecting a decrease in gene expression and an increase in chromatin accessibility at early time points followed by a return to baseline. In SKM-GN, BMYB was the most significantly enriched TF, with 18.42% enrichment among DAR-DEG pairs vs. 3.67% general enrichment (p-value = 0.00042). Most DAR-DEG pair targets of BMYB exhibited a positive correlation in training response (Fig. 4c), including Hjv (Fig. 4h), Fkbp5 (Fig. 4i), Ppard (Fig. 4j), Cited4 (Fig. S18c), Cab39 (Fig. S18d), Sik1 (Fig. S18e). Phlda3 was the only target that exhibited a negative correlation between gene and peak (Fig. S18f). A female-specific response early in EET is shared among Hjv, Fkbp5, and Ppard.

The circadian clock factor BMAL1 was the most enriched TF in liver with binding sites present in 14.47% of DAR-DEG pairs vs. 5.42% general enrichment (p-value = 0.0026). Most DAR-DEG pair targets of BMAL1 were positively correlated in training response (Fig. 4d), of which Abhd2 (Fig. 4k,l) and G6pc1 (Fig. S18g,h,i) were correlated to multiple neighboring DARs. Other positively correlated targets include Chp1 and Onecut1 (Fig. S18j,k), while Lpar3 was negatively correlated with a neighboring DAR (Fig. S18l). MAZ was identified as significantly enriched for DAR-DEG pairs in both SKM-GN (10.52% vs. 3.43%, p-value = 0.0398) (Fig. S18m) and lung (26.67% vs. 3.43%, p-value = 0.0014) (Fig. S18n). MAZ targets in SKM-GN included both positively and negatively correlated DAR-DEG pairs, including targets Kcna7 and Eif4g3, which were positively correlated with neighboring DARs (Fig. S18o,p) and Srsf2, which was negatively correlated (Fig. S18q). MAZ targets in lung were all negatively correlated between gene and peak. They also overlapped with the gene targets of SP2 in lung, except for Hspb6 (Fig. S18r), likely because of similarities in MAZ and SP2 binding sites.

Correlated DMR-DEG pairs enriched for key TFs in adipose and lung tissue
We identified DMR-DEG pairs located within 500kb of each other and whose L2FC training responses were either positively or negatively correlated across time points and sexes (t-test of Pearson correlation < 0.05, Supplementary Data S3). The majority of DMR-DESG pairs were found in WAT-SC, 194 of 245 in total, while 26 were found in lung, 9 in heart, 7 in SKSM-GN, 4 in liver, 3 in kidney, and 2 in hippocampus. The heavy skew towards WAT-SC is accounted for by the increased numbers of both DEGs and DMRs found in the tissue. We identified 14 TFs whose motifs were significantly enriched among the DMRs in correlated neighboring DMR-DEG pairs in each tissue (Fig. 4M). Seven enriched TFs were found in WAT-SC, while 5 were found in lung, and two in SKM-GN. SP2 was the most enriched TF in SKM-GN with motifs found in 57.14% of DMR-DEG pairs vs. 13.65% general enrichment (p-value = 0.00863). As was seen with DAR-DEG pairs in lung, SP2 targets were all negatively correlated between gene expression and DNA methylation in their training response (Fig. 4N). Zinc finger Zdhhc1 was negatively correlated with a DMR that included an SP2 binding site (Fig. S19d), while Slc26a10 was negatively correlated with two DMRs containing SP2 binding sites. NF1-halfsite was the most significantly enriched TF in WAT-SC with motifs found in 21.64% of DMR-DEG pairs vs. 8.08% general enrichment (p-value = 0.000027) (Fig. S19a). Targets included Erg (positive correlation, Fig. S19e), Tand Kcnj15 (negative correlation, Fig. S19f). The DMR chr2-375476_cluster1 was correlated with five DEGs, one positively correlated: Syt11 (Fig. S18g), and four negatively correlated: Crabp2 (Fig. S19h), Paqr6 (Fig. S19i), Tmem79 (Fig. S19j) and Rab25. NF1-halfsite was also significantly enriched in lung with motifs found in 33.33% of DMR-DEG pairs vs. 8.08% general enrichment (p-value = 0.0124). The TF AP-2gamma was also identified as significantly enriched in lung with motifs found in 33.33% of DMR-DEG pairs vs. 9.13% general enrichment (p-value = 0.0189). There was considerable overlap in DMR-DEG targets for both NF1-halfsite and AP-2gamma including the DMR chr7-140730_cluster4 which is correlated with 4 DEGs, three positive correlations: Agap2 (Fig. S19k), B4galnt1 (Fig. S19l), Arhgap9 (Fig. S19m), and one negative correlation: Avil (Fig. S19n). DMR chr4-156530_cluster11 was identified as a target of AP-2 gamma, but not NF1-halfsite. This DMR was negatively correlated with two transmembrane protein encoding DEGs: Tmem176a (Fig. S19o), and Tmem176b (Fig. S19p).

Characterization of TF expression responses to EET
As putative transcriptional regulators were inferred from DAR-DEG and DMR-DEG correlations in a restricted number of tissues, we sought to independently characterize TF expression responses to EET per tissue over the 8-week training period. We measured the RNA levels of all TF-encoding genes and assessed their protein abundances and phosphorylation levels based on mass spectrometry data from a subset of 6 tissues. Various subsets of TFs exhibited significant changes at the transcriptome (Fig 5a), proteome (Fig 5b), and phosphoproteome (Fig 5c) levels in each tissue. Most significant TF changes were tissue-specific (Fig. S20); however, overlaps at the transcriptomic level were identified between SKM-GN and heart, and between BAT and WAT-SC, and at the phosphoproteomics level between heart, liver and lung. BAT and WAT-SC had the largest and most significant changes in TF gene expression (Fig 5a), which included immune-related TF genes such as Irf8 and Pou2f2 responding in both tissues and Irf1 and Spi1 specifically in BAT. Transcript levels of Egr1 decreased significantly across multiple tissues, including SKM-GN, heart, kidney, and lung. Fellow early response gene Fos also decreased across the same tissues, reaching significance in kidney. Myb expression was significantly altered both in HIPPOC, where transcript levels increased in response to training in both sexes, and in kidney, where transcript levels increased in females but decreased in males.

WAT-SC and lung exhibited the most significant changes in TF protein abundance (Fig. 5b). Select TFs showed significant protein level changes in multiple tissues, including RORC which decreased in female subjects in both lung and kidney. PBX1 decreased in WAT-SC while PBX2 decreased in kidney, and ATF1 and ATF3 increased in male subjects in WAT-SC and heart, respectively. NR4A1 (NUR77) and MEF2A levels decreased in SKM-GN at week 8 of training, while MEF2C levels increased at week 2 of training.

MEF2C showed significant protein phosphorylation changes in SKM-GN, lung, and WAT-SC (Fig 5c). Other MEF family members, MEF2A and MEF2D, also exhibited significant changes in phosphorylation in lung and SKM-GN, respectively. A number of TFs had multiple significant phosphosite changes in a single tissue, including two NF1 phosphosites in heart, two STAT3 sites and two ATF2 sites in liver, two RFX5 and two MEF2C sites in lung, and two USF2 sites and four TRPS1 sites in WAT-SC. The majority of significant EET responses in TFs were ome-specific, with a few exceptions primarily in WAT-SC (Fig. S21). SREBF1, IRF8, FOXO3, FOXK1, and TCF7L2 all had significant EET responses in WAT-SC and were positively correlated in their transcriptomic and proteomic responses. We also identified positive correlations in RNA-Protein EET response in STAT4 in lung, RXRA and STAT3 in liver, and a positive correlation in RNA-Phosphoprotein EET response in CEBP6 in heart.

Now that we determined which TFs were found to have significant training responses at either a transcriptome, proteome or phosphoproteome level, we wanted to find which of their motifs were enriched among the promoter regions of DEGs, identifying potential training-induced regulatory relationships. We found a number of TFs whose motifs were enriched among DEGs (Fig 5d, Fig S22-S23), including eight statistically significant motif enrichments.

SF1 (Fig. S24a), SIX1 (Fig. S24b) and SIX2 all exhibited a significant training response at the transcriptome level in SKM-GN and their motifs were significantly enriched in target DEGs. The FOS motif was not significantly enriched in target DEGs in kidney; however, four of its target DEGs exhibited a similar EET response to FOS itself (Fig. S24c). A JUND motif was enriched in target DEGs in SKM-GN (Fig. S24d), heart (Fig. S24e), and WAT-SC (Fig. S24f), a PBX1 element was enriched in lung (Fig. S24g) and RAR:RXR was enriched in BAT (Fig. S24h). MEF2C had significant EET responses in SKM-GN at both the proteome and phosphoproteome level and its motif was significantly enriched in target DEGs (Fig. S25a,e). IRF:BATF was significantly enriched in DEG targets in lung while the heterodimer also exhibited a significant EET response at the proteome level (Fig. S25b). Other TFs with significant EET responses in protein abundance and motif enrichment in DEG targets included NR4A1 (Nur77) in SKM-GN (Fig. S25c), and PBX1 in WAT-SC (Fig. S25d). Dll4 was found to be a DEG target of PBX1 in both WAT-SC and lung. MEF2A responded significantly to EET at the phosphoproteome level and its motif was significantly enriched in target DEGs in heart (Fig. S25f), while binding sites for ATF7 in kidney (Fig. S25g) and STAT1 in liver (Fig. S25h) were also enriched in DEG targets.

SIX1 DEG targets in SKM-GN exhibit a range of functions and training responses, including collagen gene Col3a1, and muscle-contraction associated gene Lmod1, both of which decreased in response to training. Protein modification-associated genes Golga4 and Art1 increased over training, while ubiquitin gene Usp2 decreased. Malic enzyme Me3, which is involved in the oxidative decarboxylation of malate to pyruvate, is consistently expressed at higher levels during training. MEF2C protein levels significantly increased at week 2 of training and then returned to baseline, while phosphorylation increased significantly throughout the eight weeks of training. Among MEF2C target genes, clock gene Per1 demonstrated the highest increase in expression at the onset of training, while Sema6c, Ankh, Ptpn1, and Phkg1 exhibited decreased expression following training. Dystrophin (Dmd), a critical protein for muscle fiber integrity24, was the most negatively correlated with MEF2C protein level changes, but positively correlated with MEF2C phosphorylation changes. IRF:BATF targets in lung predominantly exhibited decreased expression following training, including tubulin gene Tuba1c, mitochondrial biogenesis-associated gene Perm1, and actin cytoskeleton organizational gene Cfl1. NR4A1 (NUR77) significantly decreased in protein level in SKM-GN, as did the majority of its DEG targets including heat shock protein Hspa1l and dual specificity phosphatase Dupd1.

We calculated the sample-level correlation in EET response between significant TFs and their target genes, and identified TFs whose DEG targets were most likely to be correlated to TF response (Fig. S26). NR4A1 and SIX2 in SKM-GN and JUND in heart were significantly correlated with DEG targets (Fig. S26a,b); however, the strongest correlations were found in lung (Fig. S26e). TFs PU.1 and IRF8 were found to be the most significantly correlated with DEG targets, with IRF4 and RUNX also achieving significance. IRF8, RUNX, and Spi1, the gene encoding the PU.1 TF, all share similar patterns of EET response, namely a decrease in expression most pronounced in female subjects and early training time points. Indeed, this pattern is positively correlated with the decreased proportion of monocytes identified in lung by cell type decomposition analysis (Fig. 5e). IRF8 is known to act on monocyte differentiation, and CBFB and RUNX have been identified as upstream regulators of IRF8 activity25. Indeed, CBFB knockouts in mice blocked IRF8 expression and monocyte differentiation, with monocyte activity being restored upon reintroduction of IRF8. Previous mouse work showed that RUNX1 also acts as an upstream regulator of PU.126. These regulatory relationships have not been previously characterized in EET response. However, our data reveal positive sample-level correlations in EET response between IRF8 and CBFB protein amounts (Fig. 5f), IRF8 and RUNX1 protein amounts (Fig. S27a), and between RUNX1 protein and PU.1 (Spi1) transcript levels (Fig. S27b). The majority of PU.1 DEG targets (Fig. S27c) and IRF8 DEG targets (Fig. S27d) show similar correlations with their respective TFs. In light of prior knowledge and our current findings, we propose a pathway involving the activation of IRF8 by RUNX1 and CBFB, the activation of PU.1 by RUNX1, and, via their gene targets, the induction of differentiation of myeloid progenitors into monocytes (Fig. 5g). In our study of EET response, we found this pathway inhibited at the level of RUNX1 and CBFB (Fig. 5f,S27a), potentially producing a cascade of inhibition that leads to a decrease in monocyte abundance in the lung in females following exercise training.

DARs vs. DMRs vs. DEGaPs show distinct TF motif enrichment patterns that differentially correlate with TF gene expression
The lack of nearby DARs or DMRs for the majority of DEGs within each tissue (see Fig. 3a,b) led us to hypothesize that DARs, DMRs, and DEGaPs may mediate different paths of transcriptional regulation: i) DARs and DMRs coordinating a combination of gene-proximal and gene-distal cis-regulatory interactions, ii) a combination of statically open cis-regulatory elements (DEGaPs) and changes in TF behavior influencing differential gene expression. To address this, we analyzed TF binding site enrichment at either DARs or DEGaPs relative to all open chromatin peaks in each tissue, or at DMRs relative to all methylated CpG sites in each tissue (Fig. 6a,b,S28; see Methods).

Due to their very low number of DARs, WAT and HIPPOC were removed from the analysis (Fig. 6a). Motif enrichment patterns varied greatly across the six remaining tissues, reflecting a high degree of tissue-specificity. Motifs for both FOX and KLF families of TFs were over-represented among SKM-GN DARs. KLFs are zinc-finger TFs that have been associated with myogenesis and muscle fusion via their recruitment to Muscle Creatine Kinase (MCK) promoters27. Motifs enriched in DARs included SIX2 and MEF2C in heart, COUP-TFII in kidney, and SP2 in liver. HOXA10 and HOXD12 motifs were enriched in lung, while both IRF8 and IRF3 motifs were enriched in BAT. Tissue-specific patterns of motif enrichment were maintained when measuring the frequency of motif presence in DARs across tissues (Fig S29a).

A similar pattern of tissue-specific motif enrichment was observed among DMRs for each tissue (Fig S28a), albeit with different TF motifs than in DARs. A MYOG motif was the most significantly enriched among DMRs in SKM-GN, along with binding motifs for two neuronal TFs NEUROD1 and NEUROG2. Another motif for nervous system-associated TF NKX2.2 was the most significantly enriched in heart. A PRDM9 motif was the most enriched in HIPPOC, as was PAX5 in kidney, and multiple HOX TF motifs in liver including HOXA10 and HOXD12, which were also enriched among DARs in lung. TF motif enrichment significance among DMRs overlapped more between lung, BAT and WAT-SC with ELK4, ELK1, ETV1, and ETV4 being enriched across the three tissues. GLIS3 motif enrichment was restricted to WAT-SC. As with DARs, the frequency of motif presence in DMRs across tissues maintained tissue specificity and even more sharply distinguished lung, BAT and WAT-SC tissues (Fig S29b).

Motif enrichment patterns among DEGaPs differed considerably from those in DARs and DMRs (Fig. 6b, S29c). Indeed, pairwise motif enrichment comparisons between DARs, DMRs, and DEGaPs per tissue were weakly correlated (<0.39; Fig S30, Fig S31, Fig S32). Motif enrichment significance in DEGaPs was greater in lung, BAT, and WAT-SC, presumably due to their higher proportions of DEGs (see Fig. 1b). MEF2 TF motifs were enriched across lung, BAT, and WAT-SC, as well as in SKM and heart, forming one cluster; on the other hand, ETS and EPLF TF motifs were more exclusively enriched in lung, BAT, and WAT-SC and formed another cluster. MEF2 TFs are typically involved in muscle tissue regeneration28. ETS and ELF TFs are associated with the regulation of immunity29,30, suggesting that they may be related to the immune cell type composition changes occurring in adipose tissues. With respect to tissue-specific motif enrichment, sites for immediate early gene-encoded JUN and FOS TFs were enriched in HIPPOC (Fig S29c), HNF1 and ZIC families of TFs in kidney, and FOX and SOX families of TFs in liver.

To establish correlations between the motifs enriched at either DARs, DMRs, or DEGaPs and the relative expression levels of the corresponding TFs, we examined TF gene expression patterns in control tissues (Fig S33a, Fig S33b, Fig S33c, respectively) and L2FC in TF gene expression following EET (Fig S34a, Fig S34b, Fig S34c, respectively). We found a stronger correlation between motif enrichment and TF control gene expression levels among DEGaPs (Fig S35g-i) than among DARs (Fig S35a-c), and consistently found a negative correlation consistently among DMRs (Fig S35d-f). Conversely, there was no correlation between motif enrichment and L2FC in TF gene expression, be it among DARs (Fig. S29a,S34a), DMRs (Fig. S29b,S34b), or DEGaPs (Fig. S29c,S34c). These findings suggest that the influence of DEGaPs on DEG expression is more dependent upon the expression level of the associated TFs, whereas DARs and DMRs can influence DEGs directly, supporting the concept that DEGaPs vs. DARs and DMRs mediate distinct paths of transcriptional regulation.

DARs, DMRs, and DEGaPs show cross-tissue motif enrichment conservation at specific genomic regions
The different genomic distributions of DARs and DMRs (see Fig. 1h,j) also suggested that some regions may contribute more to the regulation of EET gene responses than others. Thus, following the mapping of tissue DARs, DMRs, and DEGaPs to various genomic regions (i.e. promoter, upstream, downstream, 3’UTR, 5’UTR, intron, exon, distal intergenic), we examined motif enrichment within each subcategory of DARs, DMRs, and DEGaPs. We found that tissues shared conserved patterns of TF motif enrichment at the DARs mapped to the proximal promoter, downstream, and 3’UTR regions (Fig 6c). Proximal promoter DARs showed strong enrichment for SP and KLF TF motifs, MAZ and PITX1. Downstream DARs were enriched for binding sites for NPAS and BMAL1, two core circadian clock TFs, as well as NKX TFs. The NF1 subcluster also showed enrichment among 3’ UTR DARs, while THRB was most highly enriched among 3’ UTR DARs.

In contrast, TF motifs were most highly enriched among downstream DMRs relative to other genomic regions (Fig S28b), including multiple KLF TFs, MYOG, MAZ, and EGR1. ETV2 and ETS1 motifs were enriched among proximal promoter DMRs, while BM2 and PIT1 were enriched among upstream DMRs. Another cluster of TF motifs were enriched among exon DMRs including c-MYC, HIF2A, and MYOD. TF motif enrichment among DEGaPs was most significant in intronic regions, one of the regions most populated with DEGaPs (Fig 6d). However, there was no enrichment among distal intergenic DEGaPs despite their large proportion of total DEGaPs. TF motifs enriched among intronic DEGaPs included SCL, SMAD3, THRB and ERRA. Significant enrichment was also found among exonic DEGaPs, which was highly correlated to intronic DEGaP enrichment. While also highly enriched among intronic DEGaPs, the most enriched TF motifs among promoter DEGaPs included MAZ and multiple KLF family members. Overall, genomic features appeared to be a driving force in clustering the enriched binding motifs among DEGaPs, DARs, and DMRs.

Differential TF motif enrichment at DEGaPs of upregulated vs. downregulated DEGs
An associated MoTrPAC paper6, found that the biological pathways that were over-represented among up-regulated vs. down-regulated DEGs in each tissue after 8 weeks of training were distinct, suggesting that different sets of TFs were involved in their regulation. We sought to predict the sets of TFs that preferably bind to up- versus down-regulated genes. After identifying the DEGaPs related to up-regulated DEGs and those linked to down-regulated DEGs in each tissue, we refined each DEGaP category into distal intergenic, intron, and promoter peak subsets (the most prevalent genomic features among peaks, as shown in Fig. 1g), and determined their respective patterns of TF motif enrichment. We then generated z scores for each TF enrichment across the peak sets, and applied hierarchical clustering to the TFs. We found that in most tissues, clustering was mainly driven by differences in TF motif enrichment between the promoter peaks associated with up-regulated DEGs vs. those associated with down-regulated DEGs (Fig S36). Adipose tissues had the greatest motif enrichment overlap between promoter peaks associated with up-regulation and those associated with down-regulation.

Restricting the TF motifs enriched within each subset of promoter peaks (i.e. promoter peaks associated with either up- or down-regulation) to those conserved across all four time points (Fig. S37), revealed some overlap in motif enrichment between the tissues (Fig. 7a,b). For instance, CLOCK, BHLHE41, and MYC (c-Myc) were among the TF motifs enriched across WAT-SC, heart, lung, SKM-GN, and kidney at the promoter peaks associated with up-regulation; another set of TF motifs enriched primarily in heart and lung included ZFX and EBF2, while ATF2 was enriched in SKM-GN and kidney. Among the promoter peaks associated with downregulation, there was less consistent enrichment across tissues. However, the NKX2 motif was enriched across SKM-GN, kidney, liver, and lung, and the SOX9 binding site was enriched among 6 of the 8 tissues. Twenty-seven TFs were shared between the 50 TF motifs associated with up-regulation (Fig. 7a) and the 39 associated with down-regulation (Fig. 7b), suggesting their involvement in both positive and negative regulation depending on the tissue. Heart, SKM-GN, and kidney shared the most enriched TF motifs associated with up-regulation and lung, liver, and SKM-GN shared the most enriched motifs associated with down-regulation (Fig 7c,d). Overall, while gene responses and TF expression responses were largely tissue-specific, tissues showed an overlap in the sets of TF motifs associated with either up-regulated or down-regulated DEGs, suggesting that specific patterns of EET-induced regulation could be conserved across tissues.

Connecting TFs to DEGs associated with EET-induced phenotypic changes
Exercise elicits various phenotypic changes such as increased aerobic capacity (VO2max) and reduced body fat mass, all of which reflect physiological adaptations. As several parameters were clinically measured throughout the EET period, we sought to: i) assess possible correlations between phenotypic alterations and gene expression responses to training, ii) infer key transcriptional regulators based on their motif enrichment at the accessible promoters of DEGs that correlated with phenotypic changes.

Parameters measured included weight, body fat mass, body lean mass, body water, lactate levels, and VO2 max (Fig. S38). Body weight and fat mass were significantly lower in 8-week trained rats than in controls (Fig S38a,b), with the greatest difference observed in males. VO2 max significantly increased in both females and males in response to EET (Fig S38c). Body lean mass also increased in both sexes, though significance was only reached in females (Fig S38d). Lactate (Fig S38e) and body water (Fig S38f) showed no significant changes relative to control, though a significant decrease in lactate levels was observed between week 4 and week 8 in males. Positively correlated parameters included VO2 max vs. body lean mass, and body weight vs. body fat mass (Fig 8a).

In most tissues, DEGs formed 2 separate clusters based on their positive or negative correlation with a measured parameter change (Fig S39). Correlation with body fat involved the highest proportion of DEGs in most tissues (Fig. 8b), except for BAT and HIPPOC.

To identify key regulators of the DEGs that were associated with phenotypic changes, we analyzed motif enrichment at the corresponding promoter DEGaPs. In SKM-GSN, the promoter DEGaPs of genes that were positively correlated with changes in VO2 max were significantly enriched for NR5A2 and ERRG motifs (Fig. 8c). The VO2 max-correlated DEG targets of those two TFs and of AP-2GAMMA and AP-2ALPHA, whose motif enrichment did not reach significance, showed considerable overlap (Fig 8d). The positive correlation between three of the DEGs targeted by either NR5A2, ERRG, or/and AP-2GAMMA (Me3, Rora, and Lgi3) and VO2 max change is depicted in Fig. 8e-g.

In SKM-GN, we also found a positive correlation between DEGs and body weight changes and identified the MyoD motif as significantly enriched at the corresponding promoter DEGaPs (Fig. 8h). While an enrichment for the COUP-TFII and ERRA motifs was also detected, it did not reach significance (Fig 8h). The Chd7, Rnf13, and RGD1562029 DEGs showed concurrent enrichment for the MyoD, COUP-TFII, and ERRA motifs (Fig 8i). The positive correlation between Chd7 expression and body weight change is illustrated in Fig. 8j.

As SP2, MAZ and BMYB have been identified as likely major regulators of the DARs-associated DEGs (DAR-DEG pairs; see Fig. 4), we asked whether those DEGs correlated with phenotypic changes. Oas2 expression in lung, a target of both SP2 and MAZ, was negatively correlated with body fat change (Fig. 8k), as was Sall2 in SKM-GN, another identified MAZ target (Fig. 8l). Three BMYB targets in SKM-GN were significantly correlated with VO2 max change, including Cited4 (Fig. 8m) and Pfkfb3 (Fig. 8n) positively correlated, and Sik1 (Fig. 8o) negatively correlated. In the significant body fat and body weight correlations, we noted sex differences both in EET-dependent RNA levels and in phenotypic changes, as week 8-trained males deviated more from controls than week 8-trained females (Fig. 8j-l). In view of the characterization of TF expression responses to EET, we next evaluated correlations between TF expression changes, at either the RNA or protein level (see Fig. 5) and phenotypic changes. We observed the strongest correlations between TF RNAs and body fat changes (Fig S40a), with Pknox1 in SKM-GN being the most positively correlated, and Srebf1 in the kidney being the most negatively correlated. Srebf1 in the kidney was also the most positively correlated with body lean mass changes. Rora was the most positively correlated with VO2 max changes. Zeb1 and Zeb2 shared similar correlation patterns, namely a positive correlation with body fat in WAT-SC. Gata2 had diverging correlation patterns depending on the tissue; it was positively correlated with body weight and body fat in WAT-SC, and negatively correlated with both measures in the lung. Similar to the correlations with TF RNA level-changes, the strongest correlations at the protein level were with body fat changes. Remarkably, liver showed the highest frequency of strong correlations (Fig S40b). In WAT-SC, JUN was the most positively correlated with body fat, while TCFL2 and EHF were the most negatively correlated. Similar to their respective transcripts, ZEB1 and ZEB2 proteins in WAT-SC both exhibited a positive correlation with body fat, while ZEB2 in lung was also positively correlated with body fat. IRF3 was the most positively correlated with body lean changes in SKM-GN, and CLOCK was negatively correlated with body fat in WAT-SC. Connecting EET-modulated genes and proteins, and specifically TFs, with correlated phenotypic changes suggest potential mechanistic roles. In particular, MAZ and SMAD3 target gene correlation with body fat changes reinforces their proposed functional role in the training response.

Methods
Animal procedures necessary for generating the data used by this paper from the MoTrPAC PASS1B study5 were approved by the University of Iowa’s Institutional Animal Care and Use Committee.

Inclusion and Ethics Statement
Research responsibilities and authorship were agreed amongst the collaborators during the research and writing processes. As this was a rat-based study of exercise, there were no international local collaborators to include in this study; however, we endeavored to include all participating researchers in the author list, either individually named or as part of the MoTrPAC consortium. Animal welfare regulations were strictly followed and animal procedures were approved by the institutional animal care and use committee at the University of Iowa.

Animal study design

Animal care
Male and female Fischer 344 (F344) inbred rats from the National Institute on Aging were transported to the University of Iowa a minimum of 4 weeks prior to starting exercise training. Rats were housed with the same sex, 2 per cage (146.4 square inches of floor space) in ventilated racks (Thoren Maxi-Miser IVC Caging System) on Tekland 7093 Shredded Aspen bedding and fed the Lab Diet 5L79 pelleted diet.

Rats were acclimated to a reverse dark-light cycle with lights off at 9:00am and lights on at 9:00pm, with temperature controlled at 68-77°F and humidity of 25-55%. For these studies we will use Zeitgeber Time (ZT) to refer to time of day relative to the time of lights off/lights on with lights off at ZT12. All experimental interventions and husbandry were performed during the active, dark phase of the rats under red light. All animal procedures were approved by the Institutional Animal Care and Use Committee at the University of Iowa.

Treadmill familiarization
Prior to exercise training, rats were acclimated to treadmill exercise on a Panlab 5-lane rat treadmill (Harvard Instruments, Model LE8710RTS). Day 1-2 consisted of static treadmill exploration for 10 minutes. Days 3-5 consisted of running at 6m/min at 0° incline for 10 minutes, speed was increased to 10m/min between Days 6-12. On Day 12, rat running behavior was scored for compliance with running at 12m/min for 5 minutes at a 10° incline. Based on running behavior, rats received a score from 1-4 with 4 being the highest score. Rats that were assigned a score of 1 were removed from the study. 25 male and 25 female compliant rats were randomized to a control or one of 4 training groups, resulting in 5 male and 5 female rats per group.

Progressive exercise training protocol
At 6-months of age, 1, 2, 4 or 8 weeks of exercise training began. Control rats were placed on a static treadmill for 15 min per day, 5 times per week. Exercise training consisted of a progressive training protocol 5 times per week at ZT13-20, to target 70% VO2max (see below). Training times each day for each rat subject are included in Supplemental Data S4. Week 1 sessions started at 13m/min for males and 16m/min for females at 5° for 20 minutes, with duration increased by one minute each day until reaching 50 min on day 31 of training. The treadmill grade was increased from 5° to 10° at the start of week 3. The treadmill speed increased at the start of week 2 (15m/min males, 18m/min females), 4 (18m/min males, 21m/min females), 5 (20m/min males, 23m/min females), 6 (23m/min males, 26m/min females), and 7 (25m/min males, 28m/min females) and was fixed for the final 10 days of training. Rats performing less than 4 days of training per week were removed from the study and euthanized as described below.

Body composition measurements
Body composition (lean tissue%, fat tissue% and body fluid) was measured for all rats 13 days prior to the start of training and 5 days prior to euthanasia in the 4 and 8-week training groups using the minispec LF90II Body Composition Rat and Mice Analyzer (Bruker, 6.2 MHz Time Domain Nuclear Magnetic Resonance (TD-NMR) system). VO2max was determined prior to commencing training in all rats, and during the last week of training for rats in the 4- and 8-week exercise groups, in a single-lane enclosed treadmill (Columbus Instruments Metabolic Modular Treadmill), with rats acclimated two days prior to testing. For testing, the rat was placed in the treadmill and testing began once oxygen consumption stabilized. The testing protocol consisted of a 15 min warm up at 9 m/min and 0° incline. The incline was increased to 10° and treadmill speed was increased by 1.8 m/min every 2 minutes72. During the test, electric shocks were used if the rat stopped running and sat on the shock area. Testing stopped when the rat sat on the shock area 3 consecutive times and did not respond to increased shock. Blood was then taken from the tail to measure lactate. VO2max was determined as a leveling off of oxygen uptake, despite increased workload, a respiratory exchange ratio above 1.05, and an unhaemolyzed blood lactate concentration ≥6 mM.

Euthanasia and tissue collection
Rats in the 1 week training group were euthanized and had tissues collected for omic analysis after 1 week of training, while rats in the 2 week training group were euthanized and had tissues collected for omic analysis after 2 weeks of training. A similar pattern is found for the 4 week and 8 week training groups. The 8 week control group of rats were euthanized and had tissues collected after eight weeks, similar to the 8 week training group. On the day of euthanasia, food was removed at ZT11.5, 3 hours before tissue collections which took place between ZT14.5-17.5, a minimum of 48 hours post their last exercise bout. Rats were deeply anesthetized with approximately 1-2% isoflurane in oxygen, and gastrocnemius, white adipose, liver, lung, and heart were removed under anesthesia. Following removal of the heart, the rat was decapitated using a guillotine. The brain was removed from the skull and hippocampus dissected. The remaining tissues ( kidney, brown adipose, and hippocampus) were dissected post death. All tissues were cleaned of excess connective/fat tissue and immediately flash-frozen in liquid nitrogen, placed in cryovials and stored at -80°C. Rat tissues were archived and cryopulverized at the MoTrPAC Biospecimens Repository, until distributed to Chemical Analysis Sites for respective assays5. RNA-seq, ATAC-seq and RRBS analyses were conducted in all eight tissues and time points. Proteomics and phosphoproteomics were conducted in all time points for six of the eight tissues: SKM-GN, heart, kidney, liver, lung, and WAT-SC.

Data production and quantification
ATAC-seq data generation
Nuclei from 30 mg white adipose, 15 mg brown adipose, and 10 mg of other tissue samples were extracted using the Omni-ATAC protocol with modifications73. The white adipose, brown adipose, and hippocampus tissues were processed using no-douncing nuclei extraction. Cryopulverized tissue powder was incubated in the homogenization buffer for 10 min at 4°C, tubes inverted every 2-3 minutes. Other tissue powder was incubated in the homogenization buffer for 5 minutes on ice and dounced 10 times using pestle A and 20 times with pestle B. Nuclei were stained with DAPI and counted using an automated cell counter. 50,000 nuclei (or max. 500 µl nuclei) were added to 1 ml ATAC-RSB buffer and spun at 1000 g for 10 minutes, and the supernatant was removed.

The nuclei pellet was resuspended in 50 µl of transposition mixture and incubated at 37°C for 30 minutes with 1000 rpm shaking. The transposed DNA was purified using Qiagen MinElute Purification kits (Qiagen # 28006), and amplified using NEBnext High-Fidelity 2x PCR Master Mix (NEB, M0541L) and custom indexed primers73. 1.8x SPRIselect beads were used to clean the PCR reaction and remove primer dimers. The ATAC-seq libraries were sequenced on a NovaSeq 6000 using 2x50bp with 35 million pairs of reads per sample.

ATAC-seq data processing and normalization
Reads were demultiplexed with bcl2fastq2 (v2.20.0) and processed with the ENCODE ATAC-seq pipeline (v1.7.0) (https://github.com/ENCODE-DCC/atac-seq-pipeline)74. Samples from the same sex and training group were analyzed as biological replicates. Reads were aligned to genome rn6.0.9675 with Bowtie 2 (v2.3.4.3)76. Duplicate and mitochondrial chromosome mapped reads were removed. Peaks were called using MACS2 (v2.2.4)77, both from reads from each sample and pooled reads from all biological replicates. Pooled peaks were compared with the peaks called for each replicate individually using irreproducible discovery rate78 and thresholded to generate an optimal set of peaks. Optimal peaks from all workflows were concatenated, trimmed to 200 base pairs around the summit, sorted and merged with bedtools (v2.29.0)79 to generate a master peak list. This peak list was intersected with the filtered alignments from each sample using bedtools coverage to generate a peak by sample matrix of raw counts. Peaks from non-autosomal chromosomes were removed. Peaks that did not have at least 10 read counts in four samples in a tissue were removed. Filtered raw counts were then quantile-normalized with limma-voom80. This version of the normalized data was used for downstream analyses.

ATAC-seq peak assignment to genomic features
Accessible regions identified using ATAC-seq were assigned to one of the nine terms of genomic features using Ensembl Rn6 GTF (gene annotation file) and function annotatePeak of package ChIPseeker81 (v1.8.6).

Nine genomic features are defined as:
Promoter (<=1 kb) (proximal promoter): within +/- 1 kb from the transcription start site (TSS);
Promoter (1-2 kb): +/- 1 to 2 kb from the TSS;
Upstream (<5kb) 2-5 kb upstream from the TSS;
Downstream (< 5 kb): within 5kb downstream of the transcription end site (TES);
5' UTR (5’ untranslated region); Exon; Intron; 3' UTR (3’ untranslated region);
Distal Intergenic: regions > 5kb downstream of TES or > 5kb upstream from next TSS;
Overlaps Gene: overlaps with gene annotation, but not in any terms above.

All ATAC-seq identified accessible regions were assigned to the closest genomic feature of a genome. Differentially expressed gene associated peaks (DEGaPs) are defined as all accessible regions assigned to the differentially expressed gene.

DNA methylation data generation
Rat tissues were disrupted in GenFind v2 lysis buffer (Beckman Coulter, Indianapolis, IN) with a tissue bead ruptor (Omni International, Kennesaw, GA). Genomic DNA was extracted in a Biomek FxP automation workstation (Beckman Coulter, Chaska, MN). DNA sample quantification was done by Qubit HS dsDNA assay (Thermo Fisher Scientific) and quality was determined by Nanodrop A260/280 and A260/230 ratios. Reduced representation bisulfite sequencing (RRBS)82 libraries were generated with the Ovation® RRBS Methyl-Seq kit from Tecan Genomics (Baldwin Park, CA). Quantity of the libraries was determined by the Qubit HS ds DNA assay (Thermo Fisher Scientific) and quality evaluation was done by either Bioanalyzer High Sensitivity DNA Chip or Fragment Analyzer High Sensitivity NGS Assay (Agilent Technologies, Santa Clara, CA). Libraries were sequenced on a NovaSeq 6000 platform (Illumina, San Diego, CA) using paired-end 100 base-pair run configuration.

RRBS data processing and normalization
bcl2fastq (version 2.20) was used to demultiplex reads with options --use-bases-mask Y*,I8Y*,I*,Y* --mask-short-adapter-reads 0 --minimum-trimmed-read-length 0 (Illumina, San Diego, CA). FastQC (v0.11.8) was used to calculate pre-alignment QC metrics83 and reads were indexed and aligned to the Ensembl Rattus norvegicus (rn6) genome using Bismark (v0.20.0)84. Bowtie 2 (v2.3.4.3) was used to quantify the percent of reads that mapped to globin, rRNA and phix sequence contaminants and spike-ins76. Chromosome mapping percentages were calculated with SAMtools (v1.3.1)85 and NuGEN’s “nodup.py” script quantified PCR duplicates.

CpG sites were selected for downstream analysis if they exhibited methylation coverage >= 10 in all samples. Individual CpG sites were divided into 500 base-pair windows and clustered with the Markov clustering algorithm R package MCL86. Quantile normalization of sites/clusters was conducted separately on each tissue using R package preprocessCore87.

RNA-seq data generation
Rat tissue powders were further disrupted using Agencourt RNAdvance tissue lysis buffer (Beckman Coulter, Brea, CA) with a tissue bead ruptor (Omni International, Kennesaw, GA, # 19-040E). The total RNA was quantified using NanoDrop (ThermoFisher Scientific, # ND-ONE-W) and Ribogreen assay (ThermoFisher Scientific). Total RNA quality was determined by Fragment Analyzer Standard Sensitivity RNA analysis (Agilent Technologies, Santa Clara, CA).

200-300 ng total RNA was used for library generation. The Universal Plus mRNA-Seq kit from NuGEN/Tecan (# 9133) was used to select polyadenylated RNA. The generated sequencing libraries contain dual barcodes (i7 and i5) and UMIs (unique molecular identifiers) to accurately quantify the transcript levels. The RNA-seq libraries were sequenced on a NovaSeq 6000 using 2x100 bp with 35 million pairs of reads per sample.

RNA-seq data processing and normalization
Reads were demultiplexed with bcl2fastq2 (v2.20.0). Adapters were trimmed with cutadapt (v1.18). STAR (v2.7.0d) was used to index and align reads to genome rn6.0.96 and gene annotations75. Bowtie 2 (v2.3.4.3) was used to index and align reads to globin, rRNA, and phix sequences in order to quantify the percent of reads that mapped to these contaminants and spike-ins76. UMIs were used to accurately quantify PCR duplicates with NuGEN’s “nudup.py” script (https://github.com/tecangenomics/nudup). QC metrics from every stage of the quantification pipeline were compiled, in part with multiQC (v1.6)88. Lowly expressed genes (having 0.5 or fewer counts per million in all but one sample) were removed and normalization was performed separately in each tissue. These filtered raw counts were used as input for differential analysis with DESeq289. To generate normalized sample-level data, filtered gene counts were TMM-normalized using edgeR::calcNormFactors, followed by conversion to log counts per million with edgeR::cpm90. The same normalization technique was used on the 8 week control samples of each tissue for cross-tissue comparisons.

Proteomics data generation
Liquid chromatography tandem mass spectrometry (LC-MS/MS) was conducted on six tissues: heart and liver at the Broad Institute and skeletal muscle, kidney, lung, and white adipose at Pacific Northwest National Laboratory (PNNL). Sample processing followed a modified version of a previous protocol91. Peptides were labeled using tandem mass tag (TMT)92 and samples were grouped into sex- and training time point-based TMT11 multiplexes. Multiplex samples were fractionated by high pH reversed phase separation. Heart and liver samples underwent online separation with a nanoflow Proxeon EASY-nLC 1200 UHPLC system (Thermo Fisher Scientific), and then analyzed with a Q-Exactive Plus mass spectrometer (Thermo Fisher Scientific). The remaining tissues’ samples underwent online separation with a nanoAcquity M-Class UHPELC system (Waters), and analyzed with a Q Exactive HF mass spectrometer (Thermo Fisher Scientific).

Phosphoproteomics data generation
Phosphopeptide enrichment was performed through immobilized metal affinity chromatography (IMAC)5. Phosphopeptides were eluted off IMAC beads in 3x70 µl agarose bead elution buffer, desalted with C18 stage tips, eluted with 50% ACN, and then lyophilized. Samples were reconstituted in 3% ACN / 0.1% FA for LC-MS/MS analysis. Heart and liver samples were separated by a nanoflow Proxeon EASY-nLC 1200 UHPLC system (Thermo Fisher Scientific) and analyzed with a Q-Exactive HFX mass spectrometer (Thermo Fisher Scientific). SKM-GN, WAT-SC, kidney and lung samples were separated by a Dionex Ultimate 300I0 UHPLC direct-inject system (Thermo Fisher Scientific) then analyzed with a Q-Exactive HFX mass spectrometer.

Proteomics and phosphoproteomics data processing and normalization
For heart and liver, raw MS/MS data samples were processed by a Spectrum Mill (v.7.09.215) (Agilent Technologies). For the remaining tissues, sample processing was implemented by an in-house cloud-based proteomics pipeline executed in the Google Cloud Platform5. In all tissues, MS2 spectra were processed and searched against the rat RefSeq protein database (downloaded November 2018). Log 2 TMT ratios to the common reference were used as quantitative values for all proteins and phosphosites. Principal component analysis and median protein abundance across samples were used to find sample outliers. Proteomics features that were not fully quantified in at least two plexes within a tissue and non-rat contaminants were removed. Median-centering and mean absolute deviation scaling of Log TMT ratios were done for sample normalization. Plex batch effects were removed using limma::removeBatchEffect function in R (v 3.48.0). Phosphoproteome data was not normalized to the total proteome due to the lack of complete overlap of phosphosites and total proteome features (80.5% - 89.7%).

Statistical analysis
Differential analysis
Differential analyses were performed in each tissue of each ome. Males and females in one dataset were analyzed separately. Limma with empirical Bayes variance shrinkage was used for ATAC-seq, proteomics, and phosphoproteomics data93; the edgeR pipeline for methylation analysis was used for RRBS data94; DESeq2 was used for RNA-Seq83. For all proteomics and ATAC-seq data, the input for differential analysis was normalized as described above. For RNA-Seq, the input was filtered raw counts, in accordance with the DESeq2 workflow.

To select analytes that changed over the training time course, we performed F-tests (limma, edgeR::glmQLFTest) or likelihood ratio tests (DESeq2::nbinomLRT, lrtest) to compare a full model with ome-specific technical covariates and training group as a factor variable (i.e. sedentary control, 1 week, 2 weeks, 4 weeks, 8 weeks) against a reduced model with only technical covariates. For each analyte, male- and female-specific p-values were combined using Fisher’s sum of logs meta-analysis to provide a single p-value, referred to as the training p-value. To account for false discovery rate across all statistical tests, the training p-values were adjusted across all datasets within each ome using Independent Hypothesis Weighting (IHW) with tissue as a covariate95. Training-differential features were selected at 10% IHW FDR.

We used the contrasts of each training time point versus the sex-matched sedentary controls to calculate time- and sex-specific effect sizes, their variance, and their p-values (e.g., using linear F-tests), referred to as the timewise summary statistics. Specifically, for limma models we used limma::contrasts.fit and limma::eBayes, for DESeq2 models we used DESeq2::DESeq, for edgeR models we used edgeR::glmQLFTest. Covariates were selected from assay-specific technical metrics that explained variance in the data and were not correlated with exercise training: RNA integrity number (RIN), median 5’-3’ bias, percent of reads mapping to globin, and percent of PCR duplicates as quantified with Unique Molecular Identifiers (UMIs) for RNA-Seq; fraction of reads in peaks and library preparation batch for ATAC-seq.

Genomic Distribution of Peaks and Methylation Sites
Gene and gene-region annotations for ATAC-seq peaks and methylation sites were generated through the annotatePeak function in the R package ChIPseeker81. Identified gene regions included proximal promoter (<= 1kb from TSS), distal promoter (1-2 kb from TSS), upstream (< 5kb), exon, intron, 5’ UTR, 3’ UTR, downstream (< 5kb), overlaps gene, and distal intergenic (outside the gene body or near vicinity). We then measured the frequency of peaks and methylation sites to be present in each gene region, as well as the frequency of each gene region annotation in DARs and DMRs. This was also conducted on the subsets of active peaks (we define peaks with a median normalized accessibility value > -1) and DEGaPs for each tissue. Significant changes in gene region distribution were identified by the fisher.test R function.

Tissue enrichment of differential analytes
For each tissue and ome (RNA-seq, ATAC-seq, RRBS), the differential analytes were collected. For each differential analyte in a given ome, we measured the mean baseline (8 week control samples) expression. Z-scores were generated measuring the differences in baseline expression across all 8 tissues for each differential analyte. Differential analytes were ordered and annotated by their tissue of differential behavior.

Identification of transcription factor motifs using HOMER
Transcription factor motif enrichment analysis was performed on sets of DARs, DMRs, and DEGaPs for each tissue. DARs and DMRs for motif enrichment analysis were selected for each tissue by satisfying an adjusted p-value threshold of 0.1. Similarly, DEGaPs for each tissue were selected by isolating the DEGs that satisfied an adjusted p-value threshold of 0.1, and selecting peaks annotated to the DEGs that contained a median normalized accessibility of -1. For genomic feature-specific analysis, DARs, DMRs, and DEGaPs were divided based upon their gene region annotation. The analysis was carried out by findMotifsGenome.pl (HOMER v4.11.1)96. It was performed on the ±50 bp flanking regions of the peak summits. The search lengths of the motifs were 8, 10, and 12 bp. We applied the -find flag to generate a list of all known rat motifs contained within the ±50 bp flanking regions of the summits for each peak in the ATAC-seq dataset and within the ±50 bp flanking regions of the center of each methylation site in the RRBS dataset, using the same settings as above. HOMER-identified TF motif enrichment for input lists of DARs, DMRs, and DEGaPs are calculated using a cumulative binomial distribution comparing TF motif representation among a target input list and a background list from which the target list is a subset.

DAR genomic feature TF motif enrichment analysis
Applying the output from HOMER, the top ten significantly enriched TF motifs among DARs, DMRs and DEGaPs in each tissue were selected for further downstream analysis and cross-tissue comparisons. TFs were removed from further analysis if their gene was not expressed in the tissue in which their motifs were enriched. TF motif enrichments for differentially accessible regions (DARs) divided into gene features were calculated using the Fisher test. The test compared the ratio of DARs containing the motif for a specific TF/non-DARs containing this motif in one genomic feature, and the ratio of DARs containing this motif / non-DARs containing this motif in other genomic features. p-values were adjusted and FDR cutoff = 0.1 to select significant motifs in specific genomic features. This analysis was repeated for DMRs and DEGaPs. For the top TFs enriched among the DARs, DMRs, and DEGaPs of each tissue, we calculate their baseline expression (in control samples) and generated z scores of baseline expression across all tissues. Pearson correlation was calculated for TF enrichment z score versus TF baseline expression across all tissues.

Correlations between DARs, DMRs, and DEGs
We selected DARs whose centers were within 500kb of a DEG TSS in each tissue. We then calculated the Pearson correlation of the L2FC of the DAR and the DEG for each EET sample relative to their sex control values using the cor.test function in R. We considered a DAR-DEG pair for further analysis if their training response Pearson correlation achieved a p-value less than 0.05 and the DAR contained a known motif for a TF expressed within the tissue. We identified TFs that were statistically enriched for containing binding sites within the DARs of DAR-DEG pairs using a binomial test comparing correlated targets:total correlated versus total targets:total active peaks in each tissue. We followed the same protocol to identify correlated DMR-DEG pairs and their associated enriched TFs.

TFs significantly responding to EET and enriched for DEG targets
TFs were identified that showed significant response to EET at the RNA-seq, proteomic or phosphoproteomic level. Among TFs with significant EET responses, we calculated their enrichment for the promoter regions of DEGs. Enrichment significance was determined with a binomial test comparing the number of TF promoter targets in DEGs : total TF promoter targets versus the number of DEGs : total gene count in a given tissue.

TF Enrichment among active promoter regions of DEGs with sex-consistent training responses
In each tissue, we identified subsets of DEGs that were up-regulated (L2FC > 0) in both sexes or down-regulated (L2FC < 0) in both sexes for each time point and then identified the DEGaPs associated with each set of up-regulated or down-regulated DEGs. We divided the DEGaPs by gene region, isolating the most predominant regions of distal intergenic, promoter (here a merge of the proximal and distal promoter assignments), and intron. For each TF expressed in a given tissue, we calculated the frequency of TF binding sites within each DEGaP population and then calculated z scores of binding site frequency across the DEGaP populations. As up-regulated and down-regulated promoter DEGaPs divided TF enrichment, we isolated the most enriched TFs for these two classes of DEGaPs, selecting the TFs most discerning between up-regulation and down-regulation within each tissue.

Cell-type deconvolution
Cell type deconvolution was conducted by the R package CellCODE97 using the getallSPVs function. Marker sets were generated using the IRIS (Immune Response In Silico98 and DMAP (Differentiation Map) reference datasets99. The Kruskal-wallis test was implemented with the R function kruskal.test to determine if the variability in cell type proportion across samples in a given tissue would suggest a significant training response or sex difference.

Pathway enrichment analyses
Gene set enrichment analysis (GSEA) was performed using ssGSEA2.0 with differential enrichment analysis t-scores (trained versus control) used as input. The gene set database used included canonical pathways available through the Human Matrisome database available through MatrisomeDB (http://matrisome.org/), and the MitoPathways database available through MitoCarta 3.082–84. Human gene symbols were mapped to rat orthologs before running the analysis. The ssGSEA2 function was run using parameters that avoid normalization, require at least 5 overlapping features with the gene set, and use the area under the curve as the enrichment metric (sample.norm.type = “none”, weight=0.75, correl.type = “rank”, statistic =“area.under.RES”, output.score.type = “NES”, min.overlap=5). (cite the landscape paper )

Overrepresentation analysis of significant DEGs was performed for each tissue separately using the R package gprofiler2:gost using Gene Ontology Biological Process, Reactome, WikiPathways and KEGG databases against background gene sets defined as a set of all genes expressed in each tissue, taken from the rna-seq data100. Top 10 non-redundant pathway enrichments based on Benjamini-Hochberg (BH) corrected false discovery rate for each tissue are displayed as bubble plots with sizes indicating the number of significant genes relative to the pathway size (number of genes in that pathway) and colors indicating the significance (BH-corrected p-value). At least 10 genes were required to be differential in a pathway with a maximum of 200 genes.

Correlations between Phenotypic Measures and DEGs
Phenotypic measures were calculated at weeks 4 and 8 of EET and in week 8 controls. Measures were presented as changes between time point and original baseline measurements in each rat. For each phenotypic measure and DEG combination, we calculated the Pearson correlation between the change in phenotypic measure between baseline and a given time point, and the gene expression of the DEG at the time point for each animal subject. We isolated the DEGs that exhibited > 0.5 or < -0.5 correlations with each phenotypic measure in each tissue and selected the DEGaPs annotated to the promoter region of the DEGs. TF motif enrichment significance among a set of positively or negatively correlated DEG’s promoter DEGaPs in a tissue were determined by an exact binomial test comparing the frequency of enrichment among phenotype-correlated DEGs versus general enrichment among the promoter DEGaPs in the tissue.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this work, some of the rat datasets generated by the MoTrPCA consortium are analyzed to model the regulatory programs\nduring endurance exercise training. Specifically, the authors focus on the role played by transcription factors (TFs) in gene\nexpression regulation and DNA methylation. This study is complementary (or a continuation) of other works on the same\ndatasets as:\n- MoTrPAC Study Group et al. Temporal dynamics of the multi-omic response to endurance exercise training across tissues.\nPreprint at https://doi.org/10.1101/2022.09.21.508770 (2022).\n- Nair, V. D., et al. Molecular adaptations in response to exercise training are associated with tissue-specific transcriptomic\nand epigenomic signatures. Cell Genomics. In press.\nThe amount and variety of data obtained by this consortium are amazing: different omic modalities, eight tissues, etc.\nHowever, the datasets analyzed in this work present two potential limitations: (i) Diagonal multi-omic dataset, since it seems\nthat the different omics have not been measured on the same rats or samples. (ii) Sample size.\nRegarding (i), it is not clear from the data description that samples for RNA-seq, ATAC-seq and RRBS were taken from the\nsame rats. I assumed that from the subsequent multi-omic analysis. Measuring all the omics on the same sample/rat\nprovides a much better scenario for integrative analysis of the multiple omics considered since we can correlate rat\nmeasurements across omics to infer regulatory relationships. On the contrary, the authors mostly correlate average log2FC\nmeasurements across omics, which hampers the identification of specific regulations per sex or training period. As the\nexperimental design is not optimal in this sense, it is not possible to apply the different methodologies designed for vertical\n(statistical) integration that provide multi-layer regulatory networks.\nRegarding (ii), the sample size of an experiment is a critical parameter to obtain robust and trustable results from the\nstatistical analysis. The authors should indicate from the very beginning what this sample size is since it is difficult to deduce\nfrom their dataset description. Again, I assumed they have 5 rats at most for each time point and sex. In addition, the signal-\nto-noise ratio may differ between omics and consequently, each omic may have different statistical power while having the\nsame sample size. The authors should the power of each dataset and it is possible to estimate power from the multi-omic\nperspective (see, for instance, https://www.nature.com/articles/s41467-020-16937-8).\nConsidering these two limitations, the authors should prove the validity of their results by experimental validation since the\nmain contribution of their work is not the dataset (published somewhere else) or the multi-omic integration strategy (which is\nquite limited by the experimental design and mostly based on correlations from few observations), but the regulatory patterns\nor relationships found, which are quite general.\nSpecific comments\nIn general, the description of data or statistical and bioinformatics analyses is not detailed enough to understand what the\nauthors did. Moreover, it was not easy to find this information in the links provided to access the scripts. For the sake of\nreproducibility, the manuscript would benefit from more clarity in the explanations.\nDifferential analysis: According to the authors, for each tissue and ome, they performed a separate analysis for each sex and\nthen combined the p-values. I did not understand why they did not include sex, training time and the interaction between\nthem in the same model to assess the effects better while keeping a higher sample size. Since they presumably have a\nbalanced experimental design, this should be the most optimal way to proceed.\nFigure 2: Please explain how these Z-scores were computed.\nDAR-DEG associations: The authors associated each DAR with their closest gene and then only considered DEGs. Do they\nassociate a DAR to a DEG if the differential change was found in the same week or for the same sex in both omics or simple\nto DEG in any comparison? If the latter is true, does this make sense?\nProteomics and phosphoproteomics: Adding these two omes to the picture seems promising since it can corroborate if a\ngiven TF found from RNA data is confirmed by the abundance of its corresponding protein and whether it is active according\nto its phosphorylation. However, the authors analyse these three omes separately and find different altered TFs from each of\nthem. Is it possible to somehow connect all this information in a more confirmative manner by trying to correlate the values\nfor the same TF across the three omics? It would also be interesting to correlate the TF values with their target genes to\ncorroborate potential regulations found by TF binding sites enrichment? By the way, SIX1 is discussed in this section, but I\ncould not find it in Fig. 5A, B or C.\nPathway enrichment: The authors should provide more details about this analysis, in which cases they apply an over-\nrepresentation analysis or a Gene Set Enrichment Analysis (GSEA). For GSEA, which is the score used to “rank” the omic\nfeatures?", "reviewer_2": "This study presents ATAC, DNAm, RNA-seq in 8 rat tissues during 8 weeks endurance training. This paper is a companion\nto these 2 articles:\n• MoTrPAC Study Group et al. Temporal dynamics of the multi-omic response to endurance exercise training across tissues.\n2022.09.21.508770 Preprint at https://doi.org/10.1101/2022.09.21.508770 (2022).\n• Nair, V. D., et al. Molecular adaptations in response to exercise training are associated with tissue-specific transcriptomic\nand epigenomic signatures. Cell Genomics. In press.\nThe first paper presents the whole “resource” of the MoTrPAC consortium, involving multiomics, while this paper is\nspecifically focused on TF dynamics using the ATAC, DNAm and RNA-seq data.\nNovelty: the dataset is very comprehensive and valuable for the field. However, it has already been independently published\nby the authors in other studies. This manuscript is rather an in-depth exploration of an already-published dataset. This\nexploration is nonetheless extensive and wide-ranging, and the authors extract several meaningful observations from the\nomic data.\nMethodology: The manuscript seems methodologically robust and uses up-to-date pipelines. The authors also share all of\nthe code necessary to reproduce the analyses. The results are dense and, at times, a bit confusing to the reader.\nGeneral comment:\nGeneral heatmaps of the papers: the authors should consider removing the cell borders (which they color in gray) for\nheatmaps with a lot of cells, because when zooming out (such as in 100% size), the heatmaps are often seen completely\ngray because of the high number of borders. I think it is done for some heatmaps with lots of cells, but I think it should be\ndone for everyone.\nResults comments:\nSection “Characterization of epigenetic and transcriptional responses to endurance training”\n1. The authors should show general dimensionality reductions of omics data (by tissue or all-together, depending on tissue\neffect). Or if these are more described in the companion paper, please cite here.\n2. Which are the “control” groups? 0 week housing? Or multiple week housing without training?\n3. Fig 1C: it is conceptually a bit different to display the percentage of shared “detected” genes or “peaks” (because these are\ndifferentially expressed or accessible across tissues), than “Methyl sites” because the detection of CpG sites depends on the\nRRBS coverage, not on anything else. Or do the authors mean to say “methylated sites”?\n4. In Fig. S2-S3-S4: I understand the rationale of looking at log2FC distributions in the sense of comparing the “effect size” of\nexercise across different tissues. However, for comparing across different omics, the authors should check if logFC in\ngeneral is associated to variability in the omic measurements (e.g. baseline levels of variance could be higher across\nMETHYLATION sites than across RNA counts or ATAC peaks, and this could be related to seeing higher logFC changes at\nDMRs than at DARs, etc. )\n5. Figure S2, S3, S4 in general are a bit confusing. The authors talk about logFCs but which is the specific comparison\nbeing done? When they state “comparison across sexes and time points” do they mean: control vs training (any timepoint),\ntaking into account sex in the model)? Maybe, they could explore some interaction analyses to better talk about\n“consistency” across sexes or timepoints.\n6. Figure S5. Is it possible to cluster the columns of this heatmap? Maybe by discretizing the enrichments to -1 (significant\nnegative NES), 0 (non significant), 1 (significant positive) or something like that. It could be nice to see the results, though if\nthey are not informative the original plot as it is is also applicable.\n7. “Differential analysis was conducted in male and female samples separately, suggesting that DEGs in WAT-SC are\nindependent of sex-specific cell type composition changes.” I don’t understand this. Does this mean that, only for WAT-SC,\nthe models used are separately for sex? Or just for this part and this is “data not shown”?\n8. For genomic distributions of DARs, DMRs… don’t the authors think that it could be informative to separate gains vs losses\nof accessibility or methylation?\n9. Figure 2 is very interesting. It is a nice way of showing how tissue specificity of epigenomic patterns can determine or pre-\ncondition the response to exercise. This could be extended even more in results or discussion. Here, again, it maybe could\nbe of interest to explore (particularly for the DNA methylation) taking into account direction of change: if up- or down-DMRs\ntend to be hypo or hyper-methylated. Because here, if one sees that “Conversely, DMRs in a given tissue were more\nhypomethylated in the tissue with the significant training response”, one would tend to think that “probably, most exercise\nchanges are hypermethylation”, because low-methylation regions can only go up, and very high-methylated regions can\nonly go down.\nSection: “Identification of distal correlated epigenetic regulation despite few DARs and DMRs\nmapping to adjacent DEGs”\n1. It is very interesting in figure 3E that for WAT-SC, specifically, the distribution of distance between DAR and DMR is quite\ndifferent. They tend to be located further away from each other, right? Can the authors discuss this or provide hypotheses\nhere?\n2. Figures 3F-H are indeed very interesting and serve to show a lot of patterns. 3G is slightly more bimodal in positive\ncorrelations, there’s a bit of an enrichment of positive correlation between far-away DMRs and DEGs. For 3H, I find it a bit\nharder to see this trend in the main picture. Can there be some statistical support for this?\nSection: “DAR-DEG pairs are associated with distinct pathways in each tissue and MAZ and\nSMAD3 represent key regulatory TFs”\n1. Do the authors think that here, in figure S11, the results are not as strongly tissue specific as for the prior intersections\nbetween DERs, DARs, DEGs ? There are some common pathways. To see this better, could you order the bubble plot (S11)\nby pathways first appearing in BAT, then HEART, and so on? Also, the color scale on the plot is confusing. Yellow usually\ntends to be “more significant” than read, but here it is reversed. Or maybe color by log10(adj. p).\nIf the authors agree, it may be better to mention that there are here some patterns shared by tissues. So maybe tissue\nindependent responses are more at the pathway rather than the specific gene level\n2. The Fig. 4 part is a bit harder to understand for a non-expert. What is the represented correlation? I understand that the\nauthors first do enrichment and identify certain TF binding sites in the DARs of DAR-DEG pairs (and this could maybe go in\nmain figure). But then, what is the specific correlation between a region and a gene or TF in the heatmap? Please clarify this\nsection in general as it is a bit confusing. Also for the DMR part.\nSection: “Characterization of TF expression responses to EET”\n1. “Various subsets of TFs exhibited significant changes at the transcriptome (Fig 5a), proteome (Fig 5b), and\nphosphoproteome (Fig 5c) levels in each tissue.”, so is this done with a multiple testing control in all genes? Only in all Tfs?\n2. It is increasingly being reported that protein and transcriptomic changes may not be very correlated. Do the authors have\nany findings on this? E.g. comparisons between the trends in change of specific significant TFs in RNA vs protein or\nviceversa? Or do the authors find clear correlations?\nSection: “DARs vs. DMRs vs. DEGaPs show distinct TF motif enrichment patterns that differentially correlate with TF gene\nexpression”\n1. The authors look at DARs, DMRs and DEGaPs (DEG-associated peaks). “To address this,\nwe analyzed TF binding site enrichment at either DARs, DMRs, or DEGaPs relative to all open\nchromatin peaks in each tissue”. Comparing DMRs to all open chromatin peaks is not very intuitive. I understand using all\npeaks as background for DARs, but for DMRs, what other alternatives can there be? (e.g. permutation sampling similar\nregions from the genome, etc.) Especially because RRBS has a very particular bias in the distribution of the profiled CpG\nsites. This can be a major issue, the authors should look at this in more detail.\n2. It is very interesting that DEGaP enrichment is very different to DAR enrichment. I agree that this strategy by the authors is\neffective to study the different putative regulation mechanisms that they propose.\n3. It is a bit strange to cite Fig S21-23 before S18-20\nSection: “DARs, DMRs, and DEGaPs show cross-tissue motif enrichment conservation at\nspecific genomic regions”\n1. In Fig. 6d-f, what is the background for doing enrichment in genomic regions for DARs, DMRs, DEGaPs ? Again, this is a\nvery important factor to then observe enrichments. If possible, please consider validation strategies using other backgrounds\n(such as permutation testing)\nSection: “Identification of DEGs associated with EET-induced phenotypic changes”\n1. What about doing correlations with immune cell compositions?", "reviewer_3": "In the manuscript entitled ”Multi-omic identification of key transcriptional regulatory programs during endurance exercise\ntraining”, Smith and co-workers have performed impressive omics analyses in multiple organs in male and female inbred\nrats exposed to long-term treadmill running. A very large amount of exercise induced changes, associations and correlations\nare presented which will most probably have a significant impact on future research related to exercise physiology.\nGeneral comments:\nThis is an interesting and ambitious study which has the potential to have a significant impact on future studies related to the\neffects of physical exercise. The ms includes a large number results and it would be helpful if the biological importance of\nthese results in different the tissues were discussed in more general terms. For instance, results from both genders are\npresented, and it would be helpful if the biological importance of gender-related differences were discussed. Multiple organs\nare analyzed which is a strength of the study. However, organs are presented as separate tissues although there is reason\nto believe that there is significant exercise induced communication between organs which will affect regulation of gene\nexpression, i.e., organ-to-organ communication in response to the exercise. This needs to be analyzed and discussed in\nmore detail.\nSpecific comments:\nTo improve the readability of the ms it is suggested that the list of used abbreviations is made complete and presented in\nalphabetic order.\nThe effect of circadian rhythm is discussed in the ms and all experiments were performed during the active dark phase under\nred light in rats acclimatized to a reverse dark-light cycle. Please add the time of the day when rats were trained on the\ntreadmill in the Methods section.", "reviewer_4": "The authors performed multi OMIC muscle analyses in rats to uncover TFs mediating gene responses to 1- 8 weeks of\nendurance training across multiple tissues. This study is leveraging a large-scale consortium data. Overall, the study is\ninteresting and is utilising multiple rat tissues and holistic approach to exercise. However, I noticed some fundamental\nissues with the study design, and collective analyses. My comments are as follow:\n1. Title: Please add the word rats to the title.\n2. Abstract: Please, again, use the model you used-rat model in the abstract. There is a very big distinct and leap between\nrats and humans. Please add the number of rats and their sex to the abstract.\n3. Introduction:\n- “Our companion multi-tissue analysis of the molecular response dynamics during endurance exercise training found that\nthe majority of differentially regulated genes are tissue specific whereas a small proportion are shared across multiple\ntissues (5,6)”\nI fail to understand how is this study significantly different the references 5 and 6. This is very weird that the authors has very\nsimilar papers published/pre prints. Is this study only adding the RRBS and TF to the mix?\n- Please mention which tissues did you analyse in the introduction.\n4. Methods:\n- The methods section describing the training is very vague and not clear to me.\nI am confused with the training method that the authors used. In the first figure it is not clear in which time points did you\ncollect the samples from the animals. ‘At 6-months of age, 1, 2, 4 or 8 weeks of exercise training began’ what was the length\nof the exercise training? How did you match the control rats in terms of sample collection? Did you collect samples from\nsome of the animals at 1,2, 4 or 8? Which time points did the authors used for the OMIC analyses and was it standardised\nbetween tissues?\n- If you used time and sex in the statistical model, can you explain what it means- does it mean that your outcome (say TF for\nexample) will be the same if you apply 2 weeks of training and 8 weeks of training in rats? Would you expect the same\nDNAm sites and TF to up/down regulate?\n- There is no mention of the number of animals and number in each sex (I might have missed so apologies if I did).\n- Cell-type deconvolution. Was it used only on DEGs? Have you attempted cell type deconvolution on the DNAm? You have\ndescribed many layers of OMICS and I just wonder why is the assumption that only DEGs will be cell-type specific. I am\naware of the Teschendorff method for DNAm in humans (PMID: 35277705) but not aware of a method for rats- have you\nthought about that?\n- RRBS data. Can you clarify if this was a whole genome sequencing or targeted sequencing and if so, did you use the\nsame method for the all tissues you generated.\n5. Results:\n- The results section is generally well described.\n- Again, this is a key point that might mask the entire analyses. The length of training is confusing and if you have analysed\nresponses to training after 1 week or 2 weeks vs 8 weeks, I anticipate a very different response. Fundamentally, I do not\nthink you can claim 1 or 2 weeks training is sufficient to understand to elicit an endurance response that can uncover the\nendurance-related pathways, even if you accounted for the training length in your statistical analyses. “To understand the\nepigenetic and transcriptional response mechanisms elicited during eight weeks of EET, we analyzed ATAC-seq, RNA-seq\nand RRBS profiles generated in skeletal muscle (gastrocnemius; SKM-GN), heart, hippocampus (HIPPOC), kidney, liver,\nlung, brown adipose tissue (BAT), and subcutaneous white adipose tissue (WAT-SC) from rats subjected to 1, 2, 4, and 8\nweeks of training”"}, "author_response": "Response: Thank you for this comment. To clarify, we are not comparing different assays\ndone on different animals--the same samples from the same tissues from the same\nindividual rats were used for all omics assays. You raise a valid point about the log2FC\ncorrelation analysis that had previously been measured comparing overall L2FC values.\nWe updated all correlation analyses to reflect sample level comparisons instead and have\nrevised our interpretation and which important relationships to highlight, as appropriate.\nThese changes are most visible in Figures 3 and 4 in the revised manuscript.\nResponse: There were 5 male and 5 female rats for each time point with a total of 50\nanimals studied simultaneously for all omes in eight tissues. We now reiterate this in the\nresults and methods sections. The expected change in VO2max was the main determining\nfactor for selection of age, training duration and number of included rats. We now include\npower analysis for the five different omes (including RNAseq, ATACseq, RRBS, protein\nabundance and protein phosphorylation) using the multiPower package as Supplemental\nFigure S5 in all eight tissues. Considerable power is present for all omes, except for the\nphosphoproteome.\nResponse: The reviewer raises two concerns about our study. With regards to the sample\nsize, multiomic power analysis shows considerable power for the various omes. We believe\nthe ability to identify thousands of differential analytes following FDR correction further\nindicates that this study captures significant molecular responses to training. With regards\nto the other concern, the omics data were all generated from the same rats, and now\nfollowing the reviewer’s suggestion, we have ensured that all cross-ome correlation\nanalyses are conducted at the individual subject level to make optimal use of the depth of\nthis data for each studied rat. Furthermore, we have identified interactions in human ChIP-\nseq data between a number of enriched TFs and their target genes identified in this study,\nas shown below in Response Figure 1A-D. We believe that these extra analyses address the\nreviewers’ concerns that precipitated the request for experimental validation.\nResponse Figure 1A. SP2 predicted targets in the manuscript and ChIP-seq data of SP2\nbinding with target genes.\nResponse Figure 1B. BMYB predicted targets in the manuscript and ChIP-seq data of BMYB\nbinding with target genes.\nResponse Figure 1C. BMAL1 predicted targets in the manuscript and ChIP-seq data of\nBMAL1 binding with target genes.\nResponse Figure 1D. MAZ predicted targets in the manuscript and ChIP-seq data of MAZ\nbinding with target genes.\nWe emphasize that the goal of this manuscript is to leverage an unparalleled large scale\nNIH 21 institution consortium exercise mapping study which has, as its goal, to generate\nreference data and perform analyses on those data. There is no mechanism or resources\nfor performing additional studies beyond the large scale reference study analyzed in this\nmanuscript. The goal of the present study is to elucidate the gene regulatory circuitry of\nexercise training.\nWe note that the reviewer’s concerns about the study design would apply equally to the\nmain overview MoTrPAC consortium paper published as a cover article in Nature that\nreported the overall temporal dynamics of the responses across all omes and was based\non the same data and study design as the present manuscript. That paper went through a\nvery extensive review without a request for additional studies. Overall, we think the present\nstudy is a valuable contribution to our understanding of molecular underpinnings of the\nresponse to exercise training and it is powered at the extreme limits of what is possible for\nthe data that was generated in a landmark multiomic exercise study.\nResponse: The code used to conduct the analyses and to generate the figures in this\nmanuscript are publically available at:\nhttps://github.com/gsmith990306/MoTrPAC_PASS1B_Transcription_Factor_Paper.\nTo address the reviewer’s comment, we have expanded the analysis portion of the methods\nsection to provide greater detail and clarity for all methods used.\nResponse: While we recognize the utility of increasing sample size when conducting\ndifferential analysis, the consortium identified different variability distributions between\nmale and female rats. This prompted us to separate the sexes when performing the\ndifferential analysis. This differential analysis decision was made for the main overview\npaper for this study (https://doi.org/10.1101/2022.09.21.508770), including all tissues and\nomes, and we sought to align our analysis with that one for consistency and comparative\npurposes. The F test that was chosen for the differential analysis looks to find any\nsignificant patterns of response over the course of training within each of the sexes. In that\nway, it does incorporate training time.\nResponse: The goal of Figure 2 is to determine if differential analytes in a particular tissue\nare generally more highly enriched in that tissue than at baseline. To that end, we\nspecifically normalized the RNAseq, ATACseq and RRBS data among baseline samples\nacross the eight tissues, and calculated z scores on the resulting tissue X significant\nanalyte matrices. We have revised the methods and results to clarify further this analysis.\nResponse: Because an F test was used to identify differential analytes in each ome, we do\nnot specify time-point specific differential signals. That said, when we associate DARs or\nDMRs to DEGs, we look for any pair within 500kb such that they are both significant in their\ntraining response across time points and sexes, and are either positively or negatively\ncorrelated with one another. Thus, they may have similar or divergent patterns of response\nacross time points.\nResponse: We propose in our analysis that the activity of TFs can be regulated in different\nmanners, including at the transcriptome, proteome and phosphoproteome levels.\nFurthermore, while most TFs can be found in the RNAseq data, the proteomics data only\ncovers a limited subset, and even when TFs are discovered in the proteomics dataset, they\nare often low in abundance and it’s less likely for significance to be identified in response to\nexercise training. Because of that, even if a TF is regulated at an RNAseq level but not found\nto have a significant training response at the proteomic level, we still consider this a\npotentially important TF regulatory event to be considered in our analysis. To that end, we\nlook separately for changes in TF activity at each omic level and then connect these\nchanges to enrichments for target DEGs. That said, identifying correlated responses to\ntraining for TFs at multiple omes is a very interesting question and we have now conducted\nthis analysis. In the next six figures, we plot the correlation, and correlation significance,\nfor the training response of TFs at the RNA vs Protein level, RNA vs Phosphoprotein level and\nProtein vs Phosphoprotein level, in each tissue with proteomics data. In each plot, cross-\nome comparisons are labeled with a TF name if the adjusted p-value is less than 0.1.\nNames are colored in orange if the TF in one of the two omes compared is significantly\nresponding to training, and colored in red if the TF in both of the omes compared is\nsignificantly responding to training. Of the TFs with significantly correlated behavior across\nomes, most of them do not show a significant response to training. Among those that do,\nnone of them were enriched for motifs located within the promoters of DEGs. We have\nincluded this analysis in the results and as Supplemental Figure S21 in the manuscript.\nFigure S21A: Gastrocnemius TF Cross-Ome Correlation Plots\nFigure S21B: Heart TF Cross-Ome Correlation Plots\nFigure S21C: Kidney TF Cross-Ome Correlation Plots\nFigure S21D: Liver TF Cross-Ome Correlation Plots\nFigure S21E: Lung TF Cross-Ome Correlation Plots\nFigure S21F: White Adipose TF Cross-Ome Correlation Plots\nResponse: We chose to focus on identifying significant TFs whose binding sites were\nenriched in DEGs compared to all tissue-expressed genes. This was because the variety of\nexpected patterns we could see, including delays in target responses over the course of\ntraining, wouldn't be identified from correlation alone. That said, this is an interesting\nanalysis approach, and we have now investigated the TFs responding to exercise training at\nan RNA, protein or phosphoprotein level and their correlation with target genes, both DEGs\nand non-DEGs. Significance is calculated based on the difference in correlation between\nDEG and non-DEG targets for each TF. Some TFs have stronger correlations with their target\ngenes and it makes sense that DEGs among targets tend to have higher correlations than\nnon-DEGs (since the TFs are also significantly responding to training). Interesting stand-out\nTFs here include Nur77 in SKM-GN, JunD in HEART and WAT-SC, and PU.1 in LUNG. We\nhave included these figures and analysis in the results section and as main figure panel 5E\nand Supplemental Figure S26.\nFigure S26A: SKM-GN correlation distributions of significant TFs and their target genes\nFigure S26B: HEART correlation distributions of significant TFs and their target genes\nFigure S26C: KIDNEY correlation distributions of significant TFs and their target genes\nFigure S26D: LIVER correlation distributions of significant TFs and their target genes\nFigure 5E/S26E: LUNG correlation distributions of significant TFs and their target genes\nFigure S26F: WAT-SC correlation distributions of significant TFs and their target genes\nResponse: We appreciate the thorough investigation. In Figure 5A, we only show a subset\nof the TFs with significant training responses because of lack of space. SIX1 does\nsignificantly change, it just did not meet the cutoff for appearance in that figure panel. We\nincluded it in Fig 5D because while its training response significance wasn’t the highest\namong significant TFs, its binding motif was still enriched at the proximal promoter of\nDEGs. We have included SIX1 in Fig 5A (and any other cases of missing TFs) to avoid further\nconfusion.\nResponse: We thank the reviewer for pointing this out. We now provide additional details in\nthe methods section and the associated figure legends. The pathway enrichment analysis\nper sex and timepoint displayed in figure S9 was performed using GSEA. For GSEA,\ndifferential enrichment analysis t-scores (trained versus control) were used as input.\nFeature-level data was summarized into gene-level data using the maximum absolute t-\nscore. The gene set database used included canonical pathways available through the\nHuman Matrisome database available through MatrisomeDB (http://matrisome.org/),\nand\nthe MitoPathways database available through MitoCarta. The pathway enrichment analysis\nfor DEGs associated with DARs, displayed in Figure S17, was an overrepresentation\nanalysis using the R package gprofiler2:gost against Gene Ontology Biological Process,\nReactome, WikiPathways and KEGG databases against background gene sets defined as a\nset of all genes expressed in each tissue, taken from the rna-seq data.\nResponse: Thank you for the suggestion, we have removed the cell borders from the\nheatmaps in the paper.\nResponse: We have now generated PCA plots in each tissue for the RNAseq, ATACseq and\nRRBS data, included in Supplemental Figures S6-8.\nResponse: The control group is 8 weeks housing without training, the control rats were\nplaced on non-moving treadmills 5 days a week to match the handling stimulation with the\ntrained rats. We have clarified this in the beginning of the results section.\nResponse: We apologize for the lack of clarity. The term “Methyl sites” refers to methylated\nCpG sites. Our main point with this figure is to highlight the general tissue-specificity of\nsignificant analytes in each ome studied. The CpG sites selected for this analysis are those\nthat satisfy a base threshold for methylation within each tissue. We have clarified this in\nthe results.\nResponse: This is a very good point. We have now incorporated into each supplemental\nfigure S2-S4 L2FC plots for non-differential features in each tissue and each ome to\nmeasure the general variability present. Indeed, methylation sites show more variability\nthan ATAC peaks or RNA counts, though less than their respective differential analytes.\nResponse: The histograms in Figures S2,S3,S4, are control vs training (any timepoint),\nseparated by sex (so Female Week 2 vs Control and Male Week 8 vs Control would be two\nseparate comparisons included in the histogram). The attached heatmaps in Figures\nS2,S3,S4, however, treat the different Sex/Week comparison L2FC values as separate\ncolumns. Here we show the distinct patterns of response across training. In Figure 1D, we\nsummarize the consistency, or lack thereof, across sexes and timepoints for each ome and\ntissue.\nResponse:\nThe heatmap now located in Figure S9 was generated from a code package that cannot be\neasily altered. We think it’s also very important to arrange the columns by tissue, sex, then\ntime point for simplicity of cross-comparison. The pathways in the rows are arranged by\nsimilar enrichments patterns across tissues, sexes, and time points.\nResponse: The differential analysis was conducted separately for male and female\nsubjects for all tissues and omes. We specifically highlight WAT-SC here because of its\nunique sex-specific cell type composition changes, compared to BAT which exhibits\ntraining-specific cell type composition changes. Sex-specific cell type composition\nchanges are okay for our deeper DEG-DAR or DEG-DMR comparisons because differential\nanalytes are identified within a specific sex, so the F test wouldn’t just be picking up sex\ndifferences. This means our identification of DEGs, DARs and DMRs would be robust\nagainst sex-specific cell type composition changes and any connection between a\ncorrelated DEG-DAR pair wouldn’t just be because they are both present in a cell type that\nis changing in frequency in the tissue. Compare this to BAT, where the training is linked to\ncell type composition changes, and we cannot be as confident that a correlated DEG-DAR\npair doesn’t exist solely because they are both from the same cell type.\nResponse: This would be an interesting addition to the analysis; however, by the nature of\nour F test-based differential analysis, our differential analytes are declared for a general\nchange in response over the course of training relative to the week 8 control rats. This\nmeans that a DAR can represent a diverse array of responses, and only a limited subset\nwould represent an increase in both males and females at a given time point. To this end,\nwe looked at isolating specifically the DARs that have a L2FC > 0.1 in both males and\nfemales relative to control (labeled W8 Up-Regulated) or a L2FC < -0.1 in both males and\nfemales relative to control (labeled W8 Down-Regulated). These subsets were quite small\nproportions of our overall set of DARs.\nResponse Figure 2A: Distribution of DARs labeled as up-regulated or down-regulated at\nweek 8 (in both sexes) or neither in each tissue. Neither = DARs that either do not show any\nincrease (L2FC > 0.1) or decrease (L2FC < -0.1) in W8 or show increases or decreases that\ndisagree between males and females.\nWhen looking at peak distributions, the smaller the number of peaks, the more likely\nrandom chance can affect the overall distribution of peaks between different gene body\nregions, and so we’re less likely to get informative results by examining such small peak\nsets. That said, we did examine the distribution of gene regions in these smaller peak sets\nand we do see relatively similar results to the total DAR lists for each tissue:\nResponse Figure 2B: Distribution of DARs among gene regions. DARs are divided into up–\nregulated or down-regulated subpopulations for each tissue, as well as examined in total.\nSubsetting the DMRs for up-regulation and down-regulation at week 8 also produces much\nsmaller groups to measure gene region composition:\nResponse Figure 2C: Distribution of DMRs labeled as up-regulated or down-regulated at\nweek 8 (in both sexes) or neither in each tissue. Neither = DMRs that either do not show any\nincrease (L2FC > 0.1) or decrease (L2FC < -0.1) in W8 or show increases or decreases that\ndisagree between males and females.\nAnd these smaller W8 Up-Regulated and W8 Down-Regulated groups have the following\ngene region composition:\nResponse Figure 2D: Distribution of DMRs among gene regions. DMRs are divided into up–\nregulated or down-regulated subpopulations for each tissue, as well as examined in total.\nThe smaller samples sizes increase the variability in gene region frequencies, although the\nmajor trends are still visible overall.\nResponse: This is a very interesting idea and your instincts were quite right. As with Figure\n3 (in the response) above, we identified DEGs, DARs, and DMRs that are up-regulated in\nboth sexes at week 8 (or down-regulated in both sexes at week 8). Generally the patterns\nwe see in the total list of DEGs, DARs and DMRs are preserved (though with a little more\noverlap in expression than in the larger differential analyte sets); however, the down-\nregulated DMRs do behave a little differently. In some tissues, like gastrocnemius, they are\nstill hypomethylated relative to other tissues, but we see less of a strong pattern in other\ntissues, like brown adipose tissue. We have included this as Supplemental S13 figure in the\npaper.\nFigure S13A: Baseline tissue RNAseq expression comparison for Up-Regulated DEGs\nFigure S13B: Baseline tissue RNAseq expression comparison for Down-Regulated DEGs\nFigure S13C: Baseline tissue ATACseq accessibility comparison for Up-Regulated DARs\nFigure S13D: Baseline tissue ATACseq accessibility comparison for Down-Regulated DARs\nFigure S13E: Baseline tissue RRBS methylation comparison for Up-Regulated DMRs\nFigure S13F: Baseline tissue RRBS methylation comparison for Down-Regulated DMRs\nResponse: The main driver here is the lack of DARs in WAT-SC. We only identify 5 in the\ntissue so the odds of a DMR being close to a DAR are slim. We have pointed this out in the\nmanuscript more clearly.\nResponse: We agree that 3H does not highlight this pattern that does appear in 3G. The\ntissue-specific nature of these patterns are a little more evident in the Supplemental\nFigures S8-S10. The tissues driving this pattern in 3G are BAT and WAT-SC. The bimodality\nof BAT can be explained by the cell type composition changes in response to training. WAT-\nSC is a bit more interesting as we do not have the connection to training-induced changes\nin cell type composition. Because of the tissue-specific response here, we were less\ninclined to focus on this as a main narrative, despite the interesting correlation pattern. In\n3H, there are fewer neighboring DARs and DMRs to connect, so we believe the likelihood of\npicking up distal connections is a lot less likely. In fact, liver and lung are the primary\ndrivers for the pattern of close connections between DARs and DMRs. Because of the\nlimited numbers of neighboring DARs and DMRs, we would be hesitant to make any\nstatistical claims about the lack of pockets of distal correlation.\nResponse: There is some overlap between pathways identified for BAT and lung, but\ngenerally, the pathways identified are tissue specific again. We have clarified this better in\nthe results section, and rearranged the pathways by tissue enrichment to highlight this\ndistinction, as seen below.\nFigure S17: Re-clustered pathway enrichments\nResponse: The goal here is to identify TF binding motifs that are enriched among\ncorrelated DAR-DEG pairs, and to highlight the peaks containing their binding sites and the\nassociated correlated DEGs. Each heatmap represents the correlated peaks-genes that\ncontain the TF binding site. The DMR part is similar with regard to DMR-DEG pairs, though\nFigure 4H is an overview figure showing all TFs that show significant enrichment for DMR-\nDEG pairs among binding sites. We have redesigned this figure to be anchored in plots of TF\nenrichment significance among DAR-DEG pairs and DMR-DEG pairs, then identifying\ninteresting specific TFs and their targets to detail further.\nResponse: This is done with multiple testing control in all genes/proteins/phosphosites.\nResponse: There are differences between TFs that are changing at the RNA and protein\nlevel. In part, our proteomic data is not exhaustive, with a number of TFs not measured in\nthe proteomic analysis. Per Reviewer 1, we have conducted this analysis and, while we\nidentify select examples of correlated behavior cross-omes, the majority of significant\nresponses to training are not correlated cross-ome. We have included this as\nSupplementary Figure S21.\nResponse:\nWe apologize for the confusion and this is an error in the text after we added the RRBS data\nto our analyses. The DARs and DEGaPs are measured relative to all open chromatin peaks;\nhowever, the DMRs are measured relative to all methylated CpG sites. Thank you for\nidentifying this error and we have corrected this in the manuscript.\nResponse: We’re glad you appreciate this interesting finding.\nResponse: You are correct that the majority of Fig S18 and S19 are referenced after S21-\nS23. We have rearranged the supplemental figure panels in this section to flow better with\nthe text.\nResponse: The background here is, for DARs and DEGaPs, the TF enrichment among active\npeaks (peaks that satisfy a minimum level of accessibility in the tissue) or, for DMRs,\nmethylated CpG sites in the tissue.\nResponse: This is a good idea to incorporate into the earlier discussions of cell type\ncomposition. We have now generated heatmaps, included as Supplemental Figure S11, for\neach tissue, of the correlation of each DEG with changes in immune cell type composition.\nThe distinct patterns for brown and white adipose relative to the other tissues are quite\nstriking. We have also included a new supplemental table S1 detailing the correlation\npatterns of each tissue’s DEGs with the immune cell type composition changes. In the\nsection where we highlight correlations with phenotypic changes, we are seeking to identify\nregulatory relationships between TF and target gene that are connected to the phenotypic\nchanges in response to EET. Identifying the TFs enriched for DEGs tightly correlated with\ncell type composition changes would suggest that the TF is merely active in a cell type\nwhose behavior is changing, as opposed to changes in the regulatory behavior of the TF\nitself. This is still an interesting line of research, but for this paper, we have sought to focus\non EET-induced changes that are independent of cell type composition changes.\nResponse:\nWe have included an extra paragraph in the discussion on the tissue-specific nature of\nexercise training response in the study.\nResponse:\nWe have included an extra paragraph describing some of the key sex-specific changes\nidentified in the study. With the inclusion of the story of IRF8 acting on monocyte cell type\nproportion in lung in the response to exercise, we have highlighted a critical response to\ntraining that is also sex-specific in nature.\nResponse: Examining cross-tissue communication is a really compelling idea to\ninvestigate and an important part of the research goals for the MoTrPAC consortium across\nour different experimental studies. Tissue samples were collected from animal subjects in\nthis study 48 hours post exercise and with time periods so spread apart, it is unlikely we\ncan address this question directly, tracking the movement of molecules between tissues.\nFollow-up MoTrPAC studies on the response to an acute bout of exercise, in both rats and\nhumans, with multiple post-exercise time points in the first 48 hours post exercise will seek\nto address these questions directly. We can, however, document the shared responses to\nendurance exercise training by TFs across tissues. We found that the responses are largely\ntissue-specific; however, there are cases of cross-tissue shared response.\nFigure S20A: TF significant training responses in transcriptomics across tissues\nFigure S20B: TF significant training responses in proteomics across tissues\nFigure S20C: TF significant training responses in phosphoproteomics across tissues\nIf we quantify these overlaps, we identify in RNAseq, considerable overlap in significant TF\nresponses between SKM-GN and HEART, and between BAT and WAT-SC. Overlap is much\nmore limited in the proteomics data, while in the phosphoproteomics data, HEART has\nconsiderable overlap with LIVER and LUNG and LUNG also has considerable overlap with\nLIVER and WAT-SC. We have included these findings as supplemental figure S20 in the\nmanuscript.\nFigure S20D: Overlap of TF significant training responses in transcriptomics across tissues\nFigure S20E: Overlap of TF significant training responses in proteomics across tissues\nFigure S20F: Overlap of TF significant training responses in phosphoproteomics across\ntissues\nResponse: We made the abbreviation list more inclusive and alphabetic.\nResponse: We have included a supplemental table S4 with the time of day for each daily\ntraining bout for each animal.\nResponse: We have made the addition\nResponse: We have incorporated rat and sample number and the fact that both sexes are\nequally represented into the abstract\nResponse: The MoTrPAC rat exercise training study has generated a number of papers\nanalyzing this dataset as part of a publication package to be published concurrently, or\nclose to concurrently. A central landscape paper acts as an overview of the experiment and\nomic results, while individual companions take deeper dives into specific tissues, omes or\nanalysis questions. This paper adds a deep analysis of TFs and transcriptional regulation,\ngetting to the heart of the goals of the MoTrPAC study, seeking to identify molecular\ntransducers of exercise response.\nResponse: We have added the tissue names to the introduction.\nResponse: Rats began the study at 6 months of age and then underwent 1,2,4 or 8 weeks\nof exercise training. At the end of a given time point, rats were sacrificed and tissues\ncollected. Control rats underwent the 8 weeks of non-exercise to be the same age as the 8\nweek exercise training condition rats. Omic analysis was done on the same tissue samples\ncollected at each of the 1,2,4 and 8 week time points. We have clarified this in the\nmethods.\nResponse: We identify analytes across omes that are responding to training at different\npoints in the training process. Some would change 2 weeks post exercise, and some at 8\nweeks. One change would not necessarily suggest a change at the other. There is a subset\nof changes that are consistent across time points. Also, changes at one time point in a TF\nmight cause a change at a later time point in target genes. That is why when we look at TFs\nresponding to training, we then search for enrichment of DEGs among binding targets, even\nif the DEGs do not respond at the same time point.\nResponse: We have made this more visible in the beginning of the results section.\nCurrently the number of rats used (25 male and 25 female - 5 per control and exercise time\npoint group), is only mentioned in the middle of the methods.\nResponse: The cell type deconvolution used all genes in the dataset and compared\nagainst a reference single cell dataset from a previous study. We used RNAseq and not\nDNAm because we have a reference based on scRNAseq data. We show in this paper that,\nwhile the cell type deconvolution analysis was done using just the RNAseq data, we see\nthe effects of the cell type changes across ome. The best example is brown adipose, where\nwe see the increase in immune cell proportion post training. The strong correlation or anti-\ncorrelation across the majority of significant training response analytes in all omes\nregardless of their proximity to each other in the genome suggests that the same force of\ncell type proportion changes is impacting RNAseq, ATACseq and RRBS equally.\nResponse: This was whole genome sequencing and the same method was used for all 8\ntissues.\nResponse: We appreciate the kind comment.\nResponse: We incorporated analysis from each of the time points into our analysis. We\nwanted to track what was changing early in the training response and what took more time\nto develop. Figure 7, for example, focuses on TF enrichment to significant responses at 8\nweeks of training. Examining each time point also gives us more power to examine\ncorrelated changes in training response. We can be more confident in a shared pattern of\nexpression if it not only mirrors week 8, but also weeks 1,2,4."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed most of the reviewers' suggestions. I still have concerns about the robustness and validity of some of the results presented, not only due to the reduced sample size per time-point, but also to the low effect size of EET. This can be observed in the PCA plots in Figures S6-S8, where the main source of variability is the gender of rats. As the authors performed independent analysis for each sex, we can see in the PCA plots that, in general, in both sexes, control rats are mixed with EET rats in almost all tissues and omics. This hampers the task of finding differences between EET and control groups, as also shown in the power study (Figure S5), where the sample size required to achieve the minimum statistical power is between 100 and 1000 rats per group in all omics and tissues. However, I recognize the authors' effort to improve the manuscript and dive into the data to extract meaningful biological conclusions, while I fully understand the impossibility of performing experimental validation of their results.", "reviewer_2": "The authors have addressed most of the comments and concerns raised by the referees, and the quality of the paper has been substantially improved. I believe it will be a study of great interest for researchers working in the field.", "reviewer_4": "I thank the authors for incorporating my comments. I have no further comments/ inquiries."}, "author_response": "Response: We appreciate the concerns of the reviewer that the main driving variability in\nthe top principal components is frequently sex. This, and the discrepancy in variation\nstructure between male and female rats, led to the differential analysis being conducted\nseparately on male and female rats (with a consensus adjusted p-value based on\ncombining the results of the separate analyses). To assuage concerns about the presence\nof non-sex variability in our data, we include here plots of the association of top 10 PCs\nwith Sex and Group (exercise time points and control) for RNAseq, ATACseq and RRBS:\nResponse Figure 1: Association of RNAseq PCA results top 10 PCs for each tissue with Sex\nand Group\nResponse Figure 2: Association of ATACseq PCA results top 10 PCs for each tissue with Sex\nand Group\nResponse Figure 3: Association of RRBS PCA results top 10 PCs for each tissue with Sex\nand Group\nWhile sex is often the driving association of PC1, considerable connection is generally\nfound between other PCs and exercise group. This coupled with the thousands of\ndifferential analytes that we have identified in this study lends credence to the value at its\ncurrent sample size.\nResponse: We thank the reviewer for their thorough investigation of our research.\nResponse: We appreciate the reviewer’s thoughtful suggestions and are pleased that we\nhave satisfied their concerns.\nResponse: As with reviewer #2, we’re pleased that we have satisfied their concerns and\nappreciate their valuable suggestions and edits."}]

Decision letter