The Pusan National University-Institutional Animal Care and Use Committee approved the study design (Approval Number: PNU-2015-0864). A skeletal muscle tissue biopsy was performed on the leg of a neonatal Thoroughbred to culture horse muscle cells (HMCs). HMCs were cultured according to the previous study [11]. Briefly, primary HMCs were maintained and sub-passaged in Medium 199 (Gibco, Grand Island, NY, USA) supplemented with 10% foetal bovine serum (Invitrogen, Carlsbad, CA, USA) and 1% antibiotic-antimycotic (Gibco, USA). HMCs were incubated in a humidified atmosphere with 5% CO₂ at 37°C. Routine fluid renewals were made three times a week. At 70% to 80% confluence, cells were gently washed twice with phosphate-buffered saline and harvested using 0.05% trypsin-ethylenediaminetetraacetic acid (Welgene, Daegu, Korea) to isolate total RNA. We chose 42°C for heat stress according to the physiological condition of horse skeletal muscle tissue after exercise [12]. A value of 2% of O₂ was determined to induce the hypoxic condition according to the previous work using a hypoxic chamber [13].

Total RNA was extracted from HMCs by RNA isolation reagent (TRIzol, Invitrogen, Karlsruhe, Germany). Next, 200 μL of chloroform was added to remove cells from the organic solvent, then the mixture was shaken for 10 s, maintained at 4°C for 5 min, and then centrifuged at 4°C for 15 min. The supernatant was removed and added to a new test tube, mixed with an equal amount of isopropanol, and maintained at 4°C for 15 min to collect RNA pellets. Isopropanol was removed from the solution by centrifuging at 4°C for 15 min, then was sterilized with 75% ethanol and dissolved in RNase-free water. The purity of the extracted RNA was confirmed by measuring absorbance at 230 nm and 260 nm using a spectrophotometer (Nanodrop Technologies Inc., Wilmington, DE, USA), and only RNA with a purity (optical density value of 230 nm/260 nm) greater than 1.8 was selected and stored at −70°C for further experiments. To synthesize cDNA, 1 μg of RNA, 1 μL of oligo-dT (Invitrogen, Waltham, MA, USA) and RNase-free water were added. The RNA was denatured at 80°C for 3 min, and cDNA was synthesized using 4 μL of 5× reverse transcription (RT) buffer, 5 μL of 2 mM dNTPs, 0.5 μL of RNase inhibitor and 1 μL of M-MLV (Moloney-murine leukaemia virus) RT (Promega, Madison, WI, USA).

NCBI (http://www.ncbi.nlm.nih.gov) and the Ensembl Genome Browser (www.ensembl.org) were used to retrieve gene sequence information. The primers for amplification of the genes (Table 1) were designed using PRIMER3 software (http://bioinfo.ut.ee/primer3-0.4.0/). Real-time quantitative polymerase chain reaction (qPCR) was performed using a thermal cycler (C1000 Thermal Cycler, Bio Rad, Hercules, CA, USA) to measure the relevant expression of target genes in 25 μL of reaction solution, constructed as follows: 2 μL diluted cDNA (50 ng/μL) added to 14 μL SYBR green master mix (Bio Rad, USA) and 1 μL each of 5 pmol/μL diluted forward and reverse primers. The conditions used for the real-time qPCR were as follows: initial denaturation at 94°C for 10 min followed by 40 cycles of denaturation at 94°C for 10 s, annealing at 60°C for 10 s, and extension at 72°C for 30 s. All measurements were carried out in triplicate, and the 2^−ΔΔCt method was used to determine relative gene expression. The relative expression of target genes was normalized to glyceraldehyde-3-phospate dehydrogenase.

Cell viability was assayed by measuring blue formazan that was metabolized from MTT by mitochondrial dehydrogenase. The HMCs were re-suspended in the medium one day before H₂O₂ and cortisol treatment, at a density of 2×10⁵ cells per well in 24-well culture plates. Liquid medium was replaced with fresh medium containing dimethyl sulfoxide (DMSO) for control. The HMCs were incubated with various concentrations of H₂O₂ (Junsei, Tokyo, Japan) and cortisol (Merck, Darmstadt, Germany). MTT (5 mg/mL) was added to each well and incubated for 4 h at 37°C. The formazan product formed was dissolved by adding 200 μL DMSO to each well, and the absorbance was measured at 570 nm on a microplate reader (Ultra Multifunctional Microplate Reader, Tecan US Inc., Durham, NC, USA). All measurements were performed in triplicate and repeated at least three times.

The HMCs were plated on 60-mm culture dishes at a concentration of 5×10⁵ cells/plate and incubated at 37°C under 5% CO₂ for 4 h under each stressor. Then, cell images were acquired by using an inverted microscope.

The DEG data from the previous research was used for bioinformatics analysis [10]. As all genes were listed in Ensembl gene IDs, the horse Ensemble gene IDs were converted to official gene symbols by cross-matching with human Ensembl gene IDs and the official gene symbols. The official gene symbols of human homologues of equine genes were used for functional clustering and enrichment analyzes using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) tool. The functional categories based on co-occurrence with a set of upregulated genes in DEGs were analyzed using the Kyoto encyclopedia of genes and genomes (KEGG) pathway tool in the DAVID. The statistical significance of pathway over-representation was calculated based on the DAVID software EASE score (p value based on a modified Fisher’s exact t-test) and corrected for multiple testing using the Benjamini-Hochberg stepdown correction.

