Experiment overview

Eight four-year old unconditioned Thoroughbred horses (castrated males) were exercised to maximum heart-rate or fatigue in a standardized incremental-step exercise test [25–27] on a high-speed equine treadmill. Skeletal muscle biopsy samples were collected at three time points: at rest pre-exercise (T₀), immediately post-exercise (T₁) and four hours post-exercise (T₂). In a direct comparison microarray experiment, equine cDNA microarrays were hybridised with samples from T₀Vs T₁ and from T₀Vs T₂ for each animal.

Exercise parameters

Following warm-up, the exercise test comprised an average of six (range 5 - 7) incremental steps achieving a mean maximum velocity of 12.4 ± 0.2 m/s and a mean distance of 4,362.9 ± 102.7 m for an average duration of 8.77 ± 0.5 min. Mean maximal heart rate was 218 ± 9 beats per minute. Mean peak post-exercise (T₁) lactate concentrations were 13.3 ± 1.2 mmol/l and were significantly increased compared to pre-exercise values (P < 0.0001).

Microarray annotation and gene ontology

Of the 9,333 ESTs on the microarray 8,519 aligned to a single location on the equine genome (EquCab 2.0), 372 aligned to more than one location and the remaining 442 failed to align to any location with high confidence. Fewer than 50% (4,631) of the ESTs matched an Ensembl gene, the majority (4,166) of which had human orthologs. The human orthologs were used to create input files for gene ontology functional analyses using the DAVID software package [28, 29].

The functional representation of ESTs on the microarray relative to all genes in the Equus caballus Ensembl database that had human orthologs (66%) was assessed using 15 broad GO categories (developmental process, multicellular organismal process, biological regulation, metabolic process, cellular process, macromolecular complex, organelle part, organelle, cell part, cell, transporter activity, transcription regulator activity, molecular transducer activity, catalytic activity, binding). A similar distribution pattern among GO categories was observed for ESTs on the microarray when compared to all Ensembl genes (Additional file 1).

Immediate response to exercise

Delayed response to exercise

Validation of a panel of genes by real time qRT-PCR

