Anthropometric and physiological data

The data summarized in Table 2 confirmed that the exercise was moderate and athletes reacted as normal healthy subjects [17].

Flow cytometry analysis

miRNA and mRNA differential expression profiles

The miRNA and mRNA expression profiles in the whole blood for each time point and each athlete were determined using microarray analysis. PAXgene RNA tubes enable the isolation of intracellular RNA of circulating leukocytes, including B cells, T cells, neutrophils, monocytes, and other less abundant cell types. Furthermore, a large proportion of reticulocyte-derived globin mRNA is prepared from PAXgene blood RNA tubes, as it has been demonstrated previously [19]. Affymetrix GeneChip Human Gene 1.0 ST arrays contain both miRNA and mRNA probes and are capable of measuring about 20,000 mRNAs and 200 miRNAs.

The data are presented in Additional file 1. 298 mRNAs and 5 miRNAs were changed significantly (at least 40% change with the adjusted P-value threshold of 0.05), including hsa-miR-21-5p, hsa-miR-24-2-5p, hsa-miR-27a-5p, hsa-miR-181a-5p and hsa-miR-181b-5p. Remarkably, hsa-miR-24-2-5p is clustered with hsa-miR-27a-5p, and hsa-miR-181a-5p is clustered with hsa-miR-181b-5p. Consistently, the clustered miRNAs exhibited similar expression profiles over time (see Additional file 1 and Figures 1, 2, 3 and 4).

Pathway analysis of differentially expressed mRNAs

Pathway analysis of mRNA targets of differentially expressed miRNAs

It has been demonstrated that the same mRNAs can be targeted by more than one miRNA which provides for more efficient and specific regulation [20, 21]. We found 49 mRNAs which are known to be validated targets for 2 or even 3 differentially expressed miRNAs. They have higher potential to be involved in exercise-induced regulation. Table 6 shows the pathways enriched with some of these mRNAs. Notably, these exercise-relevant pathways (cell death, stress response, proliferation) comprise significant number of intersecting mRNAs. Figure 5 presents the identified regulatory miRNA-mRNA network for all 3 revealed pathways. Interestingly, some of these genes are known to interact with each other. Namely transcription factor MYC was reported to be functionally associated with RNA helicase DDX3X [22], apoptosis regulator BCL2 [23] and tumor suppressor BRCA1 [24]. BRCA1 in turn interacts itself with BCL2 [25] and transcription factor E2F1 [26]. The presented data support the regulatory role of identified miRNAs in response to exercise.

Dynamically regulated miRNA-mRNA networks

We monitored the transcriptome expression level before and following exercise and this allowed us to reveal dynamically regulated miRNA-mRNA networks. We used a two-step approach to identify the mRNA targets for the differentially expressed miRNAs. First, we looked for anti-correlating groups of miRNAs and mRNAs expression of which over time tended to have opposite profiles. Target mRNA degradation is one of the mechanisms of miRNA action when their perfect complementarity occurs [1, 27]. Thus, the second step of our strategy was either analysis of published data and selection of experimentally validated target mRNAs for a given miRNA or theoretical prediction of miRNA targets based on their complementarity, using one of the most popular web resources TargetScan [28]. The expression profiles and biological function of selected miRNAs and mRNAs were analysed in more detail. Based on this analysis, our final identified miRNA-mRNA pairs have a high probability of being involved in the regulation of exercise-related physiological processes.

hsa-miR-21-5p

MiRNA hsa-miR-21-5p demonstrated different expression profiles over time (adjusted P-value 0.0039) however remarkably anti-correlating with experimentally validated target mRNAs TGFBR3 (adjusted P-value 1.079E-06), PDGFD (adjusted P-value 1.28E-06) and PPM1L (adjusted P-value 6.78E-07) (Figure 1). RHOBTB3 mRNA (adjusted P-value 0.002) predicted by TargetScan to be a potential target for hsa-miR-21-2-5p behaved similarly. We observed the up-regulation of hsa-miR-21-5p one hour after exercise. The differences in kinetics can be potentially explained by individualities of each athlete.

