Front. Genet.Frontiers in GeneticsFront. Genet.1664-8021Frontiers Media S.A.165162810.3389/fgene.2025.1651628GeneticsOriginal ResearchComparative analysis of blood whole transcriptome profiles in Yili horses pre- and post-5000-meter racingSu et al.10.3389/fgene.2025.1651628SuYi1RenWanlu12MaShikun1MengJun12YaoXinkui12ZengYaqi12LiZexu1LiLuling1WangRan1WangJianwen12*1College of Animal Science, Xinjiang Agricultural University, Xinjiang, China2Xinjiang Key Laboratory of Equine Breeding and Exercise Physiology, Urumqi, China
Edited by:Carrie S. Wilson, Agricultural Research Service (USDA), United States
Reviewed by:Dongyi Bai, Inner Mongolia Agricultural University, China
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
This study employed Yili horses participating in a 5000-meter race as a model to investigate exercise-induced gene expression alterations in peripheral blood using whole transcriptome sequencing. Jugular vein blood samples from the three leading horses were collected pre- and immediately post-race, yielding 2,171 differentially expressed mRNAs (2,080 upregulated, 91 downregulated), 4,375 differentially expressed LncRNAs (4,354 upregulated), and 68 differentially expressed circRNAs (64 upregulated). GO/KEGG analyses demonstrated significant enrichment of differential mRNAs in transmembrane transport function and pivotal signaling pathways (cAMP, MAPK, PI3K-Akt). Differential lncRNAs targeted neuro-signaling pathways (e.g., Neuroactive ligand-receptor interaction, Calcium signaling) and developmental regulators (stem cell pluripotency). Source genes of circRNAs were enriched in axon guidance and immune-related T cell receptor signaling. Molecular functions converged on transporter/receptor activity (mRNA/lncRNA) and nucleic acid/GTP binding (circRNA source genes). The protein-protein interaction analysis identified ten central genes within the heat shock protein family, such as HSP90AA1 and HSPA4. Notably, significant upregulation of HCN4, IGF1, PTHR1, and FGF23 indicated their potential roles in modulating cardiac rhythm, promoting tissue repair, and maintaining calcium-phosphorus homeostasis during exercise adaptation. This study provides comprehensive overview of transcriptomic regulatory mechanisms in the blood of Yili horses, offering a molecular framework for advancing understanding of physiological adaptation to exercise and optimizing equine exercise protocols.
Yili horsesexercise stresswhole transcriptome sequencingdifferentially expressed genessignaling pathwayssection-at-acceptanceLivestock Genomics1 Introduction
As a locally developed breed in China, the Ili horse demonstrates outstanding genetic traits and athletic performance within China’s horse racing industry. However, a significant performance gap remains between this breed and internationally renowned speed horse breeds. Therefore, this study focuses on the Ili horse as the research subject, which is instrumental in accelerating the development process of domestically-bred speed horses within China. The sport horse field is a critical sector within the modern equine industry, and their competitive performance directly impacts economic outcomes. Accurate physiological monitoring is essential to maintaining the competitive health and welfare of sport horses. In recent years, substantial datasets on hematological, physiological, and biochemical parameters across various breeds have been collected (Harari et al., 2024; Giers et al., 2024; Masko et al., 2021), elucidating the dynamic metabolic, immune, and muscular responses to exercise-induced stress (Huangsaksri et al., 2024). Despite this progress, a systematic investigation into exercise-induced global gene expression alterations and the underlying regulatory networks in horses remains lacking.
Exercise, as a potent physiological stimulus, remodels systemic metabolism and adaptive responses through the activation of multifaceted molecular regulatory networks. Evidence indicates that transcriptional regulation triggered by exercise includes essential biological pathways, including energy metabolism, oxidative stress, and inflammatory signaling (Yang et al., 2013; Reitzner et al., 2024). For instance, gene expression analyses in the skeletal muscle of thoroughbreds post-intense exercise revealed marked alterations in genes associated with transcription factor activity, including oxidoreductases and protein-binding elements (Farries et al., 2020). Similarly, transcriptomic profiling of Arabian horses participating in endurance races indicated dynamic shifts in gene networks governing cell migration and tissue repair (Cappelli et al., 2018). Investigations into equine exercise physiology further demonstrate that acute exercise rapidly upregulates genes within the glycolytic pathway, TCA cycle, and oxidative phosphorylation, emphasizing the transcriptional basis of adaptive capacity (Bryan et al., 2017). While epigenetic research has identified correlations between DNA methylation patterns in sport horse blood and key signaling pathways such as PI3K-Akt (Stefaniuk-Szmuki et al., 2018), comprehensive transcriptomic analyses remain limited. Whole-transcriptome approaches enable simultaneous characterization of mRNA expression dynamics alongside the intricate regulation exerted by non-coding RNAs and alternative splicing events, offering a broader framework to elucidate the molecular basis of exercise adaptation.
In this study, peripheral blood samples were collected from Yili horses participating in a 5000-meter race, both pre- and immediately post-race, and subjected to high-throughput whole transcriptome sequencing to comprehensively profile exercise-induced gene expression alterations. Differentially expressed genes (DEGs) and associated regulatory networks were identified, enabling the enrichment of signaling pathways implicated in energy metabolism, inflammatory processes, and tissue repair. These findings contribute to a molecular-level understanding of physiological adaptation in sport horses and support the development of more effective training regimens and health management protocols.
2 Materials and methods2.1 Sample collection
For the 5000-meter speed race of Yili horses conducted in Zhaosu County, Xinjiang’s Ili Kazakh Autonomous Prefecture, all participating horses (24 in total) underwent veterinary health examinations (lameness screening) and breed verification (via passport inspection) prior to the competition. All entrants were required to arrive at the competition terrain by 18:00 on the day preceding the event. At 20:00 that same day, 5 mL jugular venous blood samples were collected from all horses using EDTA anticoagulant vacuum tubes. Within 5 min post-race, venous blood samples were again drawn from the jugular veins of the top three finishers (Group A). For this study, pre-race blood samples corresponding to these top three performers were selected as experimental specimens (Group B). Detailed information on the horses is presented in Table 1. All collected samples were immediately cryopreserved in liquid nitrogen for subsequent analysis.
Equine details.
Horse ID
Age (Y.O.)
Sex
Height at withers (cm)
Trunk length (cm)
Breed
Race results
1
4
Stallion
153
150
Yili Horse
5′23″704
2
4
Stallion
152
147
Yili Horse
5′32″530
3
4
Stallion
153
162
Yili Horse
5′41″339
2.2 Library preparation
Total RNAs were extracted from blood samples, and mRNAs were enriched using mRNA capture magnetic beads. The enriched mRNAs were subsequently fragmented and used for cDNA synthesis. Following cDNA purification with Hieff NGS® DNA Selection Beads, target fragments were screened, and the size-selected fragments were subjected to PCR amplification to construct libraries with an insert size of 350–400 bp High-quality libraries were submitted to Novogen (Beijing, China) for 150 bp paired-end sequencing on the Illumina NovaSeq 6,000 platform (Illumina, CA, USA).
2.3 Bioinformatics analysis
FastQC (fastqc_v0.11.8) was employed to assess the quality metrics of raw Illumina sequencing data. Adapter sequences and low-quality reads were removed using Fastp (fastp 0.23.1), yielding high-quality clean reads. Bismark (version 0.24.0) was subsequently applied to map the clean reads to the Equus Caballus reference genome (EquCab3.0). Reads successfully aligned to the reference genome were designated as target sequences for subsequent standardized and customized analyses.
Filtering parameters included the elimination of reads containing adapter contamination, those with terminal base quality scores below 3, or ambiguous bases (N). A sliding window approach, with a 4-base window and a quality threshold of 15, was used to truncate reads at the point where the average quality within the window fell below the threshold. Reads shorter than 36 nt post-trimming or lacking valid pairing were excluded from further analysis (Ana et al., 2016).
2.4 GO and KEGG enrichment analyses
GO and KEGG enrichment analyses of DEGs and differentially expressed LncRNA target genes were conducted using ClusterProfiler (v3.10.1) and KOBAS (v2.0). GO analysis included three functional categories: biological process (BP), cellular component (CC), and molecular function (MF). Significantly enriched GO terms and KEGG pathways were identified based on a threshold of p adj <0.05.
Based on the intersection of DEGs and known protein interaction pairs from the STRING database, a PPI network was assembled. Homologous protein interaction relationships were integrated into the network, which was subsequently visualized with Cytoscape (v3.8.0). The cytoHubba plug-in was applied to extract the top 10 core genes within the network.
2.6 RT-qPCR
To validate the transcriptome sequencing results, seven mRNAs were randomly selected for real-time quantitative PCR (qPCR) analysis (Hutchinson, 2015). Extract total RNA from whole-blood samples using the TRIZOL reagent (Thermo Fisher Scientific, Catalog # 15596026, USA). Mix 300 μL of whole blood with 900 μL of TRIZOL thoroughly and let it stand for 5 min. Then, add 200 μL of chloroform (Sinopharm Chemical Reagent Co., Ltd., Catalog # 10006818, China), centrifuge at 12,000 rpm at 4 °C for 10 min. Collect the supernatant, add 500 μL of isopropyl alcohol (Sinopharm Chemical Reagent Co., Ltd., Catalog # 80109218, China) to precipitate the RNA, wash it with 75% ethanol, and dissolve it in nuclease-free water (Servicebio, Catalog #G4700, China). Measure the RNA concentration and purity using the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA) with an A260/A280 ratio of 1.8–2.0. Then, 2 μg of total RNA was reversely transcribed into cDNA using HiScript® Q RT SuperMix kit (TOYOBO Life Science, Cat# FSQ-101, Japan) in a 20 μL reaction mixture under the condition of 37 °C for 15 min and 98 °C for 5 min in a Veriti thermal cycler (Thermo Fisher Scientific, Waltham, MA, USA). The amplification was carried out using ChamQ SYBR qPCR Master Mix (TOYOBO Life Science, Cat# QPK-201, Japan) on a TL-988 real-time PCR system (Xi’an Tianlong Technology Co., Ltd., China). Primers with specific specificity were designed using Primer 5 software (Table 2), and the 20 μL reaction mixture consisted of 10 μL of pre-mix solution, 0.4 μM forward and reverse primers, and 1 μL of cDNA template. The amplification procedure involved two stages: Stage 1: 95 °C for 30 s for primer denaturation; Stage 2 (40 cycles): 95 °C for 15 s for denaturation, followed by 60 °C for 30 s for annealing/extension. Stage 3 (Melting Curve): From 65 °C to 95 °C, collect one fluorescence signal per 0.5 °C increment. All samples have three technical repeats. Gene relative expression levels were calculated using the 2−ΔΔCT method, with GAPDH serving as the reference gene for normalization (Supplementary Table S1).
RT-qPCR primers.
Gene name
Primer sequences (3'→5')
Product size (bp)
GAPDH
F: CATCAAATGGGGCGATGCTG
158
R: GGTTCACGCCCATCACAAAC
CACYBP
F:AGTCCCCACTGAGAATGTGC
127
R: GCCTTCCACAGAGATGGGTT
DNAJB1
F: CGTCGGACGAGGAGATCAAG
160
R: CCGAAGCGGTCGAAAATGTC
HSP90AA1
F: CACAGGTGAGACCAAGGACC
105
R: CAGTACTCGTCGATCGGCTC
LOC100054696
F: CGATGTTTTGGGGGTCAAACC
81
R: GGGTGGTACTTCAAGGCCAG
LOC100072672
F: CAACCCGGAGGACAAGTACC
103
R: ACCGCAAGGCATTTCATCT
PTGES3
F: ACATGGGTGGTGATGAGGATG
117
R: CCGGTGATGGTAACATTCCTT
SLC5A3
F: GGAAGCGCTGCAATGAACAA
136
R: TTGATGAAGCCTGGCCTGTT
3 Results3.1 Sequencing quality analysis
Blood transcriptome data were obtained from six samples, comprising early-stage (Group B) and late-stage (Group A) time points of the race. Quality control metrics, summarized in Table 3, indicated that Group B and Group A yielded 238,441,444 and 246,501,422 clean reads, respectively. Q20 and Q30 values exceeded 98.00% and 94.00%, while GC content ranged from 46.35% to 49.92%,Mapped reads ≥90% for all samples. Both groups demonstrated sequencing pass rates above 90%, confirming the robustness and reliability of the datasets for further analyses.
Sequencing quality and read counts.
Sample
Raw reads
Clean reads
Error rate
Q20 (%)
Q30 (%)
GC content/(%)
Mapped reads
B1
78,769,258
76,747,064
0.01
98.54
96.12
49.92
71,729,558 (93.46%)
B2
86,281,582
84,453,076
0.01
98.1
94.88
49.61
76,454,194 (90.53%)
B3
78,833,206
77,241,304
0.01
98.67
96.42
49.35
72,705,768 (94.13%)
A1
83,030,990
81,319,316
0.01
98.61
96.32
47.76
75,915,707 (93.36%)
A2
85,764,908
82,977,310
0.01
98.43
95.95
46.69
74,818,193 (90.17%)
A3
84,442,934
82,204,796
0.01
98.66
96.35
49.41
76,141,631 (92.62%)
3.2 Gene expression profiling
The FPKM distribution, visualized via box plots (Figure 1A), reflected consistent global gene expression across samples. The x-axis denoted sample identifiers, while the y-axis represented log10(FPKM). Minimal dispersion and comparable expression ranges among samples supported the data’s reproducibility. Pearson correlation analysis (Figure 1B) further indicated clear intergroup distinction, with low intra-group variability, validating the appropriateness and reliability of sample selection for comparative transcriptomic analysis.
mRNA expression level pre- and post-race. (A) mRNA expression level; (B) Pearson correlation coefficient diagram. Note: (A) illustrate the distribution of gene expression levels across six biological replicates (B1-B3, A1-A3), with box plots showing the transformed expression values after logarithmic transformation by base 10 (FPKM +1). The vertical axis represents the transformed expression values, while the boxes represent interquartile range (IQR). The lines extend beyond 1.5 times IQR for outliers. (B) Display a heatmap depicting sample-wise correlations in expression, with color depth indicating the strength of Pearson correlation coefficient (R2) (ranging from light purple 0.7 to dark blue 1.0). Diagonal lines indicate self-correlation within the same samples (R2 = 1.0).
Panel A shows box plots of FPKM distribution across six samples (B1, B2, B3, A1, A2, A3) with log-transformed values. Panel B is a heat map of Pearson correlation coefficients between these samples, with values ranging from 0.657 to 1.0, indicating varying levels of correlation.3.3 Analysis of differential mRNAs, LncRNAs, and circRNAs
Differential expression was determined using |log2(Fold Change)| > 0 and Padj <0.05 as selection thresholds. Comparative analysis between Group A and Group B identified 2,171 differentially expressed mRNAs, including HCN4, IGF1, PTHR1, and FGF23, with 2,080 upregulated and 91 downregulated (Figure 2A). A total of 4,375 differentially expressed LncRNAs were detected, comprising 4,354 upregulated and 21 downregulated transcripts (Figure 3A). Additionally, 68 circRNAs displayed differential expression, with 4 upregulated and 64 downregulated (Figure 4A). The cluster analysis results demonstrated that the differential mRNAs, LncRNAs, and circRNAs in blood samples before and after the competition, as shown in Figures 2B, 3B, 4B, were highly reproducible, revealing significant differences between the groups.
Volcano plot (A) and heatmap (B) illustrating mRNA expression differences between group B and group A.
Volcano plot and heatmap showing gene expression analysis. Panel A is a volcano plot with log2 fold change against negative log10 adjusted p-values, highlighting UP (red), DOWN (green), and NO change (blue) in genes. Panel B is a heatmap displaying hierarchical clustering of gene expression levels before and after treatment, with red indicating upregulation and green indicating downregulation. Groups are labeled 'After' in blue and 'Before' in pink.
Volcano plot (A) and heatmap (B) illustrating LncRNA expression differences between group B and group A.
Panel A shows a volcano plot with points indicating gene expression changes: red for upregulated (4354), green for downregulated (21), and blue for unchanged (8592). Panel B is a heatmap illustrating hierarchical clustering of gene expression, with red indicating higher and green lower expression levels. The samples are categorized into two groups, After and Before, color-coded in pink and blue.
Volcano plot (A) and heatmap (B) illustrating circRNA expression differences between group B and group A.
Panel A displays a volcano plot with data points in blue, green, and red, showing changes in expression based on fold change and p-value. Blue represents no change, green is downregulated, and red is upregulated. Panel B is a heat map with hierarchical clustering, comparing expression levels before and after treatment. Colors represent expression levels, with red indicating high and green indicating low expression.3.4 GO and KEGG analyses of differential mRNAs, LncRNAs, and circRNAs
GO enrichment analysis and database annotation (Figure 5A) identified significant (P < 0.05) enrichment of differential mRNAs in 17 BP and 26 MF categories. Among the GO terms, transmembrane transport, transporter activity, and transmembrane transporter activity were prominently represented. KEGG pathway analysis (Figure 6A) further revealed that DEGs were enriched in signaling cascades including cAMP, MAPK, and PI3K-Akt pathways.
Three vertical bar charts labeled A, B, and C depict data categorized into Biological Processes (BP) in red, Cellular Components (CC) in green, and Molecular Functions (MF) in blue. The y-axis shows negative log-transformed p-values, while the x-axis lists descriptions. Each chart highlights varying data patterns across categories, with specific values indicated above the bars. Legends clarify the color-coded categories.
Three dot plots labeled A, B, and C show gene enrichment analysis results. Each plot displays pathways on the y-axis and gene ratio on the x-axis. Dot size indicates count, and color represents adjusted p-values, with a spectrum from blue to red, where blue represents higher p-values and red represents lower p-values. Plot A highlights pathways like cholinergic synapse and MAPK signaling. Plot B shows pathways such as axon guidance and Cushing syndrome. Plot C highlights pathways like focal adhesion and proteoglycans in cancer.
GO enrichment analysis was subsequently performed on differential LncRNAs (Figure 5B), revealing significant enrichment (P < 0.05) across 11 BP, 1 CC, and 27 MF categories. The enriched BPs were the G-protein-coupled receptor signaling pathway, multicellular organismal process, and cell adhesion. The sole enriched CC was the extracellular region. Within MF, enriched terms included transporter activity, molecular function regulator, signal transducer activity, molecular transducer activity, and signaling receptor activity. KEGG pathway analysis (Figure 6B) further demonstrated significant enrichment of differential LncRNAs in pathways including Neuroactive ligand-receptor interaction, cAMP signaling, Calcium signaling, regulation of stem cell pluripotency, and cGMP-PKG signaling.
GO-based functional annotation of the differentially expressed circRNAs (Figure 5C) identified enrichment across 68 corresponding source genes. In the BF category, annotations were primarily associated with nucleic acid metabolic processes and regulation of biosynthetic activity. CC terms predominantly involved cellular and organelle-related localizations. Within the MF category, functional associations were concentrated in GTP binding, DNA binding, and related molecular interactions. KEGG pathway enrichment analysis further revealed significant representation in signaling pathways such as Axon guidance and T cell receptor signaling (Figure 6C). GO and KEGG enrichment can be found in the supplementary material (Supplementary Table S3; Supplementary Table S4).
3.5 PPI network analysis
A PPI network was established based on DEGs identified in the blood of Yili horses pre- and post-race (Figure 7), enabling the identification of the top 10 core genes: HSP90AA1, HSPA4, HSP90AB1, STIP1, KDM6B, HSPA4L, HSPA8, PTGES3, HSPD1, and TCP1. Each exhibited markedly elevated expression.
PPI network of DEGs and identification of core genes. Note:The outer blue nodes indicate differentially expressed genes, while the inner yellow nodes (n = 10) denote hub genes, with lines representing significant co-expression relationships (Pearson |r| > 0.9).
Network diagram showing interconnected nodes. Central yellow nodes labeled with genes like HSP90AA1 and STIP1 are surrounded by blue nodes representing other genes. Lines indicate interactions among nodes, suggesting a complex biological network.3.6 RT-qPCR verification
The fluorescence quantification outcomes of the seven mRNAs, as depicted in Figure 8, exhibited concordance with the RNA-seq results, confirming the reliability of the transcriptomic data analysis. This consistency validated the RNA-seq dataset as a sound basis for subsequent investigations (Supplementary Table S1).
Validation of RNA-seq data using RT-qPcR for Yili horses, Data are expressed as mean t standard error of mean (SEM). (A) Indicate Log2Fcobtained by RNA-seq and RT-qPCR of DEGs. (B) Represent the relative expressions of DEGs, in RT-qPcR.
Bar graphs compare gene expression levels. Chart A shows log2 fold change between RT-qPCR and RNA-seq. Chart B displays relative expression in groups B and A. Genes listed on the x-axis include C1H100RF128, CACYBP, DNAJB1, HSP90AA1, LOC100054696, LOC100072672, PTGES3, and SLC5A3.4 Discussion
During racing, equine physiological adaptation requires coordinated modulation of multiple signaling pathways to maintain systemic homeostasis. In this study, whole transcriptome sequencing of pre- and post-race blood samples from Yili horses identified 2,171 differentially expressed mRNAs. These transcripts represent potential targets for regulating exercise-related responses in horses. Notably, HCN4, IGF1, PTHR1, and FGF23 emerged as shared differential mRNAs, with all four exhibiting significant post-race upregulation relative to the remaining 2,077 transcripts. This expression pattern suggests a potential link between elevated mRNAs expression and enhanced exercise capacity. As a “central pacemaker” within the sinoatrial node (SAN), HCN4 responds to increased cAMP induced by physical activity via its C-terminal cyclic nucleotide binding domain (CNBD) (Porro et al., 2020), enhancing the pacemaker current (If) to elevate heart rate and maintain cardiac output stability during 5000-meter races. This mechanism is further substantiated by the accumulation of cAMP signal pathways and their direct relevance to the protection against exercise-induced arrhythmias (Lin et al., 2024). HCN4 plays a central role in regulating cardiac rhythmicity. Studies in rabbit models have demonstrated that increased HCN4 expression supports sinoatrial node recovery and mitigates arrhythmogenic events, highlighting its contribution to heart rate modulation and myocardial energy balance (Tuerxun et al., 2024). Yukihiro et al. (Saito et al., 2022) demonstrated that targeted integration of the HCN4 gene into the AAVS1 locus via TALEN-mediated editing induced HCN4 channel protein expression and elevated IF current, confirming the role of HCN4 overexpression in enhancing the pacing capacity of cardiomyocytes. Post-race analysis revealed a marked upregulation of HCN4 expression, suggesting that exercise may trigger activation of the HCN4 gene, thereby increasing channel density on the cardiomyocyte membrane and modulating cardiac pacing to support enhanced physical performance. Nevertheless, an excessive elevation in heart rate may compromise cardiac output or induce arrhythmogenic events, ultimately impairing exercise capacity and heightening cardiovascular risk. IGF1 contributes to the synthesis of proteoglycans, supports cellular viability, and attenuates cytokine-driven catabolic processes (Li et al., 2024; Jin et al., 2025; Menetrey et al., 2000). IGF1 elevation facilitates dystrophin synthesis by activating the PI3K-Akt signaling pathway, thereby reducing muscular membrane damage (Song et al., 2021). The corresponding protein plays a regulatory role in maintaining cellular metabolic homeostasis and promoting tissue repair mechanisms (Christopoulos et al., 2015). In this study, post-race upregulation of the DEG IGF1 in Yili horses suggests a potential involvement in tissue repair. Vigorous exercise during racing is known to induce muscle cell damage (Flick et al., 2021; Wilson et al., 2023); in response, elevated IGF1 expression may contribute to cellular repair and preservation of cytoskeletal integrity (Mir et al., 2017). Activation of the cAMP signaling pathway is triggered when PTHR1 binds to PTH or PTHrP (Singh et al., 2005; Gesty-Palmer et al., 2006), promoting renal calcium reabsorption while concurrently suppressing phosphate reabsorption (Liu et al., 2012). This cascade also enhances bone resorption and intestinal calcium uptake, contributing to the maintenance of systemic calcium levels (Fu et al., 2021). During high-intensity exercise, increased calcium demand for muscle contraction and energy metabolism may stimulate PTH secretion, which, through PTH1R, enhances renal calcium reabsorption and skeletal calcium release to sustain calcium homeostasis (Goltzman and Hendy, 2015). FGF23, predominantly secreted by osteocytes and osteoblasts, regulates phosphate balance by modulating renal phosphate clearance and vitamin D metabolism (Lee et al., 2014; Liu et al., 2015; Hondares et al., 2014). By binding to the Klotho protein and FGFR receptor, FGF23 initiates downstream signaling in target tissues including the kidneys, parathyroid glands, and bones (Wang et al., 2024; Shacham et al., 2023; Javier et al., 2025). According to Li et al. (Li et al., 2016), exercise-induced FGF23 modulated skeletal muscle oxidative stress by limiting ROS overproduction and enhancing mitochondrial function to improve exercise performance. In the present analysis, significant upregulation of FGF23 in post-race blood samples implicates its potential role in optimizing exercise performance through metabolic regulation in horses.
GO enrichment analysis revealed significant enrichment of differentially expressed mRNAs in BP and MF categories in the blood of Yili horses before and after racing. KEGG pathway analysis indicated marked enrichment in cAMP, MAPK, and PI3K-Akt signaling pathways. The cAMP pathway mediates intracellular signal transduction by modulating cAMP synthesis and degradation in response to extracellular stimuli, thereby influencing metabolic activity, gene expression, cellular differentiation, and survival (Glebov-McCloud et al., 2024; Zhang et al., 2024). Intense physical exertion may trigger an adrenaline surge that activates the cAMP-PKA axis, triggering glycogen phosphorylation, thereby enhancing muscle glycogen breakdown rate. Upon ligand binding to G protein-coupled receptors (GPCRs), adenylate cyclase (AC) is activated, catalyzing the conversion of ATP into cAMP and promoting energy production (Akizuki et al., 2021; Heckman et al., 2018). The MAPK pathway contributes to cellular proliferation and differentiation and is potentially associated with inflammatory signaling (Jia et al., 2025). The Mitogen-activated Protein Kinase (MAPK) pathway modulates cytokine generation (such as Interleukin-6) in cells during intense exercise in horses, mitigating muscle damage. The Phosphatidylinositol-3-Kinase/AKT cascade regulates cellular survival, growth, and metabolic function (Goyal et al., 2022). In the Iliyun horse race, IGF1 elevation is facilitated by stimulating proteoglycan synthesis to accelerate muscle repair, and by enhancing GLUT4 translocation to promote glycogen replenishment for metabolic adaptation. These pathways likely support energy homeostasis during intensive exercise and modulate post-race inflammatory responses in horses.
The PPI network constructed from DEGs in Yili horse blood pre- and post-race revealed that the top 10 core genes were all associated with heat shock proteins, exhibiting marked upregulation. HSP90AA1 contributes to cellular stress responses and the maintenance of homeostasis (Goyal et al., 2022), while HSPA4, regulated by heat shock factor, is upregulated under stressors such as elevated temperature to maintain intracellular homeostasis (Li et al., 2019). Additionally, KDM6B is implicated in the modulation of inflammatory processes (Yang et al., 2024). These core genes may serve as candidate regulators influencing Yili horse race performance, although the underlying molecular mechanisms remain to be elucidated.
A total of 4,375 differentially expressed LncRNAs were detected in the peripheral blood of Yili horses pre- and post-race. The high expression levels observed across both time points suggested that these LncRNAs were not merely transcriptional byproducts of mRNA but may exert specific regulatory roles. GO enrichment of the predicted target genes revealed significant associations with biological regulation and metabolic processes. KEGG pathway analysis indicated enrichment in pathways such as Neuroactive ligand-receptor interaction, cAMP signaling, Calcium signaling, regulation of stem cell pluripotency, and cGMP-PKG signaling. Neuroactive ligand-receptor interactions have been implicated in mediating rapid cellular responses via ionotropic receptors (e.g., glutamate receptors) and GPCRs (Neele et al., 2018), In the 5000-meter race, LncRNAs might enhance neuronal excitability, boost muscle contraction precision and reaction speed, thereby supporting high-intensity endurance performance. In the present analysis, 124 differential LncRNAs were enriched in this pathway, suggesting their potential role in modulating neuromuscular responsiveness of Yili horses during racing. Additionally, intracellular Ca2+ serves as a critical signaling molecule regulating muscle contraction and neurotransmitter secretion within the calcium signaling pathway (Chan et al., 2025). LncRNAs may selectively target genes such as RyR1, modulating the release of calcium ions from the sarcoplasmic reticulum, thereby balancing explosive power with endurance. The enrichment of 88 differential LncRNAs in this pathway implies their involvement in exercise-related physiological processes, potentially by modulating neurotransmitter and endocrine hormone release during race conditions.
To investigate the association between circRNA expression and exercise in Yili horses, GO and KEGG enrichment analyses were conducted on the corresponding source genes. GO annotation revealed that the differentially expressed circRNA source genes were primarily linked to BP, with six genes implicated in nucleic acid circulation pathways. These results imply a potential involvement of circRNAs in modulating immune responses during exercise in Yili horses. KEGG analysis indicated that the source genes were predominantly enriched in metabolic and signal transduction pathways, including Axon guidance and T cell receptor signaling. Axon guidance relies on a combinatorial interaction of receptors and ligands to mediate cellular adhesion, attraction, or repulsion responses (de Almeida et al., 2016), During the competition, DEcirRNAs target genes such as ROBO1/SLIT2 regulate axons of motor neurons, which may potentially affect muscle contractions. While T cell antigen receptors regulate the activation and proliferation of T lymphocytes (Seiradake et al., 2016). In this study, CircRNAs by modulating genes such as CD3E, inhibit excessive inflammatory responses post-exercise and promote tissue repair, thereby elucidating the molecular basis for the quick recovery seen in the Iliki horse after intense exercise, consistent with prior research (Masko et al., 2021), suggesting a functional contribution of the corresponding circRNAs to neural signaling and immune regulatory mechanisms.
Reverse transcription quantitative polymerase chain reaction (RT-qPCR) successfully confirmed seven mRNA trends; however, tissue-specific mechanisms such as heart muscle HCN4 were not directly validated. In this study, we did not verify the identified DELncRNAs and DEcirRNAs by qRT-PCR, primarily due to the consistency between the pathways selected in our KEGG and GO enrichment analyses and those of the mRNAs. This led us to overlook the need for qRT-PCR analysis on LncRNAs and cirRNAs. We must admit that this indeed does affect our research on DELncRNAs and DEcirRNAs; in subsequent studies, we will need to delve deeper into these aspects. Despite the small sample size (n = 3), this study’s statistical analysis ensured reliability through stringent thresholds. In the data results displayed, the gene expression patterns remained clear, demonstrating the feasibility of the findings, and upon RT-qPCR validation, the reliability of the sequencing data was confirmed. This study has successfully established reliable data results in a small sample size, proving the feasibility and stability of the research technique.
5 Conclusion
In this study, transcriptome sequencing was employed to profile the whole transcriptome of Yili horse blood collected before and after a 5000-meter race. DEGs were predominantly enriched in the cAMP, MAPK, and PI3K-Akt signaling pathways, suggesting potential impact of racing on these items and signaling pathways. HCN4, IGF-1, PTHR1, and FGF23 were annotated as exercise-related candidate genes in Yili horses. The findings offer a foundation for further exploration of candidate genes associated with equine exercise performance.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/,PRJNA1250149.
Ethics statement
The animal studies were approved by The study was conducted in accordance with the Declaration of Helsinki, and approved by the Animal Welfare and Ethics Committee of Xinjiang Agricultural University (approval code: 2024003 and date of approval: 22 April 2024). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent was obtained from the owners for the participation of their animals in this study.
Author contributions
YS: Conceptualization, Data curation, Formal Analysis, Methodology, Software, Writing – original draft, Writing – review and editing. WR: Conceptualization, Formal Analysis, Methodology, Writing – review and editing. SM: Methodology, Conceptualization, Formal Analysis, Writing – review and editing. JuM: Investigation, Resources, Supervision, Validation, Writing – review and editing. XY: Investigation, Project administration, Resources, Supervision, Writing – review and editing. YZ: Formal Analysis, Resources, Software, Writing – review and editing. ZL: Data curation, Software, Writing – review and editing. LL: Formal Analysis, Visualization, Writing – review and editing. RW: Visualization, Writing – review and editing. JiW: Conceptualization, Funding acquisition, Project administration, Writing – original draft, Writing – review and editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This research was funded by the National Natural Science Foundation of China (grant number 32302735), the Major Science and Technology Special Project of Xinjiang Uygur Autonomous Region (grant number 2022A02013-1), and the Central Guidance for Local Science and Technology Development Fund (grant number ZYYD2025JD02).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2025.1651628/full#supplementary-material
ReferencesAkizukiK.OnoA.XueH. C.KameshitaI.IshidaA.SueyoshiN. (2021). Biochemical characterization of four splice variants of mouse Ca2+/calmodulin-dependent protein kinase Iδ. J. Biochem.169, 445–458. 10.1093/jb/mvaa11733417706AnaC.PedroM.SoniaT.Gomez-CabreroD.CerveraA.McPhersonA. (2016). A survey of best practices for RNA-Seq data analysis. Genome Biol.17 (1), 13. 10.1186/s13059-016-0881-826813401BouwmanW.VerhaeghW.van DoornA.RaymakersR.van der PollT.van de StolpeA. (2024). Quantitative characterization of immune cells by measuring cellular signal transduction pathway activity. Sci. Rep-UK14, 24487. 10.1038/s41598-024-75666-w39424625BryanK.McGivneyB. A.FarriesG.McGettiganP. A.McGivneyC. L.GoughK. F. (2017). Equine skeletal muscle adaptations to exercise and training: evidence of differential regulation of autophagosomal and mitochondrial components. BMC Genomics18, 595. 10.1186/s12864-017-4007-928793853CappelliK.CapomaccioS.ViglinoA.SilvestrelliM.BeccatiF.MoscatiL. (2018). Circulating miRNAs as putative biomarkers of exercise adaptation in endurance horses. Front. Physiol.9, 429. 10.3389/fphys.2018.0042929740341ChanR. J.WalkerA.VardyJ.ChanA.OppegaardK.ConleyY. P. (2025). Perturbations in the neuroactive ligand-receptor interaction and renin angiotensin system pathways are associated with cancer-related cognitive impairment. Support Care Cancer33, 254. 10.1007/s00520-025-09317-940047999ChristopoulosP. F.MsaouelP.KoutsilierisM. (2015). The role of the insulin-like growth factor-1 system in breast cancer. Mol. Cancer14, 43. 10.1186/s12943-015-0291-725743390de AlmeidaR. M. C.ClendenonS. G.RichardsW. G.BoedigheimerM.DamoreM.RossettiS. (2016). Transcriptome analysis reveals manifold mechanisms of cyst development in ADPKD. Hum. Genomics10, 37. 10.1186/s40246-016-0095-x27871310FarriesG.BryanK.McGivneyC. L.McGettiganP. A.GoughK. F.BrowneJ. A. (2020). Expression quantitative trait loci in equine skeletal muscle reveals heritable variation in metabolism and the training responsive transcriptome. Front. Genet.10, 1215. 10.3389/fgene.2019.0121531850069FlickM.VintherA. M. L.JacobsenS.BergL. C.GimenoM.VerwilghenD. (2021). Effect of exercise on serum neutrophil gelatinase-associated lipocalin concentration in racehorses. Vet. Clin. Path50, 551–554. 10.1111/vcp.1302734779025FuX. K.ZhouB.YanQ. N.TaoC.QinL.WuX. H. (2021). Kindlin-2 regulates skeletal homeostasis by modulating PTH1R in mice. Signal Transduct. Tar5, 297. 10.1038/s41392-020-00328-y33361757Gesty-PalmerD.ChenM.ReiterE.AhnS.NelsonC. D.WangS. (2006). Distinct beta-arrestin- and G protein-dependent pathways for parathyroid hormone receptor-stimulated ERK1/2 activation. J. Biol. Chem.281, 10856–10864. 10.1074/jbc.M51338020016492667GiersJ.BartelA.KirschK.MuellerS. F.HorstmannS.GehlenH. (2024). Blood-based assessment of oxidative stress, inflammation, endocrine and metabolic adaptations in eventing horses accounting for plasma volume shift after exercise. Vet. Med. Sci.10, e1409. 10.1002/vms3.140938516822Glebov-McCloudA. G. P.SaideW. S.GaineM. E.StrackS. (2024). Protein kinase A in neurological disorders. J. Neurodev. Disord.16, 9. 10.1186/s11689-024-09525-038481146GoltzmanD.HendyG. N. (2015). The calcium-sensing receptor in bone-mechanistic and therapeutic insights. Nat. Rev. Endocrinol.11, 298–307. 10.1038/nrendo.2015.3025752283GoyalA.AgrawalA.VermaA.DubeyN. (2022). The PI3K-AKT pathway: a plausible therapeutic target in parkinson's disease. Exp. Mol. Pathol.129, 104846. 10.1016/j.yexmp.2022.10484636436571HarariS.DeretzS.Saint PriestB. D.RichardE.RicardA. (2024). Comparison of blood parameters in two genetically different groups of horses for functional longevity in show jumping. Front. Genet.15, 1455790. 10.3389/fgene.2024.145579039534078HeckmanP. R. A.BloklandA.BollenE. P. P.PrickaertsJ. (2018). Phosphodiesterase inhibition and modulation of corticostriatal and hippocampal circuits: clinical overview and translational considerations. Neurosci. Biobehav R.87, 233–254. 10.1016/j.neubiorev.2018.02.00729454746HondaresE.Gallego-EscuredoJ. M.FlachsP.FrontiniA.CereijoR.GodayA. (2014). Fibroblast growth factor-21 is expressed in neonatal and pheochromocytoma-induced adult human brown adipose tissue. Metabolism63, 312–317. 10.1016/j.metabol.2013.11.01424369918HuangsaksriO.SanigavateeK.PoochipakornC.WonghanchaoT.YalongM.ThongchamK. (2024). Physiological stress responses in horses participating in novice endurance rides. Heliyon10, e31874. 10.1016/j.heliyon.2024.e3187438845925HutchinsonJ. A. (2015). Quantification of mRNA expression by RT-qPCR. Transplantation99 (10), 2009–2011. 10.1097/TP.000000000000094826426919JavierA. J. S.KennedyF. M.YiX.Wehling-HenricksM.TidballJ. G.WhiteK. E. (2025). Klotho is cardioprotective in the mdx mouse model of Duchenne muscular dystrophy. Am. J. Pathol.29, 923–940. 10.1016/j.ajpath.2024.12.01739889824JiaW.DanH.LinaH.QinzhenZ.ShaozhuangL.JiapengHu (2025). Protective effect of micheliolide against inflammation and oxidative stress in asthma through the MAPK/NF-κB signalling pathway. World Allergy Organ. J.18 (8), 101091. 10.1016/J.WAOJOU.2025.10109140718157JinR.PeiH. L.YueF.ZhangX. D.ZhangZ. C.XuY. (2025). Network pharmacology combined with metabolomics reveals the mechanism of yangxuerongjin pill against type 2 diabetic peripheral neuropathy in rats. Drug Des. Dev. Ther.19, 325–347. 10.2147/DDDT.S47314639834645LeeJ.HongS. W.ParkS. E.RheeE. J.ParkC. Y.OhK. W. (2014). Exendin-4 regulates lipid metabolism and fibroblast growth factor 21 in hepatic steatosis. Metabolism63, 1041–1048. 10.1016/j.metabol.2014.04.01124933399LiD. J.FuH.ZhaoT.NiM.ShenF. M. (2016). Exercise-stimulated FGF23 promotes exercise performance via controlling the excess reactive oxygen species production and enhancing mitochondrial function in skeletal muscle. Metabolism65, 747–756. 10.1016/j.metabol.2016.02.00927085781LiP. P.WangJ. J.ZouY.SunZ. L.ZhangM. H.GengZ. M. (2019). Interaction of Hsp90AA1 with phospholipids stabilizes membranes under stress conditions. BBA-Biomembranes1861, 457–465. 10.1016/j.bbamem.2018.11.00930517848LiJ. T.GaoH.WangP.SunC.WeiZ. L.YiX. C. (2024). Plumbagin induces G2/M arrest and apoptosis and ferroptosis via ROS/p38 MAPK pathway in human osteosarcoma cells. Alex Eng. J.103, 222–236. 10.1016/j.aej.2024.06.015LinZ. Y.LinB. W.HangC. W.LuR. H.XiongH.LiuJ. Y. (2024). A new paradigm for generating high-quality cardiac pacemaker cells from mouse pluripotent stem cells. Signal Transduct. Tar9, 230. 10.1038/s41392-024-01942-w39237509LiuF.FuP.FanW. X.GouR.HuangY. Q.QiuH. Y. (2012). Involvement of parathyroid hormone-related protein in vascular calcification of chronic haemodialysis patients. Nephrology17, 552–560. 10.1111/j.1440-1797.2012.01601.x22448974LiuJ.XuY.HuY. J.WangG. (2015). The role of fibroblast growth factor 21 in the pathogenesis of non-alcoholic fatty liver disease and implications for therapy. Metabolism64, 380–390. 10.1016/j.metabol.2014.11.00925516477MaskoM.DominoM.JasinskiT.Witkowska-PilaszewiczO. (2021). The physical activity-dependent hematological and biochemical changes in school horses in comparison to blood profiles in endurance and race horses. Animals11, 1128. 10.3390/ani1104112833920044MenetreyJ.KasemkijwattanaC.DayC. S.BoschP.VogtM.FuF. H. (2000). Growth factors improve muscle healing in vivo. J. Bone Jt. Surg. Br.82, 131–137. 10.1302/0301-620x.82b1.895410697329MirS.CaiW. K.CarlsonS. W.SaatmanK. E.AndresD. A. (2017). IGF-1 mediated neurogenesis involves a novel RIT1/Akt/Sox2 Cascade. Sci. Rep-UK7, 3283. 10.1038/s41598-017-03641-928607354NeeleA. E.GijbelsM. J. J.van der VeldenS.HoeksemaM. A.BoshuizenM. C. S.PrangeK. H. M. (2018). Myeloid Kdm6b deficiency results in advanced atherosclerosis. Atherosclerosis275, 156–165. 10.1016/j.atherosclerosis.2018.05.05229908485PorroA.ThielG.MoroniA.SaponaroA. (2020). Cyclic AMP regulation and its command in the pacemaker channel HCN4. Front. Physiol.11, 771. 10.3389/fphys.2020.0077132733276ReitznerS. M.EmanuelssonE. B.ArifM.KaczkowskiB.KwonA. T. J.MardinogluA. (2024). Molecular profiling of high-level athlete skeletal muscle after acute endurance or resistance exercise- A systems biology approach. Mol. Metab.79, 101857. 10.1016/j.molmet.2023.10185738141850SaitoY.NakamuraK.YoshidaM.SugiyamaH.AkagiS.MiyoshiT. (2022). Enhancement of pacing function by HCN4 overexpression in human pluripotent stem cell-derived cardiomyocytes. Stem Cell Res. Ther.13, 141. 10.1186/s13287-022-02818-y35365232SeiradakeE.JonesE. Y.KleinR. (2016). Structural perspectives on axon guidance. Annu. Rev. Cell Dev. Bi32, 577–608. 10.1146/annurev-cellbio-111315-12500827576119ShachamE. C.MamanN.LazarevaT.MasalhaR.MahagnaL.SelaG. (2023). Normocalcemic primary hyperparathyroidism is an early stage of primary hyperparathyroidism according to fibroblast growth factor 23 level. Front. Endocrinol.14, 1152464. 10.3389/fendo.2023.115246437065752SinghA. T.GilchristA.Voyno-YasenetskayaT.Radeff-HuangJ. M.SternP. H. (2005). G alpha12/G alpha13 subunits of heterotrimeric G proteins mediate parathyroid hormone activation of phospholipase D in UMR-106 osteoblastic cells. Endocrinology146, 2171–2175. 10.1210/en.2004-128315705779SongX. T.ZhangJ. N.ZhaoD. W.ZhaiY. F.LuQ.QiM. Y. (2021). Molecular cloning, expression, and functional features of IGF1 splice variants in sheep. Endocr. Connect.10, 980–994. 10.1530/EC-21-018134319906Stefaniuk-SzmukierM.Ropka-MolikK.PiórkowskaK.ZukowskiK.Bugno-PoniewierskaM. (2018). Transcriptomic hallmarks of bone remodelling revealed by RNA-seq profiling in blood of Arabian horses during racing training regime. Gene676, 256–262. 10.1016/j.gene.2018.07.04030021131TuerxunZ.HeY. X.NiuY. X.BaoZ.LiuX. M.YangY. C. (2024). Analysis of differentially expressed murine miRNAs in acute myocardial infarction and target genes related to heart rate. Cell Biochem. Biophys.83, 963–975. 10.1007/s12013-024-01528-x39325365WangS. S.XuQ.ZhangY.CaoH. X.WangN.HuY. F. (2024). The FGF23-Klotho axis promotes microinflammation in chronic kidney disease. Cytokine184, 156781. 10.1016/j.cyto.2024.15678139454251WilsonJ.De DonatoM.AppelbaumB.GarciaC. T.PetersS. (2023). Differential expression of innate and adaptive immune genes during acute physical exercise in American quarter horses. Animals13, 308. 10.3390/ani1302030836670847YangQ. J.YangG. J.WanL. L.LiB.LuJ.YuQ. (2013). Serum metabolic profiles reveal the effect of formoterol on cachexia in tumor-bearing mice. Mol. Biosyst.9, 3015–3025. 10.1039/c3mb70134d24056883YangJ. H.WuX. X.YouJ. H. (2024). Unveiling the potential of HSPA4: a comprehensive pan-cancer analysis of HSPA4 in diagnosis, prognosis, and immunotherapy. Aging16, 2517–2541. 10.18632/aging.20549638305786ZhangZ. Y.XuZ.WangS.JiaZ. Z.ZhouZ.WangC. (2024). Optimized new Shengmai Powder modulation of cAMP/Rap1A signaling pathway attenuates myocardial fibrosis in heart failure. Chin. Med-UK19, 30. 10.1186/s13020-024-00902-438402401