Integrated analysis of microRNA and mRNA expressions in peripheral blood leukocytes of Warmblood horses before and after exercise.

Exercise capacity is a valuable trait in horses, and it has been used as a horse selection criterion. Although exercise affects molecular homeostasis and adaptation in horses, the mechanisms underlying these effects are not fully described. This study was carried out to identify changes in the blood profiles of microRNAs (miRNAs) and mRNAs induced by exercise in horse leukocytes. Total RNAs isolated from the peripheral blood leukocytes of four Warmblood horses before and after exercise were subjected to next-generation sequencing (NGS) and microarray analyses to determine the miRNA and mRNA expression profiles, respectively. The expressions of 6 miRNAs, including 4 known and 2 novel miRNAs, were altered by exercise. The predicted target genes of the differentially expressed miRNAs identified by NGS were matched to the exercise-induced mRNAs determined by microarray analysis. Five genes (LOC100050849, LOC100054517, KHDRBS3, LOC100053996, and LOC100062720) from the microarray analysis were matched to the predicted target genes of the 6 miRNAs. The subset of mRNAs and miRNAs affected by exercise in peripheral blood leukocytes may be useful in elucidating the molecular mechanisms of exercise-associated physiology in horses.

Introduction

Horses have been domesticated for working, riding, recreation, sport, and racing, and the development of specific breeds has resulted in the selection of athletic phenotypes [21]. Exercise and performance are important factors in the selection of horses because the value of the domesticated animal is determined by its abilities to show improved performance [16]. The Warmblood breed, which has been used for dressage and show jumping, was selected based on its good athletic performance [16]. It has generally been accepted that equine genetics and genomics have a critical role in determination of the horse's capabilities [2]. Many studies have tried to identify the candidate performance genes related to genetic contribution or heritability of various equine performance traits [2,12,14].

Exercise and training induce molecular homeostasis and adaptations controlled by the transcriptional and translational regulation of genes that encode proteins [34]. Short- and long-term exercise in athletic horses has been shown to cause immediate molecular changes in a wide range of mRNA transcripts [32631]. Genes underlying an adequate stress response and exercise-adapted phenotypes include those related to energy metabolism, muscle structure and function, hypoxia, and fatty acid oxidation [1435]. Appropriate gene expression and regulation have been believed to have important roles in the determination of most exercise-related phenotypes in different equine breeds [2331]. However, little research has been conducted on the non-genetic regulatory mechanisms of exercise-related genes.

There has been an increase in studies of small RNAs (sRNAs), including the gene expression and regulation of microRNAs (miRNAs), which are approximately 22 nucleotides (nt) long, and short non-coding RNAs (ncRNAs) that have important roles in regulating target gene expression through mRNA degradation and translational inhibition [11]. In horses, miRNAs have been identified in equine major organs and blood, and they are considered to have key functions in equine physiology by regulating gene expression [1720]. Performance studies for horses have used microarray technology to determine the genes associated with exercise [26]. Recently, next-generation sequencing (NGS) has emerged as an important tool when scrutinizing sRNAs, including miRNAs. Because of its ability to generate millions of reads of predetermined lengths, NGS greatly improves the capacity to detect novel miRNAs on a genomic scale [17].

In humans, brief exercise affects the expression levels of miRNA and/or mRNA in blood cells [36]. Blood collection and a subsequent genomics approach could provide a safe and rapid method to investigate exercise-related gene expression profiles in horses [3]. Furthermore, it has been reported that gene expression changes in peripheral blood cells could reflect physiological adaptations in horses [3]. This study was performed to investigate the effect of exercise on the expression of miRNAs and mRNAs in the leukocytes of healthy Warmblood horses. The present study could be valuable in identifying the potential mechanisms involved in the molecular homeostasis and adaptation of horses and in providing information for application during horse selection and breeding.

Materials and Methods

Experimental animals

