Skeletal Muscle MicroRNA Patterns in Response to a Single Bout of Exercise in Females: Biomarkers for Subsequent Training Adaptation?

microRNAs (miRs) have been proposed as a promising new class of biomarkers in the context of training adaptation. Using microarray analysis, we studied skeletal muscle miR patterns in sedentary young healthy females (n = 6) before and after a single submaximal bout of endurance exercise (‘reference training’). Subsequently, participants were subjected to a structured training program, consisting of six weeks of moderate-intensity continuous endurance training (MICT) and six weeks of high-intensity interval training (HIIT) in randomized order. In vastus lateralis muscle, we found significant downregulation of myomiRs, specifically miR-1, 133a-3p, and -5p, -133b, and -499a-5p. Similarly, exercise-associated miRs-23a-3p, -378a-5p, -128-3p, -21-5p, -107, -27a-3p, -126-3p, and -152-3p were significantly downregulated, whereas miR-23a-5p was upregulated. Furthermore, in an untargeted approach for differential expression in response to acute exercise, we identified n = 35 miRs that were downregulated and n = 20 miRs that were upregulated by factor 4.5 or more. Remarkably, KEGG pathway analysis indicated central involvement of this set of miRs in fatty acid metabolism. To reproduce these data in a larger cohort of all-female subjects (n = 29), qPCR analysis was carried out on n = 15 miRs selected from the microarray, which confirmed their differential expression. Furthermore, the acute response, i.e., the difference between miR concentrations before and after the reference training, was correlated with changes in maximum oxygen uptake (V̇O2max) in response to the training program. Here, we found that miRs-199a-3p and -19b-3p might be suitable acute-response candidates that correlate with individual degrees of training adaptation in females.

Due to their widespread functions in the regulation of skeletal muscle plasticity, a central role in exercise adaptation has been attributed to myomiRs. Furthermore, various other miRs are known to regulate central aspects of muscle biology, such as cell proliferation and differentiation, protein synthesis and metabolism, which have been implicated in exercise adaptation.
However, despite the fact that several studies on circulating miRs and exercise adaptation in different settings have been published [25][26][27][28][29][30], little is known on skeletal muscle patterns in response to resistance and, particularly, endurance exercise. In addition, some of the data were only obtained in animal models; thus, their transferability to humans is questionable. Moreover, exercise protocols are highly heterogeneous, ranging from single bouts to training of several months' duration. Finally, subjects' characteristics are diverse, and, as in most areas of exercise biology, females are underrepresented in the relevant studies [31][32][33][34][35].
Here, we analyzed skeletal muscle acute responses to a single bout of non-exhaustive cycling exercise in a group of young, healthy, sedentary females. Our aims were (I) to characterize skeletal muscle miR patterns in response to a single bout of exercise, using both a targeted and an untargeted approach, (II) to analyze the potential of miR microarray analysis of a small sample as an untargeted prescreening tool to identify candidate miRs for further testing in larger, 'targeted' qPCR studies, and (III) to test whether 'acute-response' miR patterns might be promising biomarkers in the context of developing individualized exercise regimens.