Gene network analysis using GeneMANIA (https://genemania.org) was conducted with the top six genes from the DEG data. Database information including co-expression, co-localization, shared protein domains, and predicted functional relationships between genes were collected from GeneMANIA.

All data for qRT-PCR were expressed as mean±standard deviation from three independent experiments. GraphPad Prism software (GraphPad Software, La Jolla, CA, USA) was used to evaluate the data. The statistical significance (* p< 0.05, ** p<0.01, or *** p<0.001) was assessed by two sample student’s t-test.

We used the upregulated DEGs data derived from our previous study [10]. We cut off the list based on a false discovery rate (FDR) of <0.01. Then, a set of genes was sorted depending on Log2fold change (Log₂FC) value (Log₂FC>1). Finally, we obtained 1,383 DEGs that included 1,186 known and 197 unknown upregulated genes. To analyze pathways and gene networks that relate to the list of upregulated genes, we used the DAVID bioinformatics tool with the known genes. The 1,186 genes upregulated after horse exercise in skeletal muscle tissue were involved in the multiple biological pathways derived from the database of the KEGG pathway (Figure 1). The pathways were related to inflammatory responses such as cytokine–cytokine receptor interaction, MAPK signalling pathway, tumor necrosis factor signalling pathway, chemokine signalling pathway, NF-kappa B signalling pathway, toll-like receptor signalling pathway, and Nucleotide-Binding Oligomerization Domain-like receptor pathway. Among the 1,186 known upregulated genes, 17 genes (1.4%) were significantly upregulated after exercise (Log₂FC>7). To deepen the gene networks, we selected the top six upregulated genes: interleukin 6 (IL6), IL8, chemokine (C-X-C motif) ligand 6 (CXCL6), A disintegrin and metalloproteinase with thrombospondin motifs 4 (ADAMTS4), heat shock 70 kDa protein 6 (HSPA6), and selectin E (SELE) (Table 2). We then performed GeneMANIA analysis to determine the relationships among the six genes. The network consisted of 26 genes, including the six selected genes and 20 additional genes identified by GeneMANIA (Figure 2A). Co-expression (82.16%), co-localization (9.67%), shared protein domains (4.76%), and predicted functional relationships (3.41%) were confirmed by literature searching. The pathway result consisted of 26 genes, including six selected genes and 20 additional genes (Figure 2B). The common pathways of DAVID and GeneMANIA were related to an inflammatory response such as cytokine–cytokine receptor interaction and the chemokine signalling pathway. These results indicate that the upregulated genes after horse exercise in skeletal muscle tissue strongly relate to the stress response.

After horse exercise, the horse skeletal muscle tissue is subjected to various stresses. We selected the representative exercise stressors including heat, oxidative, hypoxic, and hormone stress. Then, we established an in vitro system to validate the effect of each stressor on the transcription of the top six upregulated genes. We selected H₂O₂ and cortisol to apply oxidative and hormone stress, respectively. The effect of oxidative stress with H₂O₂ at 200 μM to 1 mM on cell viability of HMCs was assessed by the MTT assay. H₂O₂ at 800 μM reduced cell viability by approximately 75% compared to control (Figure 3A). When hormone stress with cortisol at 20 μg/mL to 80 μg/mL was applied, cortisol at 40 μg/mL reduced cell viability by approximately 75% compared to control (Figure 3B). We applied 600 μM H₂O₂ and 20 μg/mL cortisol to choose a concentration at which the cell viability was unaffected. Next, 42°C and 2% O₂ were applied for heat and hypoxic stress, respectively.

Then, we used several marker genes to confirm the induction of each stressor including heat shock 70 kDa protein 1 (Hsp72), heme oxygenase-1 (HO-1), cyclooxygenase-2 (COX-2), and hypoxia-inducible factor 1-alpha (HIF1A) for heat, oxidative, cortisol, and hypoxic stress, respectively [14–17]. Transcription of Hsp72, HO-1, and HIF1A increased significantly and COX-2 decreased (Figure 4A). Moreover, cell morphologies of each stress treatment were similar to control (Figure 4B). These results show that heat, oxidative, hypoxic and cortisol stress were induced under non-lethal conditions.

Next, we validated the effect of the exercise-related stressors on the transcriptional induction of IL6, IL8, CXCL6, ADAMTS4, HSPA6, and SELE. IL8 (* p<0.05) and HSP6A (*** p<0.001) increased while IL6 (* p<0.05), ADAMTS4 (*** p<0.001), and SELE (*** p<0.001) decreased after heat stress (Figure 5A). Interestingly, HSP6A (*** p<0.001) increased significantly, almost 400 times more than control. In oxidative stress with 600 μM of H₂O₂, IL6 (*** p<0.001), IL8 (** p<0.01), and CXCL6 (** p<0.01) significantly increased while ADAMTS4 (*** p< 0.001) and SELE (*** p<0.001) significantly decreased (Figure 5B). Hypoxic condition with 2% O₂ increased CXCL6 (* p< 0.05). However, the other genes were downregulated except HSPA6, which showed a slight increase (Figure 5C). Cortisol stress with 20 μg/mL showed decreased expression of all six genes except IL8 (* p<0.05) (Figure 5D). Although all six genes were upregulated after horse exercise in silico, these data clearly indicate that heat, oxidative, hypoxic, and hormone stress contribute differently to the transcriptional regulation of the top six upregulated genes. Also, oxidative stress led to expression similar to that seen in silico.