The present study included four horses belonging to the Seoul Horse Race Association. Informed consent was obtained from the Seoul Horse Race Association. All the experimental protocols were reviewed and approved by the Institutional Animal Care and Use Committee of Seoul National University (approval No. SNU-131218-1). The horses were of the German Warmblood type and were four geldings ranging from 7 to 15 years old (mean age, 11 years). The body weights of horses ranged from 450 to 645 kg (average 551 kg). The horses received pelleted concentrates of mixed grains (oats and corn) with timothy grass and alfalfa hay. Diets were provided to the horses at a medium intake level (1.5 kg/feeding three times each day) with free access to water. Riding and exercise training were regularly performed every other day during daytime; twice per day for approximately 30 min (usually in the morning, 10:00 am to 10:30 am; and afternoon, 2:00 pm to 2:30 pm) with mild endurance exercise. The horses were managed in individualized stables located in a single building and were closely monitored by resident veterinarians.

Blood collection and leukocyte isolation

Exercise of the horses was performed during the morning (10:00 am to 11:00 am) four hours after feeding. The horses were subjected to controlled exercises as follows: walking for 5 min, trotting for 10 min, cantering for 10 min, trotting for 10 min, cantering for 10 min, and trotting for 10 min, and walking for 5 min. The exercise intensity was higher than that of the horses during the previous exercise and training period. Immediately before and after exercise, 6 mL of blood was sampled from the jugular veins of each horse by using anti-coagulation tubes containing ethylenediaminetetraacetic acid. The samples were centrifuged at 4℃ at 1,000 × g for 10 min, and subsequently, the leukocytes from the buffy-coat area were collected by using sterilized Pasteur pipettes. To minimize the contamination of peripheral blood-derived platelets, the buffy-coat preparations were immediately washed twice times with cold phosphate-buffered saline (Welgene, Korea) at 4℃, followed by centrifugation at 700 × g. The supernatant was removed and the residual pellet was subjected to nucleic acid extraction.

Isolation of total RNA

Total RNAs including miRNAs were isolated from leukocytes by using the miRNeasy extraction kit (Qiagen, USA) following the manufacturer's recommendations. After the purification, the total RNAs were quantified by absorbance at 260 nm, and the integrity of the RNAs was assessed by using a Bioanalyzer 2100 (Agilent Technologies, USA). The high-quality of total RNAs in two horses (No. 1 and 2 horses) was subjected to NGS and microarray analyses, which was performed at the Theragen Bio Institute (Korea). For those total RNAs that were insufficient to perform both NGS and microarray analyses simultaneously due to the low amounts and/or quality of the nucleic acids, we performed NGS analysis for miRNA expression in one horse (No. 3) and transcriptome by microarray analysis in the other horse (No. 4). The mRNAs and miRNAs were analyzed separately, and their results compared. The overall scheme for the integrated analysis of miRNAs and mRNAs via NGS and microarray analysis is presented in Fig. 1.

Fig. 1

Processes in the analysis of the expression of leukocyte microRNAs and mRNAs that were changed by exercise. NGS, next-generation sequencing.

Generation of DNA library and NGS analysis

The generation and preparation of a DNA library, acquisition of raw read counts, and standard bioinformatics analyses were performed as previously described [17]. Briefly, the reads were mapped to the equine genome by using the Short Oligonucleotide Analysis Package software (Beijing Genomics Institute, China) [22] or Bowtie (Johns Hopkins University, USA) [19]. The sRNA reads were annotated using the GenBank database (National Center for Biotechnology Information [NCBI], USA) and the Rfam RNA database (Sanger Institute, UK). Sequences matching known equine sRNAs, including ribosomal RNA (rRNA), small cytoplasmic RNA (scRNA), small nucleolar RNA (snoRNA), small nuclear RNA (snRNA), and transfer RNA (tRNA), were excluded from further analyses. Reads mapped to antisense exons, antisense introns, exons, and introns, and the repeat reads were also excluded. The miRBase 21 database was used to identify known miRNAs in horses; miRNAs perfectly matched sequences in equine genome database were considered as known, and those not annotated to any RNA categories considered as novel. Mireap software (NCBI) was used to investigate the hairpin structures of the remaining miRNA precursors, dicer cleavage site, and the minimum free energy of the novel miRNAs. The names for novel miRNAs were arbitrarily assigned in numerical order following the Theragen Bio Institute's criteria for miRNA nomenclature. The sequencing data of miRNAs were deposited in the NCBI under BioProject ID PRJNA306540.

