Endurance competition affected the expression of a large number of genes, miRNAs and metabolites

¹H NMR spectra together with custom equine mRNA and miRNA microarrays were used on an experimental set of ten horses to study the effects that endurance exercise has on the blood metabolome, transcriptome and miRNome. Before delving into multi-layered data sets integration for the analysis of the regulatory network, we assessed the quality and significance of each “omic” layer individually.

A total of 54 metabolic peaks were identified, including several amino acids, energy metabolism-related metabolites, saccharides, and organic osmolytes in the plasma (Fig. 1). Among the 54 metabolic peaks identified, we observed two unassigned compounds at 3.72–3.76 and 3.85–3.88 ppm, respectively (Additional file 3).

The effect exercise had on the serum metabolome was described by multivariate statistical methods. On one hand, the principle component analysis (PCA) significantly discriminated pre and post-endurance competition samples (Monte Carlo test; p < 0.05; Additional file 4: Figure S1A). On the other hand, the orthogonal projections to latent structures (OPLS) analysis revealed that a total of 49 metabolites resonances contributed in the discrimination of the pre and post-endurance competition samples (Additional file 4: Figure S1C). Lactate, followed by methylene, N-acetyl moieties, proline, glucose and phosphocholine presented the highest contribution to class separation (Additional file 4: Figure S1C). More precisely, loadings examined on the pre-component basis as line plot emphasized that post-competition samples were characterized by higher excursion of lactate, glycerol, creatine, urea, and aromatic amino acids such as tyrosine, along with negative excursion of methylene, N-acetyl moieties, methylene esters, glucose, and phosphocholine (Additional file 4: Figure S1C).

To ensure the accuracy of the OPLS model (R ² = 0.88 and Q² = 0.92), permutations testing and cross validation were used for internal validation of the OPLS model. The explained variance, predictive capability and out-of-bag error of OPLS model remained higher (R ² = 0.57 ± 4.07; Q² = 0.94 ± 0.03 and root mean square error of prediction (RMSEP) = 0.30 ± 0.09) than those of the 1,000-permutated models (R ² = 0.53 ± 7.04; Q² = 0.16 ± 0.48 and RMSEP = 0.65 ± 0.16) to discriminate pre- and post-endurance competition samples (p < 0.001).

The relationship between metabolites and the functional map of the metabolites affected by the endurance exercise are provided in the Additional file 5: Figure S2 and Additional file 6: Figure S3, respectively. Among the most enriched pathways are those involved in glycerolipid metabolism, D-glutamine and D-glutamate metabolism, pyruvate metabolism, tricarboxylic acid (TCA) cycle and amino acids metabolism (Additional file 6: Figure S3).

Afterwards, we focused on the identification of genes and miRNAs whose expression was significantly altered at T1 (post-endurance competition) relative to T0 (baseline, pre- endurance competition), with an adjusted p-value < 0.05 after implementing Benjamini and Hochberg correction for multiple testing (Additional file 7). The application of this threshold led to the identification of 7,678 differentially expressed genes (DEGs), 3,848 of which were over-expressed and 3,830 of which were under-expressed at T1 (relative to T0). We identified 107 miRNAs differentially expressed (DEmiRNAs) when comparing pre- and post- endurance competition samples (Additional file 8). Along with orthologous human-equine DEmiRNAs, we detected 7 equine-specific DEmiRNAs.

A regulatory network of 11 metabolites, 263 genes, and 5 miRNAs has a key role in controlling the adaptation to endurance exercise

The 49 metabolite resonances that significantly discriminated between pre and post- endurance competition samples, together with the 7,678 DEGs and 107 DEmiRNAs were used to perform a multi-step approach and identify the dynamically co-regulated relationship among metabolites, genes and miRNAs that could not be detected when analyzing only single metabolites or transcripts alone (Additional file 9: Figure S4).

We performed two different complementary statistical approaches (detailed in the Methods section) to identify enriched metabolites with significantly more predictive metabolic genes among the DEGs than among other genes in the transcriptome. The enzyme-encoding genes whose protein products participate in the metabolic reaction of a given metabolite are referred to here as metabolic genes. Next, for each of the metabolic genes, we used the multiMiR package in R to build a comprehensive list of all the putative DEmiRNA regulators.