Materials and Methods
Subjects and study protocol. The analysis described in this paper was part of the 'iReAct' (individual response to physical activity) study that has previously been described [62]. Briefly, n = 42 (F: n = 30, M: n = 12) young (age 20-40 years), healthy, sedentary subjects underwent a structured training program on a bicycle ergometer. Due to the low number of males, only the group of female subjects was included in the current study, as already described [63]. Subjects' baselineVO 2 max was 30.0 +/− 3.2 mL/kg·min; min: 24.2 mL/kg·min; max: 36.9 mL/kg·min. After enrollment, they were randomized to either start with a six-week period of moderate-intensity continuous training (MICT), followed by six weeks of high-intensity interval training (HIIT), or to accomplish the training program in reverse order. In both cases, subjects exercised three times per week. At baseline (T0), after the first six weeks of training (T1) and at the end of the program (T2), diagnostics, including spiroergometry to assess maximum oxygen uptake (VO 2 max), and a standardized one-hour reference training at submaximal intensity, calculated as the average of power outputs at lactate thresholds #1 and #2, was carried out. Samples employed for assessing skeletal muscle acute response in the context of this study were taken from this reference training (at T0). Detailed descriptions of the physiological results can be found elsewhere [64,65]. Changes in miR patterns at rest with training (delta T0-T1) and potential correlations between baseline miR patterns and training adaptation, namelyVO 2 max, have previously been described [63].
Skeletal muscle biopsies. Before (pre) and two hours after the end of the reference training at baseline (post), fine-needle skeletal muscle biopsies were taken from the m. vastus lateralis as has previously been described [63]. miR isolation. Total RNA, including small RNA species, was isolated from skeletal muscle tissue using the miRNeasy Micro Kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. miR microarray analysis. Microarray analysis was carried out on samples obtained before and after the reference training from n = 6 female subjects. This group of subjects has previously been described [63] and was chosen based on maximum similarity with regard to sex, age, BMI, and baseline fitness. Affymetrix miR Array 4.0 analysis was carried out by ATLAS Biolabs, Berlin, Germany. Data were analyzed using Transcriptome Analysis Console (TAC), Version 4.0.2.15 (Thermo Fisher Scientific, Waltham, MA, USA). Based on the data format in this software, throughout this paper, individual miRs are referred to as '-3p' or '-5p' whenever possible.
Reverse transcription and qPCR analysis. Reverse transcription was carried out using the miRCURY LNA RT Kits (Qiagen, Hilden, Germany), according to the manufacturer's instructions. For qPCR analysis, the miRCURY LNA miRNA PCR Assays (Qiagen, Hilden, Germany) in conjunction with primers specific for individual miRs, were employed. Due to technical issues, qPCR analysis was not possible with samples of one of the subjects; so, only data from n = 29 participants (instead of n = 30) could be included in the analysis.
Statistical analysis. Statistical analysis was carried out using IBM SPSS Statistics Version 27.0 (IBM, Armonk, NY, USA). Data were tested for normal distribution using the Shapiro-Wilk test. For normally distributed data, mean values were compared between groups using t-tests for paired samples. Otherwise, Wilcoxon signed-rank tests were employed. Correlation analyses were carried out by means of the Spearman method, observing Cohen's guidelines for evaluating the correlation coefficient r [66]. Based on the pilot, hypothesis-generating character of our study, Bonferroni correction of significance levels was not implemented. To test for concordance between microarray and qPCR data, we applied the 'degree of concordance' γ as defined by Goodman and Kruskal. All statistical tests were two-sided, with significance defined as p ≤ 0.05 (*; significant), 0.01 (**; very significant), or 0.001 (***; extremely significant).
KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis. KEGG pathway analysis was carried out using the TarBase v7.0 database within the online platform DIANA-miRPath v3.0 [67].
Microarray expression patterns of miRs with fold changes ≥4.5 or ≤−4.5-untargeted approach. Next, we applied an untargeted strategy, screening our miR data for species with fold changes of ≥4.5 or ≤−4.5. Using this approach, we identified n = 55 miRs (n = 20 upregulated and n = 35 downregulated). Most fold changes were statistically significant (Supplementary Table S1).
Microarray expression patterns of miRs with fold changes ≥4.5 or ≤−4.5-untargeted approach. Next, we applied an untargeted strategy, screening our miR data for species with fold changes of ≥4.5 or ≤−4.5. Using this approach, we identified n = 55 miRs (n = 20 upregulated and n = 35 downregulated). Most fold changes were statistically significant (Supplementary Table S1).
Microarray expression patterns of miRs with fold changes ≥4.5 or ≤−4.5-untargeted approach. Next, we applied an untargeted strategy, screening our miR data for species with fold changes of ≥4.5 or ≤−4.5. Using this approach, we identified n = 55 miRs (n = 20 upregulated and n = 35 downregulated). Most fold changes were statistically significant (Supplementary Table S1).
Hippo signal transduction pathway (hsa04390) 1.11 × 10 −13 100 -660-5p -1226-5p -2110 -2115-5p TGFβ signal transduction pathway (hsa04350) 3.88 × 10 −5 51 -132-3p -155-5p -199a-5p -497-5p 11 steroid biosynthesis (hsa00100) 0.00016 12 -660-5p 14 protein processing in the endoplasmic reticulum (hsa04141) 0.00059 109 -197-5p -497-5p 15 Hepatitis B (hsa05161) 0.0029 88 -155-5p -497-5p 16 chronic myeloic leukemia (hsa05220) 0.0131 53 -155-5p -497-5p qPCR analysis of miR expression patterns in a larger cohort of subjects (n = 29). To assess reproducibility of miR patterns in a larger cohort, samples of n = 29 female subjects, representing almost all the females included in the study (n = 30; one sample did not yield enough high-quality RNA) were analyzed by qPCR for the n = 15 miRs described in the preceding paragraph. As shown in Table 4 and Supplementary Figure S1, we found that when this larger group of subjects was analyzed, general trends (upregulation versus downregulation) correlated between microarray and qPCR analyses for all 15 miR species analyzed, with most differences being significant or very significant. Furthermore, we found high degrees of correlation between expression patterns of individual miR species, specifically, three groups with similar regulatory patterns (group 1: miR-1-3p, miR-133a-3p, miR-133a-5p, miR-133b, miR-499a-5p, and miR-378a-5p; group 2: miR-497-5p, miR-199a-5p, and miR-27a-3p; and group 3: miR-15a-5p, miR-18a-5p, and miR-19b-3p) were identified (Table 5). Table 3. Concordance between miR microarray and qPCR data. Summary of concordant and discordant data. Overall concordant data (mean upregulation or mean downregulation in both microarray and qPCR) are labeled in green, discordant data in red. The degree of concordance γ ((C-D)/(C + D); C: number of concordant pairs, D: number of discordant pairs) indicates concordance between the two methods for specific individuals: γ = 1: complete concordance, γ = 0: no statistical correlation, γ = −1: complete discordance).  qPCR analysis of miR expression patterns in a larger cohort of subjects (n = 29). To assess reproducibility of miR patterns in a larger cohort, samples of n = 29 female subjects, representing almost all the females included in the study (n = 30; one sample did not yield enough high-quality RNA) were analyzed by qPCR for the n = 15 miRs described in the preceding paragraph. As shown in Table 4 and Supplementary Figure 1, we found that when this larger group of subjects was analyzed, general trends (upregulation versus downregulation) correlated between microarray and qPCR analyses for all 15 miR species analyzed, with most differences being significant or very significant. Further- Figure 2. Regulation of miR-199a-3p and -5p in the skeletal muscle acute response, as assessed by microarray and qPCR analysis. Spaghetti plots show skeletal muscle expression levels of miR-199a-3p and -5p, pre and post exercise, as assessed by microarray and qPCR analysis in n = 6 subjects. Identical colors represent identical subjects. Table 4. MiR acute responses in a larger cohort. In total, 15 miRs were selected from the microarray for qPCR analysis in a larger cohort of subjects (n = 29). Fold changes and p values are displayed as indicated. Overall, n = 13 miRs were downregulated, n = 1 (miR-155-5p) was upregulated, and n = 1 (miR-132-3p) was not significantly regulated. * p ≤ 0.05, ** p ≤ 0.01, n.s.: not significant. Correlation of qPCR data with ∆VO 2 max. Finally, we correlated training-induced adaptations, namely changes inVO 2 max, with miR patterns of all subjects (n = 29), based on qPCR-derived results. Subjects either started with MICT (n = 16) or HIIT (n = 13) and then switched to HIIT (n = 11) or MICT (n = 12), respectively. Subjects'VO 2 max data have previously been published [64,65] and are summarized in Supplementary Table S2, indicating gains in relativeVO 2 max throughout the training program. As shown in Figure 3, for ∆VO 2 max during the first training period, we found positive correlations with ∆miR-199a-3p (r = 0.615; p = 0.025 *) in the HIIT group and with ∆miR-19b-3p (r = 0.518; p = 0.04 *) in the MICT group. These two correlations had also been identified in the initial screening microarray analysis. Furthermore, additional correlations were found for acute changes in miR patterns and ∆VO 2 max during the second (T1/T2) or the entire (T2) training period; however, none of these reached statistical significance (data not shown).
adaptations, namely changes in VȮ2max, with miR patterns of all subjects (n = 29), based on qPCR-derived results. Subjects either started with MICT (n = 16) or HIIT (n = 13) and then switched to HIIT (n = 11) or MICT (n = 12), respectively. Subjects' VȮ2max data have previously been published [64,65] and are summarized in Supplementary Table S2, indicating gains in relative VȮ2max throughout the training program. As shown in Figure  3, for VȮ2max during the first training period, we found positive correlations with miR-199a-3p (r = 0.615; p = 0.025 *) in the HIIT group and with miR-19b-3p (r = 0.518; p = 0.04 *) in the MICT group. These two correlations had also been identified in the initial screening microarray analysis. Furthermore, additional correlations were found for acute changes in miR patterns and VȮ2max during the second (T1/T2) or the entire (T2) training period; however, none of these reached statistical significance (data not shown).