Nine genes that were found to be differentially expressed in the microarray experiment were selected for validation by real time qRT-PCR. These genes were chosen based on their involvement in muscle contraction or the response to hypoxia. Two probes (Genbank IDs: CX594334 and CX598227) that showed differential expression, but were not found within an annotated gene were also included for validation. CX594334 lies ~ 2 kb upstream of PFKFB3 and was upregulated immediately post-exercise. CX598227 lies ~ 1 kb downstream of calmodulin 1 (Calm1) and was upregulated four hours post exercise. The average gene expression of seven of the nine probes studied reached significance (P < 0.05) and six [(basic helix-loop-helix family, member e40 (BHLHE40), calmodulin 3(CALM3), HSPA1A, FOS, CX594334 and CALM1] were concordant with the microarray data (Table 6, Figure 1). A point of major concern in microarray studies is the presence of false positives within a gene list. Although the use of qRT-PCR is essential to validate the overall dataset it is not feasible to interpret the experimental findings by evaluating each gene individually. As genes function in co-operation within complex networks we report principally the expression patterns of functionally related groups of genes.

There were however, some interesting findings among the validated genes. For instance, the expression of FOS showed a high inter-sample variance in gene expression change estimates (+18.2-fold to +506.0-fold increase). The mean expression of the FOS gene mRNA transcript increased +198-fold. The biological significance of the high gene expression variance for this gene is not clear at present but warrants further investigation.

Transcripts for the HSPA1A gene and the probe CX594334 (which may represent the PFKFB3 gene) were significantly increased at T₁ (2.3-fold and 2.7-fold fold respectively). While both genes have quite different physiological roles [42, 92] both contain hypoxic response elements (HRE) and have been shown to be transcriptionally activated by the HIF-1α protein under hypoxic conditions [93].

During the recovery period, four hours post exercise, HSPA1A mRNA levels remained elevated (+5.9-fold) while CX594334 transcript levels returned to baseline (Figure 1). The BHLHE40 gene, which increased in expression +2.3-fold is a transcription factor involved in the hypoxic response, contains a HRE and is inducible in hypoxic conditions through interaction with HIF-1α [94, 95]. CALM3 and CX598227 which lies ~ 1 kb downstream of CALM1 showed directionally different changes in gene expression. CALM3 was downregulated -0.82-fold while CX598227 was upregulated +2.4-fold. This is of particular interest as little is known regarding the differential regulation of the individual genes within the Calmodulin gene family. Calmodulin is a calcium binding protein which acts as a calcium sensor [96] and plays an important role in mediating many cellular processes including muscle contraction [97, 98].

Subjects

Eight four-year old unconditioned Thoroughbred horses (castrated males), raised at the same farm for the previous 12 months and destined for National Hunt racing with the same trainer comprised the study cohort. The horses had a mean height of 165.25 cm ((± 1.44) and a mean pre-exercise weight of 565.75 kg (± 13.71). All horses participated in a standardized incremental-step exercise test [25–27] on a high-speed equine treadmill (Sato, Sato AB, Knivsta, Sweden). Before the exercise test, all horses were judged to be clinically healthy based on a veterinary examination that included a lameness assessment, resting upper airway endoscopy and basic bloodwork (complete blood count and serum biochemistry). Prior to entering the study, all of the horses had been raised together and had been kept in a grass field and fed 1.8 kg of 14% Race horses cubes (Gain horse feeds, Clonroche, Co. Wexford, Ireland) three times a day. During the study week the horses were housed in a stable and provided ab libitum access to water and fed grass hay and 2 kg of 10% Cool-n-Cooked Horse and Pony Mix (Connolly's Red Mills, Bagnelstown, Co. Carlow, Ireland) twice daily. Horses were fed approximately 3 hours and 55 minutes (235 ± 0.11 minutes) prior to the exercise test. All exercise tests were performed between 1000 - 1130 am.

Standardised exercise test

The treadmill was housed in an insulated room with temperature and humidity monitors. Prior to the exercise test, all horses were acclimatized to stand quietly and to comfortably transition gaits on the treadmill. The treadmill was set to a 6° incline for all of the exercise tests. The warm-up period consisted of 2 minutes at 2 m/s, followed by 2 minutes at 4 m/s and then 2 minutes at 6 m/s. This was then followed by an increase in treadmill velocity to 9 m/s for 60 seconds, and then a 1 m/s increase in treadmill velocity every 60 seconds until the animal was no longer able to maintain its position on the treadmill at that speed or until the heart rate reached a plateau (HR_max). Following completion of the test, the horses were quickly brought back to a walk, taken off the treadmill and washed down with cold water.

Muscle biopsy sampling

Percutaneous needle muscle biopsies [99] were obtained from the dorsal compartment of the gluteus medius muscle according to Dingboom and colleagues [100] using a 6 mm diameter, modified Bergstrom biopsy needle (Jørgen KRUUSE, Veterinary Supplies). Biopsies were taken approximately 15 cm caudodorsal (one-third of the distance) to the tuber coxae on an imaginary line drawn from the tuber coxae to the head of the tail. The biopsies were obtained at a depth of 80 mm. Each biopsy site was shaved, scrubbed with an antiseptic and desensitized by a local anaesthetic. The biopsy samples were washed with sterile PBS (BD Biosciences, San Jose, CA) and preserved in RNAlater (Ambion, UK) for 24 hours at 4°C and then stored at -20°C.

Muscle biopsy samples were collected at three time points: at rest pre-exercise (T₀), immediately post-exercise (T₁) and four hours post-exercise (T₂). Pre-exercise biopsies were collected within 93 minutes (range 68 - 93 mins) before the commencement of exercise and between 155 and 170 minutes post feeding. The immediately post-exercise (T₁) biopsies were collected within seven minutes 30 seconds (range 5 mins 45 sec - 7 mins 30 sec) following cessation of exercise and four hour post-exercise (T₂) biopsies were collected within 262 minutes (range 242 - 262 mins) following cessation of exercise.

RNA isolation and purification

Approximately 100 mg of each muscle biopsy sample was removed from RNAlater and homogenized in 3 ml TRIzol using a Kinematica Polytron Homogeniser PT 1200 C Drive unit, 230 V (AGB, Dublin, Ireland) and the aqueous and organic phases were separated using 200 μl of chloroform. Total RNA was precipitated using isopropyl alcohol (0.6 times the volume of the aqueous phase). The remaining pellet was washed once in 75% ethanol, and redissolved in 35 μl of nuclease-free water (Promega UK Ltd, Southampton, UK). Each sample was purified using the RNeasy ^® Mini kit (Qiagen Ltd, Crawley, UK) and DNase treated with RNase free DNase (Qiagen Ltd, Crawley, UK). To elute the total RNA, 35 μl of nuclease-free water were applied to the silica-gel membrane of the column to elute the total RNA, which was stored at -80°C. RNA was quantified using a NanoDrop^® ND1000 spectrophotometer V 3.5.2 (NanoDrop Technologies, Wilmington, DE) and RNA quality was subsequently assessed using the 18S/28S ratio and RNA integrity number (RIN) on an Agilent Bioanalyser with the RNA 6000 Nano LabChip kit (Agilent Technologies Ireland Ltd, Dublin, Ireland) according to manufacturer's instructions. The RNA isolated from these samples had an average RNA integrity number (RIN) of 8.43 ± 0.08 (range 8.0 - 9.0).

Microarray description and annotation

Microarray slides were printed with clones selected from a cDNA library generated using mRNA purified from the articular cartilage of a 15-month old male Thoroughbred horse [101]. Probe sets on the microarray slides were prepared and printed as previously described [102, 103].

The cDNA sequences for all annotated genes on the Equus caballus Version 2.0. (EquCab2.0) genome sequence http://www.broad.mit.edu/mammals/horse/ were downloaded from Ensembl http://www.ensembl.org, release 50. The expressed sequence tag (EST) sequences of all the probes on the array were masked to remove repeats using RepeatMasker [104] and blast searched against the cDNAs. BLAST hits were filtered to retain only hits with e values ≤ 10^-10, ≥ 50 bp long, above 95% match-target identity, and where the best hit e value was ≥ 10¹⁰ better than the next best. The EST sequences were cross-matched to horse Entrez gene IDs and to human Ensembl and Entrez gene IDs via accessions. EST matches to multiple horse or human genes were excluded. Because fewer than 50% of ESTs matched an Ensembl gene predicted gene annotations were assigned to unannotated probes of interest based on the gene located closest to the probe and homology to mammalian genes. Following a BLAST of RefSeq and protein databases search hits with an e value of < 1^-10 and hits on non-mammalian species were eliminated. If the predicted gene annotations based on location and homology did not match, the probe was not assigned a predicted annotation.

Microarray hybridisation and experimental design

Total RNA was amplified using a MessageAmpTM amplified RNA (aRNA) linear amplification kit (Ambion). 2 μg of aRNA was reverse transcribed and directly labelled with Fluor647 or Fluor555 using the SuperScript™ Plus Direct cDNA Labeling System (Bio-sciences, Dublin Ireland) according to the manufacturer's instructions. Labelled samples were combined and co-hybridised on equine cDNA microarrays using SlideHyb Glass Array Hybridisation Buffer #3 (Applied Biosystems, Cambridgeshire, UK). Microarray hybridisations were performed on an automated HS400 hybridisation station (Tecan Group Ltd. Seestrasse 103 CH-8708 Männedorf, Switzerland) with the following protocol - wash: 75°C, runs 1, wash 10 s, soak 20 s; probe injection: 85°C; denaturation: 95°C, 2 min; hybridisation: 42°C, agitation frequency medium, 4 h; hybridisation: 35°C, agitation frequency medium, 4 h; hybridisation: 30°C, agitation frequency medium, 4 h; wash: 37°C, runs 2, wash 10 s, soak 20 s; wash: 25°C, runs 2, wash 15 s, soak 30 s; wash: 25°C, runs 2, wash 20 s, soak 40 s; slide drying: 25°C 2 min.

The experimental design was a direct comparison for each animal between pre- and both post-exercise time points. Each slide was hybridised with samples from T₀Vs T₁ and from T₀Vs T₂ for each animal. Technical replicates in the form of a dye swap were performed for each investigation.

Microarray scanning and data acquisition

Hybridised and dried slides were scanned using a GenePix 4000B scanner (Molecular Devices, Berkshire, UK) and image acquisition, first-pass data analysis and filtering were carried out using the GenePix 6.0 microarray image analysis package (Molecular Devices, Berkshire, UK). As a first step of feature extraction spots that were flagged as 'poor' by the GenePix software (due to signal foreground or background contamination, shape irregularity or poor spot quality) were assigned a weight of zero and were excluded from differential expression analyses. Images of the slides were visually examined and any obvious irregularities were also flagged, assigned a weight of zero and excluded from differential expression analyses. All spots with a mean signal of less than background plus two standard deviations were flagged and were also excluded from differential expression analyses.

Microarray data analyses

All statistical analyses on the gene expression data were performed using the R language, version 2.5.1 [105] and the packages statmod and LIMMA from the Bioconductor project [106]. Robust multichip average (RMA) [107] and print tip lowess normalization [108, 109] were performed on the data before differential expression analyses were performed using the lmFit function in LIMMA. Fluor647:Fluor555 log₂ ratios were calculated for each spot on the microarray and duplicate spots were averaged. The function duplicateCorrelation [110] was used to estimate the correlation between technical replicates (dye swaps) by fitting a mixed linear model by REML individually for each probe. The function also returned a consensus correlation, which is a robust average of the individual correlations. This was used to fit a linear model to the expression data for each probe taking into account the inter-technical replicate correlation between each microarray hybridisation.

Differentially expressed targets were determined using a Bayes moderated t-test [111]. Multiple testing was addressed by controlling the false discovery rate (FDR) using the correction of Benjamini and Hochberg [112]. A probe was flagged as differentially expressed if the corrected P value was < 0.05.

Functional clustering according to gene ontology annotations

A list of Entrez IDs of human homologs of probes on the microarray was obtained in a similar manner as for microarray annotation. Using the Entrez IDs of human homologues of equine genes it was possible to use the Database for Annotation, Visualization and Integrated Discovery (DAVID) [28, 29] for functional clustering and overrepresentation analyses. DAVID was used to investigate the representation of broad gene ontology (GO) categories (Level 1) on the equine cDNA microarray relative to the whole genome.

The DAVID system was also used to cluster differentially expressed genes according to their function. For T₀vs T₂ experiments, a probe was called differentially expressed if its corrected P value was < 0.05 [112]. The enrichment of categories was assessed and compared with the proportion observed in the total population of genes on the microarray, using the Expression Analysis Systematic Explorer (EASE) tool within DAVID [113]. A different approach was used when functionally clustering differentially expressed genes from the T₀vs T₁ experiments. Although we expected there would only be small number of genes differentially expressed at T₁, there remains the possibility that some more modest but still genuine changes in gene expression may not be detected. Therefore the FatiScan [56, 114] gene enrichment test was used to analyse the transcriptional profile immediately after exercise. FatiScan is part of the Babelomics Suite Genes and tests for the asymmetrical distribution of biological labels in an ordered list of genes. Genes were ranked by differential expression and FatiScan was used to detect functional blocks (GO and KEGG pathways) that were significantly up-regulated and down-regulated immediately after exercise.

Quantitative real time RT-PCR

Selected cDNA samples from seven of the eight animals were quantified by real time qRT-PCR. One of the animals was omitted due to a shortage of RNA in the pre-exercise sample. 1 μg of total RNA from each sample was reverse transcribed into cDNA with oligo-dT primers using a SuperScript™ III first strand synthesis SuperMix kit according to the manufacturer's instructions (Invitrogen Ltd, Paisley, UK). The converted cDNA was diluted to 2.5 ng/μl working stocks and stored at -20°C for subsequent analyses. Oligonucleotide primers for real time qRT-PCR were designed using Primer3 version 3.0 http://www.primer3.sourceforge.net and commercially synthesized (MWG Biotech, Germany), details of these primers are available in additional file 5. Each reaction was carried out in a total volume of 20 μl with 2 μl of cDNA (2.5 ng/μl), 10 μl SYBR^® Green PCR Master Mix (Applied Biosystems, Cambridgeshire, UK) and 8 μl primer/H₂O. Optimal primer concentrations were determined by titrating 50, 300 and 900 nM final concentrations and disassociation curves were examined for the presence of a single product. qRT-PCR was performed using a 7500 Fast Real-Time PCR machine (Applied Biosystems, Cambridgeshire, UK).

A panel of four putative reference or 'housekeeping' genes were selected for a reference gene study. This panel comprised two frequently used reference genes (HPRT1, hypoxanthine phosphoribosyltransferase 1 gene; and GAPDH, glyceraldehyde-3-phosphate dehydrogenase gene) and two genes (NSUN6, NOL1/NOP2/Sun domain family, member 6 gene; and PIGO, phosphatidylinositol glycan anchor biosynthesis, class O gene) that were selected based on minimal variation across the time points observed in the microarray results. The panel of genes was evaluated using geNorm version 3.4 for Microsoft Excel [115]. Briefly, the gene expression stability measure 'M' for each control gene was calculated as the pairwise variation for that gene with all other tested reference genes across the exercise time-course (T₀, T₁, T₂). The candidate reference genes were ranked in order of decreasing 'M' values or increasing mRNA expression stability [85]. Based on the geNorm analyses, the NSUN6 gene was the optimal reference gene and pre-exercise (T₀) values were used to normalise the data. The 2^-ΔΔCT method (where CT is cycle threshold) was used to determine mean fold changes in gene expression [116]. The Student's t-test was used to identify significant differences in mRNA abundance between time-points.