The expression level of hsa-miR-21-5p itself is known to be stress-responsive and play an important role in heart failure [29] and renal ischemia reperfusion injury [30]. Notably the up-regulation of circulating hsa-miR-21-5p was recently reported to occur in plasma upon exercise [6]. The overall action of hsa-miR-21-5p has been described by several authors to be strongly anti-inflammatory [6, 31]. Note, that 60 min into relaxation, there was an up-regulation of hsa-miR-21-5p in all 8 subjects. This may reflect the self-protective anti-inflammatory reaction to exercise.

We identified four pairing targets for this miRNA namely TGFBR3, PDGFD, PPM1L, and RHOBTB3. TGFBR3 is a transforming growth factor (TGF)-beta type III receptor, mRNA of which is known to be up-regulated in the peripheral blood leukocytes in allograft rejection-prone recipients after intestinal transplantation thus mediating innate and adaptive inflammatory functions of leukocytes [32]. PDGFD encodes for platelet derived growth factor D [33], a member of the platelet-derived growth factor family which can regulate many cellular processes, including cell proliferation, apoptosis, transformation, migration, invasion, angiogenesis and metastasis [34]. PPM1L encodes a protein phosphatase gene [35] responsible for the regulation of stress-activated protein kinase signaling cascade and apoptosis [36]. Finally, the fourth identified target for hsa-miR-21-5p is RHOBTB3 mRNA [33] which encodes for the Rho GTPase regulating the membrane traffic of proteins [37]. Interestingly, this mRNA proved to be a blood biomarker of psychosis and shows a decreased expression level in high hallucinations states [38].

hsa-miR-24-2-5p

MiRNA hsa-miR-24-2-5p (adjusted P-value 0.00017) was up-regulated immediately after exercise, then tended to decrease during the recovery period except athlete D (Figure 2). In our study for the first time we report the reaction to exercise of this miRNA which is known to have protective effects on myocytes against myocardial ischaemia/reperfusion-induced apoptosis [39].

MYC mRNA (adjusted P-value 0.00013) which is known to be the target for this miRNA from the literature, and KCNJ2 mRNA (adjusted P-value 0.00023) predicted by TargetScan to be a potential target essentially followed anti-profile of hsa-miR-24-2-5p. The protein encoded by the MYC gene is a multifunctional, nuclear phosphoprotein that plays a role in cell cycle progression, apoptosis and cellular transformation. It functions as a transcription factor that regulates transcription of specific target genes. Interestingly, hsa-miR-24-2-5p is known to be up-regulated during hematopoietic cell terminal differentiation suppressing MYC expression [40]. Thus the increase of this miRNA we observed (Figure 2) may reflect the regulation of hematopoiesis upon exercise. Remarkably, MYC mRNA is among those which are regulated by 2 differentially expressed miRNAs, namely hsa-miR-21-5p and hsa-miR-24-2-5p (Figure 5), however its expression profile is anti-correlated with only the profile of hsa-miR-24-2-5p (Figures 1 and 2).

KCNJ2 protein is an integral membrane protein and inward-rectifier type potassium channel participating in establishing the action potential waveform and excitability of neuronal and muscle tissues [41]. This mRNA expressed in peripheral blood lymphocytes is a biomarker for Parkinson’s disease [42].

hsa-miR-27a-5p

MiRNA hsa-miR-27a-5p is clustered with hsa-miR-24-2-5p and behaved similarly to it increasing after exercise and decreasing during the recovery period except athlete D (Figure 3) with the adjusted P-value 0.00012. This miRNA was reported to promote myoblast proliferation by reducing the expression of myostatin [43].

The only mRNA target identified is predicted by TargetScan ST3GAL6. The encoded protein belongs to the sialyltransferase family and is responsible for the synthesis of selectin ligands [44].

hsa-miR-181a-5p

hsa-miR-181a-5p tended to increase after exercise and then to down-regulate during the first as well as the second period of the relaxation time (Figure 4) with an adjusted P-value of 5.83E-05. The observed differential expression of hsa-miR-181a-5p in our athletes is consistent with previously published results [7, 8]. This miRNA is characterized as a regulator of hematopoietic lineage differentiation [45] and a modulator of T cell sensitivity and selection [46]. Radom-Aizik showed up-regulation of this miRNA after 30 min interval exercise. They related it to increased T cell responsiveness and reduced susceptibility to infection due to physical activity. In our study 7 subjects showed up-regulation immediately after exercise (Figure 4).

