Abstract The central dogma of molecular biology dictates the general flow of molecular information from DNA that leads to a functional cellular outcome. In skeletal muscle fibers, the extent to which global myonuclear transcriptional alterations, accounting for epigenetic and post-transcriptional influences, contribute to an adaptive stress response is not clearly defined. In this investigation, we leveraged an integrated analysis of the myonucleus-specific DNA methylome and transcriptome, as well as myonuclear small RNA profiling to molecularly define the early phase of skeletal muscle fiber hypertrophy. The analysis of myonucleus-specific mature microRNA and other small RNA species provides new directions for exploring muscle adaptation and complemented the methylation and transcriptional information. Our integrated multi-omics interrogation revealed a coordinated myonuclear molecular landscape during muscle loading that coincides with an acute and rapid reduction of oxidative metabolism. This response may favor a biosynthesis-oriented metabolic program that supports rapid hypertrophic growth. Keywords: RNA sequencing, small RNA sequencing, RRBS, mitochondrial respiration, oxidative metabolism, epigenetics Graphical Abstract Graphical Abstract. [46]Graphical Abstract [47]Open in a new tab Introduction Skeletal muscle cells (myofibers) are syncytia containing hundreds of nuclei (myonuclei) in an expansive cytoplasmic volume.^[48]1,[49]2 A unique feature of myofibers is the remarkable capacity for growth in response to stimuli such as mechanical overload (MOV).^[50]3 Several mechanisms underlie the cellular and molecular responses to MOV that contribute to skeletal muscle hypertrophy, including increases in mammalian target of rapamycin complex 1 (mTORC1) signaling, ribosome biogenesis, and protein synthesis rates.^[51]3 In addition to these well-characterized processes, emerging evidence suggests a role for novel mechanisms, including epigenetic changes and alterations to microRNA (miRNA) expression.^[52]3 The role of miRNAs in gene repression through the RNA-induced silencing complex (RISC)-mediated targeting of transcripts in the cytoplasm is well established. However, mature miRNAs are also localized to the nucleus. Emerging work shows miRNA post-transcriptional regulatory functions within the nucleus.[53]^4–6 Several studies identified a cytoplasmic-nuclear transport mechanism for mature miRNAs^[54]7,[55]8 and direct miRNA targeting of nuclear-retained noncoding transcripts.[56]^9–11 Moreover, transcriptional gene silencing has also been reported, whereby miRNAs guide the RNA-induced transcriptional silencing (RITS) complex to promoter-associated transcripts, altering DNA methylation and chromatin structure.[57]^12–15 Specific miRNAs have since been found to play a role in the formation of epigenetic silencing complexes to induce hypermethylation-mediated transcriptional gene regulation in various cell types.[58]^16–18 However, the role of myonuclear miRNAs in skeletal muscle adaptation is not understood. We previously investigated the myonuclear DNA methylome and transcriptome, separately, in response to an acute hypertrophic stimulus in the plantaris muscles of mice.^[59]19,[60]20 Specifically, myonuclei were labeled with green fluorescent protein (GFP) before MOV and purified following 72 h of MOV by fluorescent-activated nuclei sorting (FANS). We then isolated DNA and RNA for myonuclear reduced representation bisulfite sequencing (RRBS) or myonuclear RNA sequencing (RNA-seq), respectively.^[61]19,[62]20 Due to the incompletely defined role of DNA methylation, it can be difficult to assign individual CpG sites to specific genes to infer transcriptional regulation. Also, simply comparing differentially expressed (DE) genes with hypo- and hypermethylated promoter regions typically underestimates the degree of overlap. In the current investigation, we integrated epigenetic alterations with myonuclear gene expression changes by using a novel multi-omics integration method adapted from chromatin immunoprecipitation sequencing (ChIP-seq) and RNA-seq analyses. The binding and expression target analysis (BETA) integration approach lends a holistic understanding of the CpG sites that are significantly contributing to gene expression changes.^[63]21 Furthermore, we used our conditional myonuclear labeling approach to perform small RNA-seq on purified myonuclei. Mature miRNAs were then computationally combined with the myonuclear mRNA transcriptome. Collectively, we link myofiber-specific epigenetic and post-transcriptional factors to reduced metabolism-related gene expression following 72 h of MOV; this included hypermethylation of downregulated mitochondrial respiratory chain complex genes. The downregulation of oxidative phosphorylation (OXPHOS) after acute MOV may be due to pyruvate shunting from the mitochondria in favor of aerobic glycolysis and biosynthetic precursor production for growth.^[64]22 Materials and Methods Animals All animal studies were performed in accordance with institutional guidelines and approved by the Institutional Animal Care and Use Committee of the University of Kentucky. Mice were housed in a temperature- and humidity-controlled room, maintained on a 14:10-h light-dark cycle, and food and water were provided ad libitum throughout the entire study. We used the myofiber-specific Tet-On system to label myonuclei as previously described for the HSA-rtTA^+/−;TRE-H2B-GFP^+/− (HSA-GFP) mice.^[65]23 Mice (female, ∼3 mo of age) were given doxycycline (0.5 mg/mL, 2% sucrose) for 5 d followed by a minimum 1-wk washout period prior to MOV (3 wk maximum). This strategy leads to the labeling of approximately 90%-95% of myonuclei, with negligible off-target labeling of nonmyonuclei.^[66]23 Mechanical Overload Following doxycycline treatment and a washout period, mice underwent sham surgery or synergist ablation, as previously described.^[67]19,[68]20,[69]24 Sham surgery involved the same steps of synergist ablation without removal of muscle. Following surgeries, mice resumed ambulatory cage behavior for 72 h. At this time, mice were euthanized in the morning by intraperitoneal injection of sodium pentobarbital, followed by cervical dislocation. Plantaris muscles were Dounce homogenized and subjected to FANS for purification of myonuclei according to our published protocols.^[70]19,[71]20,[72]25,[73]26 DNA and RNA from myonuclei were isolated as previously reported.^[74]19,[75]20 DNA Methylation Data Processing and Statistical Analysis We used our previously published myonuclear RRBS dataset for DNA methylation analysis.^[76]19 Quality control and adapter sequence trimming were performed using FastQC and Cutadapt, respectively, as parts of the Trim Galore wrapper. Low-quality base calls (Phred score <20) were removed prior to trimming adapter sequences. Bismark aligner was used to align the sequence reads to the bisulfite-converted GRCm39 genome prior to data processing. Coverage (.cov) files produced from Bismark aligner were used for data analysis in the methylKit R package.^[77]27 MethylKit was used to pool samples into their respective groups to maximize read coverage across the genome using a minimum read cutoff of 10 reads per base, and minimum base coverage of 1 sample per group. Differentially methylated regions (DMRs) were determined by genomic ranges for every gene promoter as defined by the mm39.bed file obtained from National Center for Biotechnology Information. Fisher’s exact test with sliding linear model correction for false discovery^[78]28 was used to qualify both differentially methylated sites and differentially methylated promoters within the dataset. Percent methylation and percent differential methylation were then obtained from methylKit following analysis. Myonuclear mRNA Sequencing Data Processing and Statistical Analysis We used our previously published myonuclear RNA-seq dataset for myonuclear transcriptome analysis.^[79]20 Quality control was performed using FastQC and Cutadapt was used to trim adapter sequences. STAR aligner was used to align sequenced reads to the GRCm39 genome. Binary alignment map (.bam) files resulting from STAR aligner were further processed using SAMtools. Bam files were indexed and unmapped reads and reads with low coverage were filtered prior to generating gene read counts. Limma^[80]29 was used to determine differential gene expression between sham and MOV plantaris myonuclear mRNA and raw P-values were adjusted for false discovery using the Benjamini–Hochberg false discovery rate (FDR) correction. Binding and Expression Target Analysis Integration for Combining Gene Expression With DNA Methylation BETA is a software that provides an integrated analysis of transcription factor binding to genomic DNA and transcript abundance using ChIP-seq and transcriptomics (RNA-seq) datasets.^[81]21 BETA takes into consideration the distance of the regulatory element relative to the transcription start site (TSS) by modeling the effect of regulation using a natural log function termed the regulatory potential (eqn [82]1), as described previously by Tang and colleagues.^[83]30 In this study, we adapted this software to integrate CpG methylomics with RNA-seq data instead. To do so, CpG islands (whether hypo- or hyper-methylated) were converted into “methylation peaks” similar to transcription factor binding peaks, which are built using the GRCm39 CpG island bedfile downloaded from the UCSC Genome Browser. Only genes differentially expressed with adjusted P < .05 from RNA sequencing analysis were included as input for gene expression. The BETA basic command was run with the following parameters: “-c 0.05—df 0.05—da 500.” graphic file with name TM0001.gif (1) In our study, a CpG island (k) within 100 kb of TSS of gene (g) is included in the calculation for the regulatory potential score (s). The distance between the CpG island and the TSS is expressed as a ratio relative to 100 kb ( Inline graphic ). The scoring is weighted based on the distance from the TSS (higher for smaller distances, lower for larger distances). BETA Integration Pathway Analysis Pathway analysis was performed using the enrichR R package. Up- and downtarget genes resulting from BETA integration of DNA methylation and RNA-seq datasets were included in independent up and down gene sets. The 2023 Gene Ontology (GO) database including biological processes, molecular functions, and cellular components was used to annotate up and down gene sets and determine enriched ontologies for each gene set. GOplot R package was used to combine log[2] fold change for each gene with their respective enriched pathways. Myonuclear Small RNA-seq The myonuclear RNA isolated from MOV plantaris muscles of mice that underwent synergist ablation (n = 4) or control plantaris muscles of mice that underwent sham surgery (n = 3) was analyzed by small RNA-seq. Library preparation was performed by Norgen Biotek (ON, Canada) using a Norgen Biotek Small RNA Library Prep Kit (Cat. 63600), and small RNA-seq was performed by Norgen Biotek using an Illumina NextSeq 500 platform. Data analysis was performed using the exceRpt small RNA-seq Pipeline (v4.6.2.) (Genboree Bioinformatics Research Laboratory). The databases miRbase version 22 for miRNAs, gtRNAdb for tRNAs, RNAdb, piRNAbank, & RNAcentral for piRNAs, and the reference genome GENCODE version 22 (mm10) were used for analyses. miRNA Target Enrichment and Network-Based Analysis Targets of differentially expressed (DE) miRNAs (adj. P < .05) were identified using miRTarBase Release 9.0.^[84]31 Differentially expressed genes from the myonuclear mRNA-seq analyses were compared to genes identified as targets of DE miRNAs, and pathway analysis was performed on overlapping genes using the EnrichR web-based tool.[85]^32–34 To identify potential miRNA targets and pathways that may have been missed by the miRTarBase database, an in-silico prediction of the miRNA target genes was performed. Target genes were detected using both miRanda and RNAhybrid tools. We used a custom python script to select shared miRNA: target gene interactions between the two software programs. Minimal free energy was set as ΔG°= −18 kJ. Each algorithm compares the 3’UTR gene region to the miRNA seed sequence and assesses for complementarity that suggests miRNA binding at the gene target. Gene ontology (sources; Reactome, Wikipathways, BioCarta) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were carried out in cluster profiler on the predicted target genes.^[86]35 The GO terms and KEGG pathways with adj. P < .05 were considered to be significantly enriched. To further assess miRNA-target regulatory network interactions that are not evident through pairwise analyses, we used the web tool MIENTURNET, which utilizes network-based visualization and analysis of miRNA-target interactions.^[87]36 The list of all DE myonuclear miRNAs was provided as input, which produced a network of miRNA-target interactions along with topological properties of all nodes. High-Resolution Respirometry For mitochondrial respiration experiments, 6–8-mo-old female C57BL/6 J mice were randomized to synergist ablation (n = 7) or sham surgery (n = 7). After 72 h, mice were euthanized via CO[2] inhalation followed by cervical dislocation. Plantaris muscles were removed, cut lengthwise into small samples, and immediately placed in ice-cold biopsy preservation solution (BIOPS) relaxing and preservation solution (10 m m Ca-EGTA buffer, 0.1 µm free calcium, 20 m m imidazole, 20 m m taurine, 50 m m K-MES, 0.5 m m DTT, 6.56 m m MgCl2, 5.77 m m adenosine triphosphate (ATP), 15 m m phosphocreatine, pH 7.1). Fiber bundles were mechanically separated using forceps in a petri dish on ice, followed by transfer intro BIOPS containing saponin at a final concentration of 50 μg/mL. Samples were agitated on ice for 30 min, then transferred to MiR05 mitochondrial respiration medium (Oroboros Instruments, Innsbruck, Austria) for a 10 min wash. High-resolution oxygen consumption measurements were conducted at 37°C in MiR05 supplemented with 20 m m creatine monohydrate using an Oroboros O2k Oxygraph (Oroboros Instruments). Fiber bundles of ∼2 mg wet mass were added to O2k chambers, and a substrate inhibitor titration (SUIT) protocol was performed as follows: 10 m m pyruvate, 2 m m malate, and 10 m m glutamate (complex I leak); 4 m m ADP (to initiate OXPHOS); 10 µm cytochrome c (to test integrity of the mitochondrial membrane); 10 m m succinate (to initiate electron flow through complex II); 250 n m titrations (up to 1.5 μm) of carbonyl cyanide p-trifluoromethoxyphenylhydrazone (FCCP) (maximum respiration/ET capacity); 10 µm rotenone (to inhibit complex I); and 5 µm antimycin A (to inhibit complex III for background correction). All samples had a cytochrome c response < 15%. Oxygen flux (JO2), expressed as pmol/s was normalized to sample wet weight, and differences between sham and MOV samples were assessed using an unpaired t-test, with significance set at P < .05. Results MOV Robustly Alters the Contribution of the Myonuclear DNA Methylome to Activate and Repress Gene Expression Mechanical overload induced hypermethylation of at least 1 CpG in 3988 promoters, 5092 exons, 7343 introns, and 3242 intergenic CpG sites and hypomethylation of 17 812 promoters, 13 722 exons, 12 508 introns, and 16 895 intergenic CpG sites across the genome ([88]Figure 1A and [89]B). Binding and expression target analysis revealed 116 upregulated and 203 downregulated gene targets associated with significant differentially methylated sites ([90]Figure 1C). Gene set enrichment of upregulated and downregulated targets from BETA indicated biases toward upregulation of growth-associated gene ontologies and downregulation of mitochondrial and metabolic gene ontologies ([91]Figure 1D). Specifically, downregulated targets included genes encoding mitochondrial respiratory chain complex subunits (ie, NADH: ubiquinone oxidoreductase subunits A9, A10, B9, Ndufa9, Ndufa10, Ndufb9 and cytochrome c oxidase subunits 4I1 and 5B, Cox4i1, Cox5b). Upregulated targets included growth factor genes, including nerve growth factor (Ngf), heparin-binding EGF-like growth factor (Hbegf), and fibroblast growth factor receptor 4 (Fgfr4). Examples of RNA-seq and RRBS integration are provided for a metabolism-related gene (6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3, Pfkfb3, [92]Figure S1A) and a cytoskeleton-related gene (obscurin like cytoskeletal adaptor 1, Obsl1, [93]Figure S1B). Both RNA-seq read coverage data (top) and RRBS % methylation per CpG site (bottom) are shown. Figure 1. [94]Figure 1. [95]Open in a new tab Myonuclear DNA methylation directly downregulates expression of genes encoding mitochondrial respiratory chain proteins. Integrative analysis combining reduced representation bisulfite sequencing (RRBS) and RNA sequencing (RNA-seq) data provides insight into the genes directly regulated by DNA methylation following 72 h of mechanical overload (MOV)-induced muscle hypertrophy. (A) Overview of myonuclear DNA methylome showing the number of significantly differentially methylated CpG sites falling within the promoter (2000 bp upstream of transcription start site), exonic, and intronic regions as well as CpG islands; MOV (n = 3) versus sham (n = 2). (B) Schematic showing the limitations of methylation analysis where significantly differentially methylated CpG sites are assigned to a specific gene because of promoter region. CpG sites may fall outside of the promoters and still participate in the regulation of genes in its vicinity. (C) BETA integration analysis of myonuclear RRBS and RNA-seq data comparing up- and downregulated genes relative to background (significance indicated by P-values in parentheses). (D) Chord plot showing up- and downregulated genes identified by BETA integration corresponding to significantly enriched gene ontologies; legend represent log[2] fold difference shown in red/blue heatmap for RNA-seq data [MOV (n = 4) versus sham (n = 3)]. Small RNA-seq Reveals a Modest Increase in the Proportion of Myonuclear miRNAs After 72 h of MOV We investigated the relative biotype distribution of small RNAs in myonuclei using small RNA-seq ([96]Figure 2A and [97]B). Although miRNA abundance dominated proportions of small RNAs present in myonuclei, small RNA biotypes other than miRNAs were also detected. The fraction of transcriptome-aligned reads mapping to different biotypes of RNAs in sham and MOV myonuclei varied significantly with respect to miRNAs, protein coding RNAs (mRNAs), and circular RNAs ([98]Figure 2C). MOV myonuclei had a higher proportion of miRNAs compared to sham (44.74% versus 39.69%, P = .018), while protein-coding RNA abundance was lower in MOV compared to sham myonuclei (11.17% versus 24.38%, P = .016). Detection of mRNA species is consistent with other studies showing uniquely mapped protein-coding RNAs by small RNA-seq.^[99]37,[100]38 Figure 2. [101]Figure 2. [102]Open in a new tab Differentially expressed (DE) miRNAs and prediction of biological processes targeted by DE miRNAs. The fraction of transcriptome-aligned reads mapping to different biotypes of RNAs in (A) sham (n = 3) and (B) MOV (n = 4) myonuclei. Small RNA biotypes other than miRNAs are detected. (C) Differences in relative biotype distribution between sham and overload tested using an independent t-test. snoRNA: small nucleolar RNA, tRNA: transfer RNA, piRNA: Piwi-interacting RNA, snRNA: small nuclear RNA, and lincRNA: long intergenic noncoding RNAs. *P < .05. (D) Volcano plot indicating mean expression level for each miRNA, red dots represent significant DE miRNAs (adj. P < .05). (E) Network-based pathway enrichment analysis of genes identified in silico as myonuclear-miRNA targets of myonuclear DE genes. Arrows represent relationships and color represents the adj. P-value. DE Myonuclear miRNAs After 72 h of MOV Target Genes Associated With Mitochondrial Function and Metabolism A total of 48 miRNAs were differentially expressed (21 upregulated, 27 downregulated) in myonuclei from MOV muscle compared to sham ([103]Figure 2D). Of the top 10 miRNAs by expression level in myonuclei, 6 were differentially expressed with MOV relative to sham ([104]Table S1). Targets of DE miRNAs identified using miRTarBase that corresponded with DE myonuclear mRNAs were found to be significantly enriched for several metabolism-related pathways, including TCA cycle and respiratory electron transport, pyruvate metabolism, and ATP synthesis ( [105]Figure S2A). Using the in-silico target prediction tools, miRanda and RNAhybrid, a total of 13,651 genes were identified as potential targets of all the DE myonuclear miRNAs. Given that a single miRNA is predicted to have many target genes, there was significant overlap between DE miRNA target genes. Of this list, 1618 genes were identified as DE following 72 h of MOV ([106]Table S2). Gene ontology analysis of these genes indicated enrichment of pathways for generation of precursor metabolites and energy positive regulation of the cell, adhesion, and energy derivation by oxidation of organic compounds ([107]Table S3). Further analysis of these pathways showed abundant activity in the central metabolic pathways, including cellular respiration and mitochondrial ATP synthesis coupled with electron transport ([108]Figure 2E). Network-Based Analysis Reveals miRNA-Based Regulation of Metabolic and Ubiquitin Ligase Genes The MIENTURNET web tool^[109]36 was used to perform regulatory network analysis of DE myonuclear miRNAs. This tool identified 73 genes that were predicted to be significant miRNA network-based targets ([110]Table S4). Of these genes, 12 were DE in our myonuclear mRNA-seq, including peroxisome proliferator-activated receptor alpha (Ppara) and solute carrier family 16 member 13 (Slc16a1), which are both involved in lipid metabolism^[111]39,[112]40 ([113]Figure 3A). These genes, which were significantly downregulated in the myonuclear mRNA-seq analysis, were identified as predicted targets of several DE myonuclear miRNAs. Interestingly, the miRNAs targeting these downregulated genes were also significantly lower in MOV compared to sham ([114]Figure S2B and C), which may represent regulatory functions of nuclear miRNAs beyond post-transcriptional gene silencing. Other significant miRNA network-based DE targets in our myonuclear mRNA-seq included ubiquitin protein ligase E3 component n-recognin 3 (Ubr3) and WW domain containing E3 ubiquitin protein ligase 1 (Wwp1), both of which encode E3 ubiquitin-ligase molecules ([115]Figure 3A).^[116]41,[117]42 Both genes were significantly downregulated in MOV compared to sham in the myonuclear mRNA-seq analysis. The miRNAs targeting these genes were significantly upregulated in myonuclei in MOV relative to sham ([118]Figure S2D). Figure 3. [119]Figure 3. [120]Open in a new tab Downregulation of oxidative metabolism genes following mechanical overload accompanied by reduced mitochondrial respiration of permeabilized plantaris muscle fiber bundles. (A) Interaction graph of selected miRNAs (blue circles) and their target genes (yellow circles) built in MIENTURNET. (B) Overview of substrate inhibitor titration (SUIT) protocol showing mean ± SD of oxygen flux (JO2) for each measurement state in sham (n = 7) and MOV (n = 7) samples [complex I leak; complex I oxidative phosphorylation (OXPHOS); complex II OXPHOS; combined complex I + II OXPHOS; ET Capacity: electron transfer capacity (maximum, noncoupled state)]. Data in (B) shown in (C-G) with individual data points and statistical analysis (unpaired t-test) between groups. Ns: not significant, *P < .05, **P < .01. Reduced Mitochondrial Respiration Following 72 h of MOV To assess whether the signature of downregulated myonuclear metabolic pathways following acute MOV coincides with an effect on mitochondrial function, we assessed mitochondrial respiration in permeabilized plantaris muscle fibers following 72 h of sham or MOV ([121]Figure 3B). Oxygen flux (JO2), a direct measure of oxygen consumption rate, was significantly lower in MOV compared to sham during several measurement states (complex I leak: P = .011; combined complex I + II-linked OXPHOS: P = .009; and ET capacity: P = .027) ([122]Figure 3C–E). Although OXPHOS was lower in MOV compared to sham during complex I and complex II, differences were not statistically significant (P = .060 and .108, respectively) ([123]Figure 3F and G). Discussion Skeletal muscle myofibers form a multinucleated syncytium with the ability to quickly enhance cellular size and DNA content in response to an acute MOV stimulus.^[124]43,[125]44 Resident myonuclei dramatically upregulate transcriptional output during early stages of MOV-induced hypertrophy (72 h).^[126]45 We previously interrogated myonuclear DNA methylation via RRBS and found evidence of growth-related regulation by myonuclei following 72 h of MOV. This signature was characterized by hypomethylation of mTOR signaling genes as well as Myc, a master regulator of ribosome biogenesis.^[127]19 Furthermore, we performed myonuclear RNA-seq, which identified enrichment of cell cycle regulators including Myc, and elevations in extracellular matrix (ECM) and RHO-GTPase genes in response to 72 h of MOV.^[128]20 In this work, we significantly expand our prior work to reveal the coordinated methylation and transcriptional regulation of metabolic genes with MOV that corresponds to lower mitochondrial respiration, and is complemented by the myonuclear miRNA landscape. Although there is not always a direct relationship between mRNA levels, protein levels, and ultimately phenotype and function in skeletal muscle,^[129]46 we show here that early molecular responses during rapid hypertrophy correspond well with the measurable cellular outcome of mitochondrial respriation. In general, significantly hypermethylated CpG sites that fall within the promoter region tend to repress a target gene and vice versa for the hypomethylated promoter CpG sites. Oftentimes, CpG sites fall into regions of the genome that may be the promoter of one gene but in close proximity of another gene upstream. In such cases, the role of the CpG sites in the promoter of one gene may be similar or different compared to its role downstream of another gene. Also, the overall number of differentially methylated sites within a region may impact the magnitude of gene expression regulation as well, even if that function is not yet mechanistically defined for a given gene. As such, using an integrated analysis by combining the methylomics data with the transcriptomics data facilitates a more comprehensive understanding of epigenomic regulation. One approach that allows for mismatched sample sizes across analyses is BETA, which was first developed for integrating ChIP-seq and RNA-seq datasets. In the current investigation, after integration of the RRBS data with the myonuclear RNA-seq, our data suggest that epigenetic modifications in myonuclear DNA may drive the observed reductions in oxidative metabolism gene expression. Integrative analysis specifically points to downregulated genes induced by DNA hypermethylation, and these genes are enriched for pathways representing mitochondrial respiratory chain complexes I, III, and IV and mitochondrial organization. Furthermore, upregulated GOs included growth factor activity and ECM organization, suggesting that growth-associated alterations in myonuclear gene expression may be linked to changes in DNA methylation as well. Together, upregulated and downregulated targets from BETA analysis reveal that epigenetic modifications to the DNA methylome contribute to modulation of the myonuclear transcriptome that aligns with altered muscle metabolism. A novel finding of the present work is the identification of mature miRNAs within myonuclei and our observation that the relative distribution of myonuclear miRNAs is increased after 72 h of MOV. Notably, miRNAs are transcribed into primary transcripts (pri-miRNAs) and processed into precursor miRNAs (pre-miRNAs) in the nucleus by Drosha, followed by export to the cytoplasm, where they are processed into mature miRNAs by Dicer.^[130]47 Traditionally, microRNA function has been thought to be restricted to the cytoplasm; however, deep sequencing of nuclear and cytoplasmic pools independently showed that a majority of miRNAs were imported into the nucleus,^[131]48 and nuclear-cytoplasmic shuttling machinery has also been identified.^[132]7,[133]49 In the current study, we identified several DE miRNAs by small RNA-seq of isolated myonuclei in MOV compared to sham. Target prediction using either a standard miRNA-target interaction database or in silico prediction demonstrated that mRNA targets of DE myonuclear miRNAs were enriched for oxidative metabolism-related pathways. The miRNA-mRNA interaction network also included key metabolic regulatory genes as well as ubiquitin ligase E3 components. While some miRNAs that were elevated after MOV were predicted to target downregulated genes, we also observed examples of downregulated miRNAs predicted to target downregulated genes. This could potentially be explained by the fact that both RNA activation and transcriptional gene silencing have been demonstrated in nuclei via recruitment of chromatin modifying proteins to chromosomal DNA, altering DNA methylation.^[134]12,[135]15,[136]50 Taken together, these data suggest that in addition to classical post-transcriptional silencing by 3′UTR targeting, myonuclear miRNA-mediated adaptation to MOV may also involve epigenetic mechanisms that contribute to the transcriptional control of myonuclear genes implicated in skeletal muscle metabolism as well as protein degradation. Our data provide further support for the contribution of miRNA-mediated mechanisms to the molecular responses to exercise.^[137]46 It is important to note that while cytoplasmic and nonmyonuclear miRNAs may also play a role in the hypertrophic response to MOV,^[138]51 our current investigation focused on myonuclei, which represent only ∼50%-60% of the total nuclei in murine skeletal muscle.^[139]25 Interestingly, miRNA alterations in myonuclei may not coincide with changes in their miRNA distribution. For example, we found that levels of miR-1, which decrease with MOV in whole muscle,^[140]52,[141]53 were not significantly different in myonuclei between MOV and sham. This suggests that while miR-1 transcription may not change in myonuclei in response to MOV, the decrease in miR-1 in whole muscle may reflect export via exosomes, as previously reported.^[142]52 We found that downregulated oxidative metabolism signatures after 72 h of MOV were accompanied by functional reductions in mitochondrial respiration. In addition to reduced maximal respiration (ET capacity) in the noncoupled state, respiratory capacity of mitochondria was also significantly reduced in MOV muscle both with saturating levels of ADP (OXPHOS) and during respiration supported by endogenous substrates (leak). Although this observation may appear intriguing considering skeletal muscle hypertrophy and associated increases in protein synthesis are bioenergetically costly processes,^[143]54 others report similar reductions in oxidative metabolism-related gene expression^[144]55 as well as reductions in oxidative phosphorylation proteins and peroxisome proliferator acivated receptor gamme (PPARG) coactivator 1 alpha (Pgc-1α) levels in response to MOV.^[145]56 These data, along with evidence for increased glucose uptake and utilization in growing muscle cells,^[146]57,[147]58 have led to the hypothesis that hypertrophying muscle fibers reprogram their metabolism in a manner characterized by “aerobic glycolysis,” commonly known as the “Warburg effect.”^[148]22 Increased glucose utilization may allow for increased nucleotide and rRNA synthesis, which may facilitate synergist ablation-induced muscle growth by increasing translational capacity.^[149]59 Further support for this metabolic reprogramming comes from reports indicating that glucose uptake, lactate secretion, and the pentose phosphate pathway increased in mice after MOV induced by synergist ablation.[150]^60–62 In C2C12 myotubes, glucose incorporates into muscle protein and RNA, which is enhanced during hypertrophy and controlled by phosphorylation of AKT and 3-phosphoglycerate dehydrogenase (PHGDH),^[151]63 a serine-synthesizing enzyme. In our data, Phgdh mRNA is upregulated in myonuclei during MOV (adj. P = .02) and has previously been linked to rapid developmental muscle hypertrophy in vivo,^[152]64 suggesting a high rate of glucose usage and biosynthetic metabolism during MOV. Thus, reduced mitochondrial respiration following acute MOV may be a consequence of inhibited pyruvate entry into mitochondria rather than dysfunctional mitochondria.^[153]65 Interestingly, 72 h of MOV reduces myonuclear Myh2 mRNA content (oxidative fast-twitch myosin 2A), whereas prolonged MOV (8 wk) ultimately induces an aerobic fiber-type shift toward myosin 2A fibers alongside 2A hypertrophy.^[154]66 Despite these observations, the mechanisms regulating energy metabolism reprogramming in rapidly hypertrophying skeletal myofibers remained unclear. The present work furthers our understanding of metabolic gene regulation in response to MOV. Our study had limitations that are worth addressing. First, our sample size for the integrative analysis using BETA was limited and exclusive to females, which likely reduced the power to identify significant genes and subsequently the significant pathways associated with MOV. Still, analytical integration of multi-omics can enhance statistical power and synergistic insights compared to single-omics approaches.[155]^67–70 Another limitation of the study is the supraphysiological nature of the synergist ablation model, which allows for very rapid hypertrophy, limiting the generalizability of the current findings to humans where muscle hypertrophy is typically slower. Worth noting, however, is that functional hypertrophy of type 2 muscle fibers in humans can be as high as 8%-14% per week during short-term tapering following heavy training.[156]^71–74 Synergist ablation induces inflammation associated with the surgical procedure.^[157]75,[158]76 Since inflammatory stress and signaling can result in reduced mitochondrial respiration,^[159]77 it is possible that this contributes to the synergist ablation-induced effects on mitochondrial metabolism. On balance, the epigenetic and transcriptional data presented here suggest that the reduction in mitochondrial metabolism is a “programmed” response, and not necessarily due to acute mitochondrial damage by inflammation. By integration of DNA methylation and myonuclear RNA-seq data, we demonstrate that epigenetic modifications in myonuclear DNA may play a role in the downregulation of nuclear oxidative metabolism gene expression with MOV. Furthermore, small RNA-seq profiling revealed a hitherto undescribed involvement of myonuclear miRNAs in the regulation of myonuclear oxidative metabolism-related gene expression. Our data provide new avenues for exploration of the molecular mechanisms driving epigenetic modifications in myonuclei. Supplementary Material zqad062_Supplemental_Files [160]zqad062_supplemental_files.zip^ (1.5MB, zip) Acknowledgments