Figure 3.
Correlations between skeletal muscle miR acute responses and training-associated changes in VȮ2max between T0 and T1. Regression analysis of acute responses of miR-199a-3p (left panel) and miR-19b-3p (right panel), as assessed by qPCR, and VȮ2max between T0 and T1. Triangles represent subjects assigned to HIIT during the first training period, and circles represent those who started with MICT. Colored symbols represent subjects included in the initial microarray screening. Spearman correlation coefficients r and p values are indicated. * p ≤ 0.05. Figure 3. Correlations between skeletal muscle miR acute responses and training-associated changes inVO 2 max between T0 and T1. Regression analysis of acute responses of miR-199a-3p (left panel) and miR-19b-3p (right panel), as assessed by qPCR, and ∆VO 2 max between T0 and T1. Triangles represent subjects assigned to HIIT during the first training period, and circles represent those who started with MICT. Colored symbols represent subjects included in the initial microarray screening. Spearman correlation coefficients r and p values are indicated. * p ≤ 0.05.

Discussion
Our data demonstrate differential expression of multiple miRs in skeletal muscle tissue in response to a single submaximal bout of exercise.
Nevertheless, it is important to be aware of the fact that miRs form a complex network, with each miR targeting a multitude of genes and with each gene product being targeted by multiple miRs [110]. Thus, in the future, functional analyses, e.g., in knockout mice, are necessary to understand the complex interplay between skeletal muscle miRs and their targets in the context of training adaptation.
We could demonstrate strong inter-individual variability with regard to skeletal muscle miR acute responses, suggesting that miRs might be suitable biomarkers for individual training adaptation. Indeed, for HIIT, we could detect correlations between miR-199a-3p acute responses and training-associated ∆VO 2 max, and similar results were obtained for MICT with regard to miR-19b-3p. Nevertheless, the facts that (I) we did not apply corrections for multiple testing, due to the pilot character of our study, and (II) we were not able to establish correlations between miR acute patterns and training adaptations in the entire cohort, suggest that these correlations might not prove true in larger and/or more heterogeneous groups. Thus, in the future, more extensive studies and, specifically, direct comparisons between male and female athletes will have to be carried out. In addition, with regard to practicability, it would be desirable to switch to circulating miRs as markers, which, as described in several studies [25,26,29], will probably be even more complex, since these represent a mixture of miRs originating not only from skeletal muscle, but from a broad variety of tissues and organs.
In conclusion, our data suggest that microarray-based prescreening of a small subset of samples, followed by qPCR validation in the entire cohort, might be a valuable instrument for establishing skeletal muscle miR patterns and their modulation by exercise. Firstly, comparison of microarray and qPCR data yielded 73.33% concordance of individual trends (upregulation versus downregulation) in the six subjects initially analyzed. Secondly, when microarray data for miRs with fold changes of 4.5 or more were validated by qPCR in the entire cohort (n = 29), we even found that 100% of the general trends (means), i.e., average up-or downregulation as detected in the microarray analysis, could be confirmed for all miRs evaluated, with the restriction that, due to a high proportion of unspecific signals, qPCR data for miR-23a-5p, -4330, -4743-5p, and -7151-3p could not be analyzed.
In addition, in the future, miR-based screening of subjects prior to starting a training regimen might be a means to predict training responses, specifically in clinical settings but also in competitive sports. In particular, such a prediction might be helpful in preidentifying subjects for which an insufficient response to 'standard' training protocols is likely. These individuals might then directly be allocated to alternative, yet-to-be-defined and 'non-standard' training regimens, from which they might potentially benefit to a higher degree.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/biom13060884/s1, Figure S1: Changes in miR concentrations in skeletal muscle tissue in response to a single bout of exercise; Table S1: miRs found to be regulated 4.5-fold or more in microarray analysis (n = 6); Table S2:VO 2 max data (means and standard deviations) of all subjects included in qPCR analysis.  Informed Consent Statement: Written informed consent was obtained from all subjects involved in the study.
Data Availability Statement: All data are available from the authors upon reasonable request.