Differentially expressed miRNAs (DEMs)

The miRNA expression levels in the pre-exercise samples were compared to those in the post-exercise samples by plotting the log2-ratio figure and scatter plot. To compare the expression level of an miRNA across samples, read counts obtained from NGS were normalized as reads per kilobase million (RPKM) [8].

The known and novel DEMs were determined based on the criteria of a greater than 2-fold difference and a p value less than 0.05 between the pre- and post-exercise samples. The predicted targets of the DEMs were identified by using computational tools, such as TargetScan (Whitehead Institute for Biomedical Research, USA) and miRanda (Computational Biology Center of Memorial Sloan-Kettering Cancer Center, USA). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed by using the predicted target genes with the Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics tool (Laboratory of Human Retrovirology and Immunoinformatics, USA) with low functional annotation clustering, a false discovery rate (FDR) < 5, and a p value < 0.05.

Microarray analysis of mRNAs and the identification of differentially expressed genes (DEGs)

Cyanine 3-labeled complementary RNAs (cRNAs) were generated by using the Low RNA Input Linear Amplification Kit (Agilent Technologies) according to the manufacturer's recommendations. The quality and quantity of the labeled cRNAs were determined by using a Nanodrop spectrophotometer (Nanodrop Technologies, USA). The labeled cRNAs were subjected to the one color, 4X44K Horse Gene Expression Microarray Kit (Agilent Technologies) and the Gene Expression Hybridization Kit (Agilent Technologies) for 17 h at 65℃. The hybridized microarray chips were washed with the Gene Expression Wash Buffer Kit (Agilent Technologies).

The microarray chip was scanned by using a DNA microarray scanner (Agilent Technologies), and the raw signal density was acquired by using Feature Extraction software (Agilent Technologies). The threshold raw signal values were set at 1.0, and all raw signal values were normalized by using a percentile shift (75th percentile). The normalized microarray expression data were analyzed by using GeneSpring GX12 software (Agilent Technologies). The DEGs were classified by using Benjamini-Hochberg's FDR method [4]. A list of the DEGs between the pre-exercise and post-exercise samples was compiled by using the criterion of log2 (fold change) ≥ 1 or ≤ −1, and then, the DEGs were subjected to further analyses by using the functional annotation tools in DAVID with a p value < 0.05. Microarray results were deposited in the NCBI under accession number (GSE76310).

Integrated analysis of DEGs and DEMs

To identify the networks and relationship between miRNAs and mRNAs, the DEGs identified by the microarray analysis were matched to their predicted target DEMs as determined by NGS.

Results

NGS analysis of sRNAs

Total RNAs were isolated from the leukocytes of horses before and after exercise, and those RNAs of a high quality and a sufficient amount were subjected to NGS analysis. The total number of read counts (range, 2,511,025–24,070,381) was obtained from the cDNA libraries of sRNAs. Clean reads with more than approximately 85% high quality were subjected to further analyses by using bioinformatics tools. The resulting sRNAs ranged from 18 to 30 nt in length, with the majority having a length of 21 to 23 nt. The comparative results of each sRNA sequence comprised 9% to 20% of common reads, which comprised more than 90% of all reads. The sRNAs were compared with the equine genome database. The unique sRNAs were mapped to the horse genome without repeats, and they comprised 7.36% to 12.25% of the total reads. The rRNA comprised 44.36% to 71.20% of the sRNAs. The total sRNAs with repeats accounted for 54.84% to 57.13% of all reads, and the scRNA comprised 47.75% to 68.69% of the sRNAs. The majority of the sRNAs that were matched with sequences in NCBI GenBank and Rfam databases were ncRNAs, including exons, introns, repeats, rRNA, snRNA, snoRNA, scRNA, srpRNA, and tRNA. The miRNAs from the total sRNAs accounted for 0.59% to 1.21% of the reads without repeats and 52.41% to 55.35% of the total reads.