In particular, we found a total of 11 differentially expressed metabolites (Additional file 10) that were putatively regulated by a total of 263 differentially expressed genes. Interestingly, we discovered that the spatiotemporal expression of those 263 metabolic genes could be extensively regulated at the transcription and post-transcription level through a total of 79 DEmiRNAs. For instance, let-7b-5p, miR-16-5p, miR-21-5p, miR-92a-3p, and miR-192-5p regulated more than ten metabolic genes each (Fig. 2).

Ultimately, we considered a total of 11 enriched metabolites correlated with the expression of 263 metabolic genes and 5 miRNAs during exercise to create the regulatory network (Fig. 2). This regulatory network underlined several well-known biochemical interactions, i.e., the relationships between the decreased glucose concentration after endurance exercise and the metabolic genes whose proteins participate in the Forkhead box protein O (FoxO) and glucagon signaling pathways to induce gluconeogenesis, and the relationship between lactate and metabolic genes encoding proteins related to glucagon signaling and propanoate metabolism (Fig. 2). Another well-known biochemical interaction was found between succinate and the expression of multiple components of the mitochondrial oxidative phosphorylation (OXPHOS) system, including multiple subunits of ATP synthase (seven Fo subunits and the ATP5A1, ATP5B and ATP5O subunit of F1), NADH: ubiquinone oxidoreductase core subunits (ND3, ND5, ND6; complex I) and cytochrome oxidase subunits (COX4I1, COX5A, COX5B, COX6A1, COX6B1). As expected, the regulatory network also highlighted the interaction between glycerol and metabolic genes that encode enzymes participating in the pentose and glucuronate pathways, as well as acetate and many downstream genes of propanoate metabolism, which included activated protein kinase (AMPK) activators such as acetyl-CoA carboxylase (ACACB) and acyl-CoA synthetase short-chain family (ACSS2, ACSS3).

These well-known biochemical interactions were mainly regulated by miR-16-5p, miR-21-5p, miR-92a-3p, miR-192-5p and let-7b-5p within the regulatory network.

