Circulating miRNAs as Putative Biomarkers of Exercise Adaptation in Endurance Horses

Endurance exercise induces metabolic adaptations and has recently been reported associated with the modulation of a particular class of small noncoding RNAs, microRNAs, that act as post-transcriptional regulators of gene expression. Released into body fluids, they termed circulating miRNAs, and they have been recognized as more effective and accurate biomarkers than classical serum markers. This study examined serum profile of miRNAs through massive parallel sequencing in response to prolonged endurance exercise in samples obtained from four competitive Arabian horses before and 2 h after the end of competition. MicroRNA identification, differential gene expression (DGE) analysis and a protein-protein interaction (PPI) network showing significantly enriched pathways of target gene clusters, were assessed and explored. Our results show modulation of more than 100 miRNAs probably arising from tissues involved in exercise responses and indicating the modulation of correlated processes as muscle remodeling, immune and inflammatory responses. Circulating miRNA high-throughput sequencing is a promising approach for sports medicine for the discovery of putative biomarkers for predicting risks related to prolonged activity and monitoring metabolic adaptations.

MicroRNAs are small noncoding RNAs (18-22 nucleotides) that act as post-transcriptional regulators of gene expression, modulating numerous biological processes, including cellular proliferation, differentiation and metabolism (Aoi and Sakuma, 2014). Their genomic context spans introns of protein coding genes, introns and exons of noncoding RNA genes and occasionally exons of protein coding genes. Their binding sites are usually in the 3 ′ UTR and rarely in the 5 ′ UTR or coding region of mRNA sequences (Finnegan and Pasquinelli, 2013). Each miRNA can regulate the expression of hundreds of genes and proteins, while a single gene may be regulated by several miRNAs (Friedman et al., 2009). They can be also released into body fluid in stable forms-termed circulating miRNAs (ci-miRNAs)-that are associated with different RNA-binding proteins or lipoprotein complexes, are actively incorporated into microvesicles and apoptotic bodies and are passively released as the result of cell death (apoptosis or necrosis), autophagy and cell rupture (Xu et al., 2015;Ghai and Wang, 2016). Circulating miRNAs function in transcriptional gene silencing but also act like hormones, in signal transduction by binding Toll-like receptors, (Fabbri et al., 2013).
This specific and controlled cell discharge indicates that released miRNAs are involved in the regulation of pathophysiological mechanisms and that ci-miRNAs are an evolutionarily conserved tool for cell-cell communication, providing fast local or systemic responses to perturbations in homeostasis (Kosaka et al., 2013). On the basis of these findings, ci-miRNAs have been proposed as biomarkers for the diagnosis and prognosis of diseases, as well as for monitoring therapeutic treatments (Ghai and Wang, 2016).
Circulating-miRNAs in the bloodstream could be also used to evaluate the response to training protocols and as molecular predictors of exercise performance in human and equine athletes (Keller et al., 2011;Lombardi et al., 2016;Mach et al., 2016;Polakovičová et al., 2016) because ci-miRNAs plasma levels depend on the type and duration of physical exercise, on the training level and the athlete's nutrition (Kirby and McCarthy, 2013;Makarova et al., 2014;van der Kolk et al., 2015).They also provide greater accuracy than classical serum/plasma biomarkers such as the proteins creatine kinase or lactate dehydrogenase (Makarova et al., 2014) as they are easily affected by physical stress, while those proteins, although consistent markers of muscle metabolism, are partially involved in the cellular response to environmental stimuli (Baggish et al., 2014). In addition, ci-miRNAs are highly stable and are resistant to variations in temperature and pH and to multiple freeze/thaw cycles, making sample storage, and handling easy (Polakovičová et al., 2016).
Various miRNA profiling methods can be used (e.g., qRT-PCR, microarray, sequencing), and they differ in terms of robustness, throughput, accuracy, sensitivity, dynamic range, cost and complexity. The major advantages of miRNAsequencing with next-generation sequencing (NGS) are the ability to detect both novel and known miRNAs, the ability to perform relative quantification of miRNA transcripts present at a low abundance and the fact that miRNA variant characterization can usually be performed in a single experiment, which is especially useful in early biomarker discovery efforts (Ghai and Wang, 2016;Lombardi et al., 2016). To date, this approach has not been applied to profiling ci-miRNA fluctuations in response to endurance exercise.
The aim of this study was to assess the serum profile of ci-miRNAs in response to prolonged endurance exercise in samples obtained from competitive horses through massive parallel sequencing. Tracking the extent of ci-miRNA modulation could improve our understanding of the molecular mechanisms involved in the whole-body adaptive response to endurance exercise and provide useful potential biomarkers to predict disease risks related to prolonged activity (i.e., poor performance), as well as track metabolic adaptations to ultimately establish efficient training programmes.