The number of known miRNAs in all samples was 229. The sequences of the unannotated sRNAs that could be mapped to the equine genome were subjected to further analyses to identify novel miRNA candidates. A total of 150 miRNAs were novel in all samples.

Exercise-induced DEMs

The expression profiles of the known and novel miRNAs in horse leukocytes before and after exercise were analyzed. The following subsets of exercise-specific miRNAs were identified: 4 known miRNAs and 2 novel miRNAs. The upregulated miRNAs after exercise included eca-miR-423-5p. The downregulated miRNAs after exercise were eca-miR-144, eca-miR-33a, and eca-miR-545. Novel miR-14-5p was present before exercise but was not detected after exercise. Novel miR-95-3p was not detected pre-exercise but was detected after exercise (Table 1). The secondary structure and mature sequence of the novel miRNAs are presented in Fig. 2.

Table 1

Differential expression analysis of known and novel miRNAs derived from NGS data

Read counts were normalized to reads per kilobase million. NGS, next-generation sequencing; miRNA, microRNA; FDR, false discovery rate.

Fig. 2

Characteristics of the novel microRNAs (miRNAs) affected by exercise in equine peripheral blood leukocytes. The structure of the two novel miRNAs. (A) Novel_mir_14-5p and (B) novel_mir_95-3p. The box indicates the mature sequence of the novel miRNA.

In silico analysis with programs such as TargetScan and miRanda revealed that the 4 known miRNAs and the 2 novel miRNAs potentially targeted 1,625 genes and 905 genes, respectively. The eca-miR-144 targeted 335 genes, eca-miR-33a targeted 353 genes, eca-miR-423-5p targeted 367 genes, eca-miR-545 targeted 570 genes, novel miR-14-5p targeted 115 genes, and novel miR-95-3p targeted 790 genes.

The KEGG pathway analysis aided in the categorization of the predicted target genes into various functional groups and resulted in a total of 13 pathways from the genes targeted by known miRNAs and 16 pathways from the genes targeted by novel miRNAs. The following significant pathways were identified based on p < 0.05: the Wnt signaling pathway, the transforming growth factor (TGF)-β signaling pathway, and the adherens junction pathway.

Identification of DEGs by microarray analysis

The exercise-related DEGs identified from three horses that exercised for one hour were based on changes in expression levels with a log2 (fold change) ≥ 1 or ≤ −1. A total of 12,422 genes were differentially regulated by exercise. In horses No. 1, 2, and 4, a total of 9,010 and 3,412 genes were differentially upand downregulated, respectively. Among these DEGs, 87 upregulated genes and 27 downregulated genes were commonly identified in those three horses (Fig. 3), and the 9 genes showing statistical significance with the expression of a log2 (fold change) ≥ 2 or ≤ −2 were present in those horses (Table 2). The functions of 114 DEGs were investigated by assessing the pathways provided in the DAVID database (Laboratory of Human Retrovirology and Immunoinformatics, USA). However, no significant pathway was identified.

Fig. 3

Differentially expressed genes from the microarray analysis before and after exercise of three horses (No. 1, 2, and 4). (A) Significantly increased genes based on log2 (fold change) > 1 or < −1. (B) Significantly decreased genes based on log2 (fold change) > 1 or < −1. Each circle in panel A and B indicate the number of genes differentially expressed in the peripheral leukocytes of horse No. 1, 2, and 4, respectively. (C) The number of differentially expressed genes in the three horses.

Table 2

A representative list of commonly up- and downregulated genes from the microarray data for three horses after exercise