The mRNAs ROPN1L (adjusted P-value 0.00024) and SLC37A3 (adjusted P-value 0.0019) were previously validated to be targets for hsa-miR-181a-5p and demonstrated pronounced anti-correlation with the miRNA expression profile. The ROPN1L gene encodes a member of the ropporin family. The encoded protein is involved in the targeting towards specific physiological substrates of Protein Kinase A, regulating glycogen, sugar, and lipid metabolism [47]. The SLC37A3 protein belongs to transmembrane sugar transporters and is responsible for sugar metabolism [48].

Ethical approval and study participants

Eight national level ski athletes took part in this study. None of them suffered from acute or chronic diseases or reported intake of medication. Participants were informed about the nature, purpose, and potential risks of the experiments and signed an informed consent statement approved by the ethics committee of Scientific Research Center Bioclinicum (Moscow, Russia).

Anthropometric measurements

Height, weight, medical historical data and resting vital signs were recorded at the time of enrolment.

Exercise test protocol

In order to determine the VO_2max values, each subject performed a treadmill test with an incremental step protocol until exhaustion as described previously [47]. VO_2max was calculated as described [49]. The anaerobic threshold was calculated using the standard V-slope method [50].

Two weeks later, athletes participated in the main exercise, consisting of running at 80% VO₂ peak for 30 min on a treadmill. The exercise was performed during the morning hours (between 8 and 11 a.m.) keeping the exact test time for each participant constant.

Lactate concentration

Lactate concentration in capillary blood was measured electrochemically using the automated analyzer Biosen C_Line (EFK Diagnostic, Germany).

Analysis of aminoacids

Sera derived from the venous blood samples were analyzed for aminoacids in Genome Analysis Centre (Helmholtz Zentrum, Munich, Germany) using the BIOCRATES AbsoluteIDQ p150 Kit (Biocrates Life Sciences AG, Innsbruck, Austria) in a combined FIA-MS/MS and LC-MS/MS assay as recommended by the manufacturer and described previously [51].

Blood sampling

Venous blood was collected at four time points during the exercise. A 20-gauge intravenous catheter was placed antiseptically into a dorsal hand vein or a vein in the distal forearm as dictated by favourable anatomy using the Seldinger technique and then secured with tape. During ME 2.5 ml of blood was collected in PAXgene blood RNA tube, 7 ml in a Serum separator tube (BD, USA) and 4.5 ml in a tube containing buffered tri-sodium citrate (BD, USA) for flow cytometry analysis at baseline (prior to exercise testing) and immediately post-exercise (within 1 min of completion of exercise testing). After 30 min of rest and after 60 min of rest following exercise testing 2.5 ml of blood was collected in PAXgene blood RNA tubes.

RNA extraction

According to the Affymetrix Manual P/N 701880 Rev. 4 total RNA was extracted using the PAXgene Blood RNA kit as recommended by the manufacturer. RNA concentrations were determined by the Nanodrop photometer (NanoDrop, USA). RNA quality was checked using the Agilent Bioanalyser 2100 System (Agilent Technologies, USA). For all samples RNA integrity number (RIN) was greater than 7.

Flow cytometry analysis

Flow cytometry analysis of blood samples was performed using fluorescently labeled antibodies against B and T cell receptors and natural killer (NK) cell markers. Cells were labeled with antibodies against CD3 (FITC), CD4 (PE), CD8 (PE), CD16 (FITC), CD56 (PE), CD19 (FITC) (Sorbent, Russia), where NK cells were distinguished from the rest of the lymphocytes via positive expression of CD56 and negative expression of CD3.

The samples were analyzed on a FACScan Calibur flow cytometer (BD Biosciences, USA) and leukocytes were gated based on forward and side scatter properties. Events in the range of 40,000–200,000 were collected depending on the occurrence of the investigated leukocyte population, and analyzed with CELLQuest Pro analysis software (BD Biosciences, USA). To ensure flow cytometric standardization, the voltage settings were updated daily using ‘Calibrate Beads’ (BD Biosciences, USA).