Ethics Statement
This study protocol has been reviewed and approved by institutional Ethics Committee of the University of Perugia (license No 2016-11). All procedures were performed in accordance with the approved guidelines. Informed consent from the owners and the approval of the Ground Jury President and the Veterinary Commission President of the event were obtained before initiation of any procedures involving the animals.

Animals
Four trained pure-breed Arabian horses (3 females and one gelding reared in Sardinia, Italy) taking part in the same 90 km endurance competition were enrolled. Three horses were 6 years old, and one was 9 years old. All horses were stabled and trained in the same training stable and were subjected to the same management practices at the stable and throughout the ride.

Training Programmes and Competition
The training activities before the competition included exercise with a horse walker every day flat work at lunge and/or in the arena, a training session of variable intensity every other day and a training session at the racetrack. Before the start of the endurance competition, the horses passed the compulsory examination performed according to the rules of the Italian Equestrian Federation (FISE). Animals were fed 1.5 h before the start of the competition with hay and concentrate pellets. In accordance with the FEI and FISE rules, all horses underwent veterinary inspections during the competition every 20-40 km, followed by 40-min periods during which the animals were provided with water, hay and a small amount of concentrate pellets. The mean recovery time (time from arrival at the entrance of the veterinary gate for intermediate veterinary inspections) was 125 s and ranged from 115 to 145 s. All horses successfully completed the race without complication; the average race speed was 17.221 km/h (range: 15.777-21.361 km/h).

Plasma Samples
Peripheral blood samples were obtained by jugular venipuncture at rest just before exercise (T0) within 2 h after the end of an acute exercise bout (T1). Blood was collected in vacutainers containing EDTA (two tube) or a serum clot activator (one tube). One EDTA tube was used to measure haematocrit (HCT), hemoglobin (HGB) concentration, red blood cell (RBC) count, and white blood cell (WBC) count by using the hematology analyser (EosBIO, Italy). The second one EDTA tube was centrifuged at 3,000 × g for 15 min to obtain plasma, which was then stored at −80 • C until ci-miRNA analysis. Blood without anticoagulant was left to coagulate and was then centrifuged (10 min at 3,000 × g); the serum obtained was stored at −80 • C until being used for creatine phosphokinase (CPK), lactate dehydrogenase (LDH), HGB, and total protein measurements. These tests were assessed by a Konelab 2001 (Sclavo, Italy) biochemical analyser using specific kits (Sentinel Diagnostics, Milan, Italy) according to the manufacturer's instructions.

Library Preparation and Next-Generation Sequencing
A total of 6 µl of extracted RNA was converted into microRNA NGS libraries using a NEBNEXT library generation kit (New England Biolabs Inc.) according to the manufacturer's instructions. Each individual RNA sample had adaptors ligated to its 3 ′ and 5 ′ ends and was reversely transcribed into cDNA. Then, the cDNA was pre-amplified with specific primers containing sample-specific indexes. After an 18-cycle pre-PCR, the libraries were purified on QiaQuick columns, and the insert efficiency was evaluated using a Bioanalyzer 2100 instrument on highsensitivity DNA chips (Agilent Inc.,). The microRNA cDNA libraries were size fractionated on a LabChip XT (Caliper Inc.,), and bands representing adaptors and 15-40 bp inserts were excised following the manufacturer's instructions. Samples were then quantified using qPCR and concentration standards. Based on the quality of the inserts and the concentration measurements, the libraries were pooled in equimolar concentrations (all concentrations of libraries to be pooled were the same). Finally, the library pools were quantified again using qPCR, and optimal concentrations of the library pools were used to generate clusters on the surface of a flow cell before sequencing using the v3 sequencing methodology according to the manufacturer's instructions (Illumina Inc.,). Samples were sequenced with an Illumina NextSeq 500 system, producing 50-nucleotide singleend reads.