Remarkably, the “omic” integration of molecular data in our study predicted some functional links between metabolites and the transcriptome that could not be characterized by analyzing separately any of the three individual data sets. One example is the functional enrichment of glutamate through metabolic genes belonging to the glutathione metabolism, which regulates the reactive oxygen species (ROS) and nitrogen oxide species (RONS) substances (e.g. the glutathione S-transferase (GST) gene family, GST-kappa 1 (GSTK1), and two microsomal GSTs (MGST1 and MGST3). In line with the exercise-induced oxidative stress, we observed a significant interaction between tyrosine and metabolic genes that affected the ubiquinone endogenous antioxidant levels. Another example of the importance of “omic” integration was represented by the interaction between urea and genes involved in the catabolism of purine and pyrimidine nucleotides from ATP to uric acid, in which free radicals are produced and cause muscle cell damage [26]. Lastly, a significant interaction that could not be characterized by analyzing separately any of the three individual data sets was the link between lactate and the activation of the hypoxia-inducible factor A (HIF1A) transcription factor, as well as lactate genes associated with the release of lipopolysaccharides (LPS)-induced pro-inflammatory cytokines, including toll-like receptor 4 (TLR4), followed by the activation of various immune-related transcription factors (TFs) such as signal transducer and activator of transcription 3 (STAT3) and the nuclear factor kappa beta (NFkβ-RelA).

Moreover, we found that the regulatory network based on the 263 metabolic genes included TFs such as retinoblastoma-Like 2 (RBL2), followed by the cooperatively transcriptional cofactors sirtuin 1 (SIRT1), E1A binding protein P300 (EP300), CREB Binding Protein (CREBBP), and B-cell lymphoma 6 protein (BCL6; Fig. 3a). Additionally, results suggested that these TFs might (i) regulate the spatiotemporal concentration patterns of glucose and lactate (Fig. 3a); and (ii) drive or repress the expression of metabolic genes and miRNAs in a feed-forward and feedback manner [27] (Fig. 3b).

The regulatory network is driven by both transcription factors and miRNAs

As suggested above, both TFs and miRNAs might regulate endurance exercise response. In order to ensure the accuracy of TFs and miRNAs regulation within the regulatory network, we decided to use the regulatory impact factor (RIF) algorithm [28, 29]. Our aim was to identify (i) putative causal regulators, and (ii) the rewired transcriptional circuits through which the TFs and miRNAs exert their regulatory impact on the transcriptome following exercise.

Among the most enriched TFs and miRNAs, we confirmed CREBBP (RIF2 = -2.51), EP300 (RIF2 = -1.55), RBL2 (RIF1 = 1.65), SIRT1 (RIF1 = 1.63), BCL6 (RIF2 = 1.48) and let-7b-5p (RIF1 = 1.51; Fig. 2). Details on the top regulators with greatest scores are fully listed in the Additional file 11.

The relationship between blood metabolome, transcriptome and miRNome that occurs during endurance exercise could be reproduced in an independent cohort of 31 horses

As a final analysis step, multiple factor analysis (MFA) was applied to an independent cohort of 31 horses to further emphasize the reliability of cross-layer “omic” analysis to predict or explain the exercise response adaptations and importantly, to confirm the interactions within our regulatory network. This independent cohort included 13 horses sampled only at T0, and 18 other horses sampled only after completing the endurance competition (T1).

At T0 (basal time), the MFA superimposed plot of each type of data and its barycenter showed particularly high levels of variability in all “omic” layers (Fig. 4a). However, the MFA projection plot of the partial representations for each “omic” layer (metabolome, transcriptome and miRNome) onto PC1 was tightly clustered (RV-coefficient = 0.51; Fig. 4b). The two metabolites showing the highest correlation (|r²| > 0.80) on the first axis were glutamine and glucose, whereas the genes most strongly correlating with PC1 (|r²| > 0.80) were RNA polymerase 21-kDA subunit (POLR2I), torsin 3A (TOR3A), vesicle-associated membrane protein 5 (VAMP5), histone binding protein 5 (COPR5), N1-acetyltransferase family member 2 (SAT2), AT-hook containing transcription factor 1(AHCTF1) and various ribosomal proteins (Fig. 4d). The eca-miR-425, miR-106b-5p, miR-138-5p, miR-148a, and miR-20a-5p were also highly linked to the first principal component (PC) (|r²| > 0.80; Fig. 4d). Genes such as E3 ligases (UBE3C) and tripartite motif containing 65 (TRIM65) presented high correlations with PC2 (|r²| > 0.90; Fig. 4e). The detailed weights of each metabolite, metabolic genes and miRNA are provided in Additional file 12.

Although the size of the cohort used at T0 was relatively small, the MFA analysis revealed that the partial representation of the mean individuals for each “omic” layer was projected together according to their capacity to finish the endurance competition (Fig. 4a). Thus, the molecular profiles at basal time of the animals that did not finish the endurance competition were projected close to each other in the space (Fig. 4a). In light of these findings, we wanted to understand whether the molecules detected at basal time could be considered as predictive biomarkers for the elimination during the endurance competition. Although cause and effect are usually difficult to decipher, we compared basal profiles of molecules that showed high correlation values (|r²| > 0.80) on the two first dimensions of the MFA with animals that finished the endurance competition and animals that were eliminated during the endurance competition because of lameness or metabolic problems. We discovered that the POLR2I, COPR5, SAT2 and UBE3C genes were clearly up regulated (corrected p < 0.05) at basal time in animals that were eliminated during the endurance competition, whereas TOR3A and TRIM65 tended to be up regulated (corrected p < 0.10) at basal time. They also presented numerically higher glutamine and glucose concentrations at basal time (Fig. 4d).

The combined analysis of metabolome, transcriptome and miRNome on post-exercise samples then highlighted an interesting co-structure between data sets (Fig. 5). Metabolome profile accounted for the highest variance on PC1, whereas the miRNome data were projected on the positive end of the PC2 (Fig. 5b). Metabolites such as lactate and fatty acid methylene moieties presented high correlation values |r²| > 0.80 with PC1, together with let-7 family, miR-16-5p, miR-17, miR-20ab-5p, miR-21-5p,miR-26b-5p and miR-98-3p (Fig. 5d). The LPS-induced immune markers, specifically TLR4, TNF alpha induced protein 3 (TNFAIP3) and interleukin-1 receptor-associated kinase-like 2 (IRAK2), were identified among the strongest variables of the PC1 (Fig. 5d), along with the increase in expression of genes related to mitochondrial metabolism (e.g. NADH dehydrogenase (ubiquinone) 1 beta sub-complex 5; NDUFB5).

Additional file 13 shows an expanded list of the genes related to LPS-induced immune responses with its MFA impact factor at T1. Biologically interesting genes with high correlation values on the PC2 were the mitochondrial ribosomal protein L35 and L22, respectively (MRPL35, MRPL22), and mitochondrial single-stranded DNA-binding protein (SSBP1). The top molecules with greatest weights on PC1 and PC2 at T1 are listed in the Additional file 13.

Animals and samples

The data sets used for metabolome, transcriptome, and miRNome analysis have already been published by Le Moyec et al. [19] and Mach et al. [8].

Le Moyec et al. [19] evaluated the effects of long endurance exercise on the plasma metabolomic profiles, whereas Mach et al. [8] identified the whole blood miRNA-mRNA relationships specifically regulated by endurance exercise in the same horses. Altogether, we selected samples from a total of 41 pure-breed or half-breed Arabian horses (13 females and 28 geldings; mean ± SEM age: 9.7 ± 1.5) participating in a 160 km endurance competition. All horses were recruited on a volunteer basis in three different French competitions (see Additional file 1). For each of these 41 horses, the combined metabolome, transcriptome, and miRNome profiles were available. To ensure sample homogeneity, the participating horses were subjected to the same management practices throughout the endurance competition and passed the International Equestrian Federation (FEI)’s compulsory examination before the start. Animals were fed 2–3 h before the start of the endurance competition with ad libitum hay and 1 kg of concentrate pellets. During the endurance competition, all the animals underwent veterinary checks every 20 to 40 km, followed by recovery periods of 40 to 50 min (in accordance with the FEI rules of endurance riding). After each vet gate check, animals were provided with ad libitum water and hay and a small amount of concentrate pellets. The median winning speed over the entire endurance competition was 15.76 ± 1.04 km/h.

In our study, the 41 horses were divided into two non-overlapping sets: the experimental set and the validation set. The validation set came from an independent cohort of animals, to ensure that the observed effects were reproducible in a broader context. The experimental set included 10 horses sampled before the 160 km endurance competition (T0, baseline) and post- endurance competition (T1).

The validation set included 13 horses sampled solely at T0 and 18 other horses sampled solely at T1. Among the horses sampled at T0, three animals finished the endurance competition, while ten horses failed a vet gate check for lameness, metabolic disorders or tiredness (see Additional file 1). All horses sampled at T1 completed the endurance competition.

Blood samples for metabolome, transcriptome, and miRNome profiling were obtained from the jugular vein at rest (Basal, T0) and/or immediately after the end of the competition (T1). Pretreatment of the blood samples was carried out immediately after the collection because field conditions provided access to refrigeration and electrical power supply. Briefly, whole blood samples from each horse were collected in sodium fluoride and oxalate tubes for metabolome profiling in order to inhibit further glycolysis that may increase lactate levels after sampling. Whole blood drawn for plasma generation was refrigerated immediately at 4 °C to minimize the metabolic activity of cells and enzymes and kept the metabolite pattern almost stable. Clotting time at 4 °C was strictly controlled for all samples to avoid cell lyses that affect metabolome components. After clotting at 4 °C, the plasma was separated from the blood cells and subsequently transported to the lab at 4 °C and frozen at -80 °C (no more than 5 h later, in all cases). Additionally, blood samples were collected in dry tubes at the end of the endurance event for the biochemical analysis. After clotting at room temperature, the tubes were centrifuged and the harvested serum was stored at 4 °C until analysis (no more than 48 later, in all cases). PAXgene Blood RNA tubes (Qiagen, Hilden, Germany) were used for transcriptome and miRNome profiling. They were kept at room temperature for no more than 5 h in all cases, and stored at -80 °C until analysis.

¹H NMR data acquisition and statistical analysis

The ¹H NMR spectra were acquired with a Bruker® Advance II spectrometer (Bruker BioSpin, Wissembourg, France) operating at 500 MHz and using a standard water-suppressed 1D spectrum (NOESY1D) sequence for the preservation of lipid signal intensities.

The methods for sample preparation, data acquisition, data quality control, spectroscopic data-pre-processing, and data pre-processing including peak alignment, scaling and normalization are broadly explained elsewhere [19]. In contrast to Le Moyec et al. [19], where in silico metabolite identification and statistics were done by the MATLAB commercial program, we used different open-source R packages, allowing automated qualitative and quantitative metabolite characterization, calculation of spatial significance and, importantly, metabolite pathway analysis.

Peaks with just one non-zero intensity (single mass events) were removed from the normalized matrix. Additionally, a metabolite peak was considered to be detectable only if it was expressed in at least 50% of the experimental samples. The metabolite identification was then performed by using structure message of metabolites acquired from other available biochemical databases, such as human metabolome database (HMD), http://www.hmdb.ca; KEGG, http://www.genome.jp/kegg/; METLIN, http://metlin.scripps.edu/; Chemical Entities of Biological Interest (http://www.ebi.ac.uk/Databases/); and Lipidmaps (http://www.lipidmaps.org/). Metabolite assignment of each peak was considered when chemical shifts of peaks in the samples were the same as in the publicly available reference databases (with a shift tolerance level of ± 0.005 ppm), in order to counter-act the effects of measurements and pre-processing variability introduced by factors such as pH values and solvents. A manual curation for identified compounds was done by an expert in horse metabolomics [19, 20]. Afterwards, the relative abundance of each metabolite was calculated as the area of the peak [50].

Exploratory analysis of the metabolites signals were first performed by PCA, which displays the internal structure of data sets in an unbiased manner and reduces data dimensionality through linear combinations of the original variables. A PCA score plot was used to reveal the presence of outliers outside the 95% significance region of ellipse (i.e., strong outliers). Interclass PCA of plasma metabolites with pre and post- endurance competition samples as an instrumental variable was also performed, based on a Monte Carlo test with 999 replicates. Unsupervised and supervised methods were performed with a toolbox to explore NMR metabolomic data sets in the R environment [51]. Then, we performed supervised projections to latent structures-Discriminant Analysis (PLS-DA) and the OPLS analysis, which integrates an orthogonal signal correction (OSC) to identify and characterize metabolic changes induced by endurance exercise. The OSC-correction approach was conducted using DeviumWeb R package; (https://github.com/dgrapov/DeviumWeb) [52].

By extracting variation from its computed partial least squares (PLS) components that is uncorrelated (orthogonal) to the responses, OPLS produces a more interpretable regression model compared to PLS [53]. Loading plots combining the reliability and correlation from the OPLS models were used to identify differential metabolites among the pre and post- endurance competition samples. In the loading plots, signals with a positive direction corresponded to metabolites that were present at high concentrations at T1. The negative direction indicated a negative direction of metabolite values at T1. A high loading score means that the metabolite in question has an excellent ability to separate pre and post-endurance competition samples. Importantly, only metabolites with loadings plots different from zero along any OPLS principal component axis were considered to have a contribution to the model. These metabolites were considered to be biological relevant and suitable for downstream analysis.

The 10-fold within model cross validation and permutation Monte Carlo testing (n = 1,000) were applied for internal validation of the Q², R² and the RMSEP values, which represent the predictability of the model, the explained variance and the error, respectively. The firmness of the model was evaluated with R² outcomes and the precision of the prediction (Q²).

Complementary to the multivariate statistical methods to identify metabolites that were statistically different between T0 and T1 samples, analysis of variance (ANOVA) was performed to delineate whether there was a significant difference between the average values of metabolites normalized intensities between pre and post- endurance competition samples. Because of the multiple testing issues, Bonferroni corrected P-values were calculated. A significance level of corrected p < 0.05 was accepted.

This data is available at the NIH Common Fund’s Data Repository and Coordinating Center (supported by NIH grant, U01-DK097430) website, http://www.metabolomicsworkbench.org), where it has been assigned a Metabolomics Workbench Project ID: (UrqK1489). The data is directly accessible at: http://dev.metabolomicsworkbench.org:22222/data/DRCCMetadata.php?Mode=Study&DataMode=AllData&StudyID=ST000503&StudyType=NMR&ResultType=2&access=UrqK1489#DataTabs

Metabolic pathway construction

Metabolic pathways were constructed according to pathway analysis on Metpa (freely available at http://metpa.metabolomics.ca). Metpa high-level functional analysis is organism specific. Because the metabolite pathways for Equus caballus do not currently exist, the enrichment analysis was performed using Homo sapiens metabolite sets.

The MetPA’s pathway’s topological analysis is based on the centrality measures of a metabolite in a given metabolic network [54]. Centrality is a local quantitative measure of the position of a node relative to the other nodes, and is often used to estimate a node’s relative importance or role in network organization [54]. Since metabolic networks are directed graphs, MetPA uses relative betweenness centrality and out degree centrality measures to calculate compound importance. The impact of a pathway is calculated as the sum of the importance measures of the matched metabolites normalized by the sum of the importance measures of all metabolites in each pathway [54].

RNA isolation, microarray experiments and mRNA and miRNA statistical analysis

Transcriptome and miRNome profiling were performed using Agilent microarrays as described earlier [8]. All the transcriptome and miRNome pre-processing, normalization and statistical analysis steps were carried out as described elsewhere [8]. In contrast to the Mach et al. [8] study, the p-values were corrected for multiple testing using the Benjamini and Hochberg method with a threshold of adjusted p < 0.05, as a compromise between the unadjusted analysis and the Bonferroni procedures.

The data sets supporting the conclusions of this article are available in the Gene Expression Omnibus (GEO) repository, [GSE73104; (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE73104].

Regulatory network of metabolites, metabolic genes and miRNAs: enrichment analysis

The underlying assumption behind the enrichment analysis is that by combining the evidence based on changes in both gene expression and metabolite concentrations, one is more likely to pinpoint the metabolites and pathways involved in the underlying biological processes during endurance exercise. The power of this approach is that such models can provide non-intuitive metabolic and physiological hypothesis [55].

First of all, we converted the equine Ensembl gene IDs to their orthologous associated Homo sapiens gene IDs through the Biomart retrieval tool (Ensembl release 83; http://www.ensembl.org/biomart/).

Then after, for each metabolite with a loading different from zero along any OPLS principal component axis, we identified all the metabolic pathways that the given metabolite participates in through the KEGGREST package in R (Homo sapiens, [hsa] organism). We retrieved the list of genes that encode proteins that participate in the reaction steps around the given metabolic pathway (using the KeggGet function from the KEGGREST package in R). The enzyme-encoding genes whose protein products participate in the metabolic reaction of a given metabolite are referred to here as metabolic genes. We considered metabolic genes as targets of a specific metabolite when they were operating in the metabolic pathways that the metabolite participates in. Each metabolite may be associated with multiple metabolic pathways, and each metabolic gene may be involved in several metabolic pathways.

To investigate the metabolites, transcripts and pathways with more biological relevance during the endurance exercise, we used two complementary statistical approaches. For each of the metabolites presenting at least one metabolic pathway with one coding metabolic gene, we performed: (i) the hypergeometric test for enrichment in target metabolic genes correlated with the metabolite; and (ii) the generalization of the hypergeometric test. Both methods were implemented on the experimental set of ten animals. For the hypergeometric test, we performed a pair-wise Pearson correlation analysis of the expression levels of DEGs and the expression levels of metabolites presenting at least one target metabolic gene. Subsequently, we subtracted the set of correlated metabolites-DEGs with r | > 0.5|; with a p < 0.05. We next looked at whether the metabolic target genes were over-represented when compared with the other genes in the transcriptome. Benjamini and Hochberg correction [56] for multiple testing was applied to the p-values obtained (false discovery rate (FDR) < 0.10).

Moreover, in order to address the hypergeometric test’s lack of power when the numbers of metabolic target genes for a given metabolite was very small, we also implemented a variant of the hypergeometric test. In this variant, we tested for enrichment in metabolic target genes by selecting DEGs regardless of the sign and value of their correlation with the metabolite. When considering a given metabolite with at least one metabolic target gene, we looked at whether their metabolic target genes presented smaller p-values than the other genes in the transcriptome using a one-sided Wilcoxon rank sum test (implemented with the “wilcox.test” function in R). The FDR was determined to correct for multiple testing. Lastly, only metabolites that were significant either in the hypergeometric test or its generalization after correction for multiple testing (FDR < 0.1) were analyzed further. The FDR was set to 0.1 in order to retain as many biological functions as possible.

Identification of miRNAs involved in the regulation of the identified metabolite-metabolic target genes regulatory network

The multiMiR package in R was used to identify potential miRNAs regulating the metabolic target genes included in the regulatory network [57]. For each miRNA expressed in our data set, we assembled a comprehensive list of all experimentally validated target genes using the multiMiR package in R. multiMiR is a comprehensive collection of predicted and validated miRNA–target interactions and their associations with diseases and drugs [58]. It contains human and murine data from 14 external databases that are categorized into three components, including three validated miRNA–target databases (miRecords, miRTarBase and TarBase), eight predicted miRNA–target databases (DIANA-microT, ElMMo, MicroCosm, miRanda, miRDB, PicTar, PITA and TargetScan), and three disease- or drug-related miRNA databases (miR2Disease, Pharmaco-miR and PhenomiR) [58]. Subsequently, we subtracted the set of miRNAs presenting at least one metabolic target gene within our regulatory network.

Regulatory impact factors analysis to unravel the transcription regulation within the regulatory network

The RIF metric [28, 29] was used to identify transcription factors and miRNAs regulators within the set of DEGs. RIF method contrasts two conditions (e.g., case vs. control) and provides a metric to each regulator considering the change in co-expression between the regulator and DEGs. In this study, we focused on genes whose expression was significantly altered at T1 relative to T0, with an adjusted p < 0.05, and we used them as putative regulators for all transcription factors reported by Vaquerizas et al. [59] and the 362 miRNAs expressed in the custom microarray (Additional file 14). The method estimated two RIF alternative measures (RIF1 and RIF2). RIF1 provides information of those regulators that are consistently more differentially co-expressed with the highly abundant and the highly differentially expressed genes. RIF2 scores each regulator considering the most altered ability to act as predictors of the abundance of DEGs.

Multiple factor analysis in an independent cohort of 31 horses

The MFA [60] was applied as an exploratory analysis of the metabolome, transcriptome, and miRNome using FactoMineR [61] R package. MFA was applied in our independent cohort of animals (n = 31) to model complex biological interaction in a holistic manner and identify potential biomarkers with a functional readout of the cellular state. MFA is an extension of PCA tailored to handle multiple data tables that measure sets of variables collected on the same individuals. MFA proceeds in three steps: first, it computes a PCA of each subset of variables and ‘normalizes’ each data table by dividing all its elements by the first singular value obtained from its PCA. In other words, it weights each data table to account from different variances among groups. Second, all the normalized data tables are aggregated into a grand data table that is analyzed via a (non-normalized) PCA that gives a set of factor scores for the observations and loadings for the variables. Third, in order to identify the contribution of each data set to the total variance, that is, to what extent each data set deviates or agrees with what the majority of data sets support, MFA projects a superimposed representation of individuals with each group of data and its barycenter. MFA also provides a measurement of similarity between the geometrical representations derived from each group of variables though a RV coefficient, which is a multivariate generalization of the Pearson correlation coefficient. For each pair of “omic” data sets, the RV-coefficient is calculated as the total co-inertia (sum of eigenvalues of co-inertia, i.e. sum of eigenvalues of the product of two cross product matrices) divided by the square root of the product of the squared total inertia (sum of the eigenvalues) from the individual analysis. As the co-structure between two data sets increases, the RV score move towards to 1. A zero RV score indicates no similarity.