Microarray analysis

RNA samples were prepared according to manufacturer’s instructions (Affymetrix Manual P/N 701880 Rev. 4) as described elsewhere [52]. The samples were hybridized on GeneChip Human Gene 1.0 ST Arrays containing both miRNA and mRNA probes (Affymetrix, USA) for 16 h at 45°C. Arrays were washed to remove non-specifically bound nucleic acids and stained on Fluidics Station 450 (Affymetrix) using FS450_0007 protocol followed by scanning on a GeneChip Scanner 3000 7G system (Affymetrix). The microarray CEL files have been deposited in the GEO database (accession GSE46075).

Microarray data processing

Microarray data was processed using bioconductor [53] xps package implementation of RMA [54]. At the first step background correction was performed based on a global model for the distribution of probe intensities [54]. Then a quantile normalization algorithm [55] (so-called probe-level normalization) was applied to the preprocessed data. Finally, fitting a robust linear model using Tukey’s median polish procedure [56] was done to convert probe intensities to the expression levels of probesets.

The statistical analysis of microarray data was performed using bioconductor [53] package limma [57]. The analysis was based on a generalized linear model [57, 58] approach. In this approach one constructs a linear data model with a structure determined by the experiment layout, and then fits this model to the actual data. The linear model is defined in terms of a so-called design matrix. The number of rows in this matrix coincides with the number of experiment samples, and the number of columns coincides with the number of factors that have an essential influence on the measurable values. The value at the i-th row and j-th column of a design matrix specifies an effect of a factor j on a sample i. Each measurable value (i.e., each probeset) in this approach is analysed independently. For each probeset a vector of its expression values E is represented in the form E = Dβ + ϵ, where D is a design matrix, β is a vector of coefficients indicating values of each factor’s actual influence on the analyzed probeset, ϵ is a vector of error, and model fitting consists in the minimization of “error term” ϵ by finding optimal coefficients β. After the coefficients β are computed for each probeset, one can test various hypotheses on the structure of considered factors. For example, in order to find probesets that are affected by a factor i, one should search for the probesets with βi statistically different from zero.

The general linear model was applied in the analysis of the studied transcriptome changes as follows. The model that takes into account both experiment time points and athletes individual features was used. Thus, the total number of factors was eleven: c1, …, c8 correspond to an expression level of each athlete in a normal state (for all samples the effect of these factors is set to 1 if a sample and a factor correspond to the same athlete, and set to 0 otherwise), and d1, d2, d3 correspond to the changes induces by the exercises, exercises and 30-minutes relaxation, and exercises and 60-minutes relaxation respectively (for all samples the effect of these factors is set to 1 if a sample and a factor correspond to the same experiment time point, and set to 0 otherwise). The total number of samples was 32: 4 samples for each athlete.

For each pair of experiment time points the detection of probesets with reliable difference between time points was performed. The probeset was considered to have a reliable difference between time points k, m if an adjusted p-value of an equality βk = βm (where βl was an actual influence of exercises at experiment time point for l = 1,2,3, and β0 = 0) was less than 0.05, and a log-fold change was greater than 0.484 (this threshold corresponds to an intensity change by more than 40%). The Benjamini-Hochberg [59] algorithm was used for multiple testing adjustment. The minimum adjusted p-value for all pairs of time points is indicated for each differentially expressed probeset in the Table of differentially expressed mRNAs and miRNAs (Additional file 1) as a statistical significance value.

Pathway analysis

Bioinformatic analysis was performed using DAVID online tool (http://david.abcc.ncifcrf.gov) as described elsewhere [60]. So all the analyzed genes were classified into several functional groups, and the groups that may be potentially associated with physiological stress were considered and listed in the tabs of excel document. P-values on the tabs are modified Fisher Exact P-Values. When members of two independent groups can fall into one of two mutually exclusive categories, Fisher Exact test is used to determine whether the proportions of those falling into each category differs by group. In DAVID annotation system, Fisher Exact test is adopted to measure the gene-enrichment in annotation terms.