Raw Data Bioinformatic Analysis
Following sequencing, we performed a data quality control step using FastQC (http://www.bioinformatics.babraham.ac.uk/ projects/fastqc). Mapping against equcab2 (Equus caballus assembly ver. 2) was carried out using Bowtie2 ver 2.2.5 (Langmead and Salzberg, 2012) in a multistep process to remove spike-in reads. The classification of reads is depicted in Figure 1 (data available in Table 2). Annotation of mapped regions was performed using reference annotations (miRBase 21).
After mapping the data and counting the relevant entries in miRBase 21, the number of known microRNAs was calculated. Comparisons of two groups considered all the data collected.

Differential Gene Expression Analyses
The differential gene expression (DGE) analysis was carried out using the edgeR statistical software package (Robinson et al., 2010) comparing the two groups applying a GLM procedure that can be seen as a generalization of a paired t-test. Normalization was carried out using the trimmed mean of M-values method based on log2-fold changes and absolute gene-wise changes Frontiers in Physiology | www.frontiersin.org in expression levels between samples (TMM normalization). Expression levels were measured as tags per million ("TPM"), and p-values for significantly differentially expressed microRNAs were estimated using an exact test on the negative binomial distribution.

Target Genes and Enrichment Analysis
For target analyses, we selected 12 miRNAs, the 6 most downregulated and the 6 most up-regulated miRNAs according to the DGE analysis from edgeR. Human orthologue miRNAs were retrieved using the miRBase software association exchanging tool, and non-matching entries were classified using the BLAST algorithm and microRNAviewer software (http://people.csail. mit.edu/akiezun/microRNAviewer/index.html).
Putative genes (predicted and/or validated) targeted by these 12 miRNAs were identified using human miRNA IDs in the miRWalk webtool (Dweep and Gretz, 2015). Only genes included in the seven available databases (miRWalk, PITA, MicroT4, RNA22, miRanda, RNAhybrid and TargetScan) were retained in the analysis. For miRNAs without gene targets, we used a validated list of genes downloaded from the latest database update (DIANA Tarbase V. 7.0, Vlachos et al., 2015). Validated gene targets were also downloaded from the miRWalk 2.0 website. A unique list that merged both predicted and validated targets was produced for each miRNA, and common genes between different lists were identified.
The Cytoscape suite (Shannon et al., 2003) was used to build a protein-protein interaction (PPI) network and the PSIQUIC service (http://www.ebi.ac.uk/Tools/webservices/psicquic/view/ home.xhtml) was used to download proteins interactions from six different public databases (ChEMBL, Spike, GeneMania, mentha, MINT, IntAct). These six networks were merged together in order to generate a unique PPI network with all the interactions. Only Homo sapiens interactions were used to build the new network, to avoid loss of information using equine data only.
The application clusterMaker 2.0 (Morris et al., 2011) with the "gLay" option was used to find different clusters inside the network according to the number of connections between nodes. Clusters were explored for their GO biological process enrichment terms using the application BiNGO, filtering for results with a corrected q-value lower than 0.05 (Benjamini Hochberg correction).
An overview of the experimental design and applied methods is depicted in Figure 2.

Hematological and Biochemical Parameters
All horses included in the study completed a race and successfully passed the final veterinary examination at the veterinary gate. Changes in average hematological and biochemical parameters measured in blood samples taken before and after the race ( Table 1) reflect the physical and metabolic effort that occurred during the endurance competition.

Sequencing Results
The data were quality checked for intensity correction, base calling and the assigning of Q-scores: the samples showed an overall high quality, with the vast majority of the sequences obtained presenting a Q-score higher than Q30. Raw data are available through BioProject with the following accession ID: PRJNA434539.
In the spike-in quality control step we observed an excellent correlation between the samples for counts corresponding to the spike-in values (R 2 > 0.99, see Supplementary Figure 1).
On average, 29 million reads were obtained per sample, and the lowest number of reads for any sample was 12.9 million. Summary statistics for reads and mapping are detailed in Table 2. After mapping the sequences and counting occurrences of relevant entries in miRBase 21, the number of known microRNAs was calculated (

Differential Gene Expression Analysis
After statistical analysis with edgeR, 102 entries (miRNA and smallRNA) were found to be differentially expressed at T1 vs. T0 at a significance of q < 0.05 and an absolute log2-fold change (logFC) >1 or <-1. Using these filters, 27 RNAs were identified as being up-regulated (logFC >1), and 75 were identified as being down-regulated (Supplementary Table 1). Table 4 shows the individual results for the 20 most differentially expressed microRNAs (10 up-regulated and 10 down-regulated). Values are expressed as the means ± standard deviation.

Target Genes and Enrichment Analysis
A subset of 12 miRNAs was created by choosing the 6 most upregulated and the 6 most down-regulated miRNAs based on the false discovery rate threshold to produce a list of unique target genes using both predicted and validated targets (Supplementary Table 2 for up-regulated and Supplementary Table 3 for downregulated). Then, common target genes were identified: only those targeted by 4 or more miRNAs were selected ( Table 5).
Another filter based on the number of interactions was applied: only entries with more than five connections were further evaluated. These genes were used to create a PPI network as explained in the materials and methods section to identify clusters of proteins within the six merged network (merged in order to generate a unique PPI network with all the interactions) considering connections between nodes as selection criteria. This was done for both up-and down-regulated miRNA-derived targets. The network built with up regulated genes contains 2621 nodes and 8600 edges, while the network built with down regulated genes contains 861 nodes and 1535 edges (Supplementary Table 5). The major protein clusters are reported in Figures 3, 4.
Gene Ontology (GO) enrichment analysis with BiNGO was performed on the clusters harboring one of the targets or where a central node was targeted by several miRNAs and recognized as relevant in our biological framework ( Table 6, and  Supplementary Table 4).

DISCUSSION
This study profiles the amounts and types of ci-miRNAs in horse athletes competing in endurance races at rest and within 2 h after the end of competition, when the majority of the significant changes in ci-miRNAs occur (Nielsen et al., 2014). We identified the modulation of a large set of miRNAs involved in the response to exercise and inferred related target genes clustered in PPI networks that are correlated with different tissues (e.g., skeletal muscle, heart, liver and blood) and processes, such as muscle remodeling, metabolic changes and the modulation of immune and inflammatory responses.
The changes in ci-miRNAs occurring in our athletes seemed to be a result of active or selective secretion rather than muscle cell damage. This is supported by increases in ci-miR-1, ci-133a, and ci-206, which are likely correlated with performance parameters and are poorly associated with biochemical markers of cardiac and/or skeletal muscle damage (Mooren et al., 2014), and by the decrease in the classical myomiR ci-miR-486 (Aoi et al., 2013;Xu et al., 2015).
Both adaptations to exercise and muscle regeneration involve myomiRs (Diniz and Wang, 2016) which act as modulators of myogenesis, mitochondrial biogenesis, hypertrophy and energy metabolism. For example miR-1 acts in the development and differentiation of smooth and skeletal muscle tissue together with miR-206 and miR-133. Overexpression of miR-1 or miR-206 promotes myogenic differentiation and the repair of exerciseinduced muscle injury (Wu, 2016) and, along with the upregulation of myomiR-133a, regulates mitochondrial biogenesis via the IGF-1-Akt (serine/threonine kinase 1) pathway through the IGF1R. In support of this, IGF1R was identified as a central node of cluster 17 (CNC 17) in our network analyses (Table 6, Figure 3). The down-regulation of mir-504 confirmed the presence of mitochondrial adaptation, which is essential in skeletal muscles involved in aerobic exercise. Namely, miR-504 acts as a negative regulator of p53, that has been implicated in regulating substrate metabolism and in exercise-induced mitochondrial biogenesis in skeletal muscle (Hu et al., 2010; Modulation is expressed in log2-fold change (logFC). Entries are sorted by FDR value. Bartlett et al., 2014). These findings highlight the essential role of the previously described miRNAs in exercise tolerance and the response to training. Cardiac muscle adaptations could involve EGFR (CNC 1, Table 6, Figure 3), plays roles in cell growth and proliferation, cardiac development, cardiac responses (adaptation/remodeling) to physiological and pathological loads. EGFR also modulates cardiomyocyte stress tolerance through the MAP kinase and phosphatidylinositol 3 kinase (PI3K)/Akt pathways (Reichelt et al., 2017). MicroRNAs regulate cardiac responses by targeting key components of these pathways. In our experiment, the upregulation of mir-1 and mir-133 indicates the inhibition of pathological cardiac hypertrophy via the suppression of target genes in the heart and neural crest . Moreover, these two miRNAs appear to be particularly responsive to aerobic exercise and positively associated with VO 2 max (Mooren et al., 2014) and are associated with an increase in whole blood content in endurance athletes and with superior cardiorespiratory fitness levels (Clauss et al., 2016;Denham and Prestes, 2016).
The occurrence of general anti-hypertrophic conditions in our athletes during endurance racing is also suggested by the decrease in miR-328 associated with atrial fibrillation and cardiac hypertrophy (Li et al., 2014). In addition, up-regulation of myomiRs, as shown in our experiment, leads to an antihypertrophic signal (Soci et al., 2016) throughout regulation of PURB. PURB (CNC 9, Table 6, Figure 3) acts on the myosin genes and is considered a repressor of slow muscle contractions (Soci et al., 2016) that acts by controlling vascular smooth muscle alpha-actin gene transcription and repressing myoblasts (Hariharan et al., 2014)..
Additionally miR-199 and the miR-99/100 family, which we found to be up-regulated and which are not strictly considered to be myomiRs, play a role in the maintenance of cardiac homeostasis and in regulating apoptosis in cardiomyocytes through their target, IGF1R (CNC 17, Table 6, Figure 3; Chen et al., 2015;Li et al., 2016). Possible adaptive changes in cardiac tissue are also suggested by the increases in miR-499 and miR-208b, both of which are expressed in skeletal and cardiac muscles (Soci et al., 2016). These miRNAs, together with miR-206, act in the gene reprogramming process, activating slow myofibres and repressing fast ones, resulting in greater muscular endurance capacity (Kirby and McCarthy, 2013) . Muscle tissue remodeling also includes an increase in sympathetic and parasympathetic innervation involving muscle and cardiac autonomic nerve remodeling along with electrophysiological changes in the heart through the overexpression of miR-206 (Jeng et al., 2009;Zhang et al., 2015).
Structural modifications in muscles also occur via the regulation of stress responses and thyroid hormone signaling pathways involved in the modulation of myosin heavy chain genes thorough the overexpression of miR-208a (Baldwin and Haddad, 2001;Grueter et al., 2012). However, target and pathway analyses suggest that morpho-functional modifications of smooth, cardiac and skeletal muscles involve not only contractile proteins but also cytoskeletal proteins. Indeed,  TAGLN, TMOD, and LASP1 (CNC 7, Table 6, Figure 3) act in stabilizing the contractile proteins and in regulating the dynamic actin cytoskeleton. Transgelin, is down-regulated by the cellular transformation involved in major cytoskeletal rearrangements, and it is implicated in smooth and cardiac muscle cell differentiation (Doll et al., 2007;Navickas et al., 2016). Tropomoduline, working in conjunction with tropomyosins, regulates the distribution and length of thin filaments in the sarcomeres of cardiac and skeletal muscle cells (Colpan et al., 2013;Rao et al., 2014). Moreover TMOD down-regulation, which occurred in our athletes, may contribute to the optimization of sarcomere length to maximize muscle activation and generate greater active forces (Kolb et al., 2016). Possible modifications of the sarcolemma may be achieved by SGCD (CNC, 13 Table 6, Figure 3), a protein that forms sarcoglycan complexes. These complexes play a key role in the mechanical stabilization of the sarcolemma, and they are modified in response to muscle contractions (Tarakci and Berger, 2016).
These results clarify the mechanisms underlying the massive structural and metabolic changes that generally occur in endurance athletes, which include an increase in mitochondria, fiber-type shifts and changes in capillary density (Kirby and McCarthy, 2013). Aerobic exercise training up-regulates miRNA expression to promote myoblast differentiation and proliferation and enables the metabolic reprogramming of skeletal and cardiac muscle cells to promote adaptations for improved physical performance (Denham and Prestes, 2016).
These anatomic changes all result in an increased ability to sustain prolonged physical activity, which also requires a greater energy supply. The reduction in ci-miR-486, which is proposed as a biomarker of advantageous adaptations to prolonged exercise, could be related to its specific uptake from blood into muscle cells (Aoi et al., 2013;Xu et al., 2015). The consequently high levels of miR-486 in these cells could result in enhanced muscle glucose uptake since the main target of this miRNA is phosphatase and tensin homolog (PTEN), a negative modulator of insulin signaling. Moreover, miR-486 can influence muscle protein turnover, acting as a modulator of the transcription factor FOXO (Xu et al., 2012(Xu et al., , 2015Kirby and McCarthy, 2013). However, independent of their specific roles and targets, the up-regulation of myomiRs in skeletal muscles after endurance exercise may be related to the decrease in muscle protein synthesis described in the post-exercise period (Rennie and Tipton, 2000), which may help to maintain cellular energy levels. Moreover, the involvement of SLC7A1 (CNC 12, Table 6, Figure 3), a non-hepatic trasponder of arginine, lysine and ornithine, indicates that amino acid uptake could also be modulated following endurance exercise.
The modulation of energy metabolism is also indicated by the involvement of SLC5a3, STX6, and IGFR1 (CNC 6 and 17, Table 6, Figure 3). The presence of SLC5a3, has recently been observed in cardiomyocytes, where it seems to act as a glucose sensor (Van Steenbergen et al., 2017), and STX6 can modify glucose uptake modulating GLUT4 translocation (Foley and Klip, Figure 3), which act as Na+/glutamate and cysteine/glutamate transporters, respectively, suggests that exercise could regulate glutamate transport across membranes, as has been proposed by other authors, possibly to alleviate the neurotoxicity of excessive glutamate release (Yang et al., 2012).
Collectively, our results show activation of immune and inflammatory responses via the modulation of miR-224 and miR-1180 through a key immune systems regulator, transcription factor κB (NF-κB) (Yang and Wang, 2016). Indeed, inflammatory stimuli, such as strenuous exercise, can induce activation of the innate immune system via miR-224 and the pentraxin 3 (Ptx3) gene, both mediators of inflammation; miR-224 is a transcriptional target of NF-κB, and the Ptx3 promoter contains an NF-κB binding site for its regulation (Rudnicki et al., 2014). The up-regulation of miR-224 leads to activation of EGF (its receptor is CNC 1, Table 6, Figure 3), transforming growth factor alpha (TGF-a), IGF (its receptor is CNC 17 Table 6, Figure 3) and tumor necrosis factor α (TNFα), inducing immune cell proliferation and migration via inflammatory response signals (Scisciani et al., 2012).
The down-regulation of miR-1180 suggests a modulation of immune cell survival and proliferation by directly targeting key inhibitors of the NF-κB signaling pathway and the apoptotic protein Bcl-2 (CNC 8, Table 6, Figure 3) (Tan et al., 2016). Moreover, PTMA (CNC 4, Table 6, Figure 3) acts in the immune response as an extracellular signaling molecule via Toll-like receptor 4 (Mosoian, 2011), decreasing the apoptotic response and increasing cell survival through activation of the serinethreonine kinase Akt (Cannavo et al., 2013).

CONCLUSIONS
The modulation of protein coding genes in response to physical activity has been widely investigated, but a new level of regulation mediated by noncoding RNAs has only recently emerged as a key player in post-transcriptional regulation that may potentially be involved in adaptations in response to exercise.
We identified a large set of modulated ci-miRNAs and target genes that clustered in PPI networks, revealing that responses to endurance exercise induce major changes in muscle remodeling, metabolic pathways and, ultimately, in the immune system and the inflammatory response mechanisms. Therefore, modulated ci-miRNAs in endurance athletes can be considered promising and reliable biomarkers of stress and/or training.
Moreover, high-throughput ci-miRNA sequencing is a promising approach for use in sports medicine. The discovery of putative biomarkers for use in predicting disease risks related to prolonged activity (i.e., poor performance) and monitoring metabolic adaptations to establish efficient training schemes could be transferred to athlete of all species, including humans.

ACKNOWLEDGMENTS
The authors are grateful to Mr. Gianluca Alunni for his excellent technical assistance with sample preparation.