There have been 28 candidate performance genes previously reported to be associated with exercise or performance in horses [35]. A total of 10 genes among the candidate 28 genes for physical performance in the horse were identified in the samples of the leukocytes of three Warmblood horses (No. 1, 2, and 4). These include upregulated genes, such as alpha-actinin-3 (ACTN3), insulin-like growth factor 1 receptor (IGF1R), adenosine monophosphate deaminase 1 (AMPD1), dopamine receptor D 1 (DRD1), dopamine receptor D 3 (DRD3), and 5-hydroxytryptamine (serotonin) transporter (5HTT), and downregulated genes, including hypoxia-inducible factor 1-alpha (HIF1A), glycogen synthase 1 (GYS1), and vascular endothelial growth factor A (VEGFA).

The microarray data were integrated to the predicted target genes of specific 6 DEMs identified by NGS. Consequently, four upregulated DEGs were targeted by two downregulated DEMs and two downregulated DEGs were targeted by a single upregulated DEM (Fig. 4). The networks of five mRNAs (LOC100050849, LOC100054517, KHDRBS3, LOC100053996, and LOC100062720) and 3 miRNAs (novel miR-95-3p, eca-miR-545, and eca-miR-144) were identified (Fig. 4).

Fig. 4

Regulation of microRNAs (miRNAs) and mRNAs is altered by exercise. (A) and (B) illustrate the relationships between the differentially expressed mRNAs in the microarray data and the mRNAs targeted by the differentially expressed miRNAs induced by exercise identified in the next-generation sequencing data.

Discussion

This report is the first study analyzing the relationships of mRNA and miRNA with exercise in Warmblood horses. In the analysis of the miRNAs, 229 known miRNAs and 150 novel miRNAs were found to be altered in circulating leukocytes after exercise. Significantly, 4 known miRNAs and 2 novel miRNAs had common changes in expression among 3 horses following one hour of exercise, suggesting these are exercise-specific miRNAs in blood leukocytes. Leukocytes have been reported to contain approximately 12,500 times as much RNA than that in platelets, and the contamination of leukocyte preparations with platelets has been considered to rarely interfere with the extracted RNA profile of leukocytes [10].

The known and novel miRNAs were eca-miR-144, eca-miR-33a, eca-miR-545, eca-miR-423-5p, novel miR-14-5p, and novel miR-95-3p. In the present study, miR-144 and miR-33a matched the exercise-specific miRNAs in animal studies. The miR-144 expression decreased in this study, and it has been observed to increase in the heart muscles of rats after one hour of swimming [23]. The differences between these two studies, such as the species, the circadian systems characteristics of diurnal and nocturnal animals, the small number of samples and its origin, and the exercise type and intensity, may be the cause of their contrasting gene expression results. Additionally, miR-33a, which was downregulated after one hour of exercise in three Warmblood horses, increased after 30 minutes of trotting in thoroughbred horses [13]. This may result from differences in the breed, the type and duration of exercise, and the history and life management previously entrained by exercise or photoperiod [29]. The functions of the novel miRNAs are presently unreported, which warrants further studies to identify the roles of the miRNAs in exercise.

The results of the target prediction for the identified DEMs were analyzed by DAVID. In the KEGG pathway analysis, links to significant pathways, including Wnt signaling pathway, TGF-β signaling pathway, and adherens junction pathway, were detected. Some studies have shown the involvement of these signaling pathways in exercise in animals.

Wnt signaling is recognized as an important regulator of bone, heart, and muscle function and has been reported to represent an exercise-induced signaling mechanism underlying skeletal muscle adaptation in rats [1]. Exercise is found to affect miRNA expression associated with TGF-β signaling [33], and the increased expression of TGF-β1 in tendon tissues after acute exercise is associated with collagen type I synthesis in tendon-related connective tissues [15]. Ventricular TGF-β1 expression has been reported to contribute to myocardial remodeling in physiological cardiac hypertrophy [5]. Regular exercise enhances regulatory T-cell function by increasing regulatory T-cell mediators TGF-β and IL-10 in mononuclear leukocytes [38]. One study reported that brief exercise is associated with adherens junction pathway on peripheral blood natural killer cells [32].

In our mRNA analysis, 87 mRNAs were significantly upregulated, and 27 mRNAs were downregulated by exercise in circulating leukocytes. Because the horses were assessed as clinicopathologically healthy and had been exposed to mild exercise or training, these genes may be associated with exercise adaptation or performance in Warmblood horses rather than other physiological conditions, such as an acute stress response to exercise. A total of 10 genes related to 28 previously reported equine performance genes were identified. Among them, ACTN3, IGF1R, AMPD1, DRD1, DRD3, and 5HTT were upregulated, and HIF1A, GYS1, and VEGFA were downregulated. However, expression of these genes failed to reveal statistical significance with exercise in this study. Interanimal variations and the limited number of horses in the present study may explain the differences between study results.

The ACTN gene, which was upregulated after exercise in this study, is present in fast and glycolytic muscle fibers and has a pivotal role in muscle strength at high speeds [6]. IGF is reported to activate exercise-induced neuronal and cognitive enhancement [9]. Adenosine monophosphate (AMP) deaminase, encoded by AMPD1, is activated in skeletal muscle after AMP accumulation during exercise [30]. The HIF1A gene is upregulated in hypoxia and under severe and repetitive oxygen stress conditions and is involved in red blood cell production, angiogenesis, and glucose metabolism [2437]. In this study, HIF1A gene expression decreased after exercise. The expression of the HIF1A gene in human leukocytes was reduced in training groups, but exercise induced upregulation in the expression of the HIF1A gene in human muscles [27]. Exercise also decreased the expression of the GYS1 gene. Muscle glycogen synthase is an important enzyme for glycogen synthesis in skeletal muscles. GYS1 is a key candidate for diagnosing polysaccharide storage myopathy in horses [25]. The VEGFA gene is involved in endothelial cell proliferation and migration in peripheral circulation. VEGF gene expression of human skeletal muscle in hypoxic or normoxic training groups is upregulated, but the expression of the gene in human leukocytes under the same conditions is downregulated [28]. The DRD1 and DRD3 genes were overexpressed. Dopamine receptors (D) pass signals from one nerve cell to an adjacent cell, and dopaminergic transmission is associated with reward dependence and the persistence of exercise [18]. The serotonin transporter transports serotonin from the synaptic cleft to the presynaptic neuron for removal and recycling, and serotonin biosynthesis in the brain results from exercise to control emotions [7].

The targeted genes from the DEMs of NGS were matched with the genes identified by microarray analysis. Five genes and three miRNAs were found. LOC100050849 is a cyclic AMP-dependent transcription factor ATF-3-like protein. LOC100054517 is a transmembrane protein 163-like protein. KHDRBS3 is a KH domain containing, RNA binding, and signal transduction associated 3-like protein. LOC100053996 is tetratricopeptide repeat 23-like protein. LOC100062720 is a potassium voltage-gated channel subfamily E member 2-like protein. A total of six regulated miRNA-mRNA networks were identified, including mRNAs known to be regulated by two of the miRNAs described. These RNAs may be potentially associated with the regulation of exercise-related physiological processes, such as respiratory, glucose metabolism, molecular adaptation with intracellular signaling, and homeostasis.

In the present study, regulatory networks for 5 mRNA genes and 3 miRNAs among the exercise-induced DEMs were identified. This study revealed regulatory networks of miRNAs and mRNAs and identified a subset of mRNAs and miRNAs in equine peripheral blood leukocytes affected by exercise. Although this study has been limited by the small number of horses sampled, the lack of methodological uniformity, and underestimation of the influence of circadian rhythms associated with a daytime scheduled exercise on the altered expression patterns of genes, our results provide background information for elucidating molecular homeostasis and adaptation in the exercise-associated physiology of horses. The small number of experimental samples and the lack of equal application of NGS and microarray analysis may have had the greatest effects on the study; thus, the results obtained from the integrated analysis of only two Warmblood horses need to be interpreted cautiously to avoid over-interpretation. A further study involving a greater number of horses with a variety of breeds under controlled conditions, including day to day variations in circadian system, would elucidate the interindividual variations in exercise-related gene expression patterns.