Adequate stress response is a critical factor during athlete horses' training and is central to our capacity to obtain better performances while safeguarding animal welfare.
In order to investigate the molecular mechanisms underlying this process, several studies have been conducted that take advantage of microarray and quantitative real-time PCR (qRT-PCR) technologies to analyse the expression of candidate genes involved in the cellular stress response.
Appropriate application of qRT-PCR, however, requires the use of reference genes whose level of expression is not affected by the test, by general physiological conditions or by inter-individual variability.
The expression of nine potential reference genes was evaluated in lymphocytes of ten endurance horses during strenuous exercise. These genes were tested by qRT-PCR and ranked according to the stability of their expression using three different methods (implemented in geNorm, NormFinder and BestKeeper). Succinate dehydrogenase complex subunit A (SDHA) and hypoxanthine phosphoribosyltransferase (HPRT) always ranked as the two most stably expressed genes. On the other hand, glyceraldehyde-3-phosphate dehydrogenase (GAPDH), transferrin receptor (TFRC) and ribosomal protein L32 (RPL32) were constantly classified as the less reliable controls.
This study underlines the importance of a careful selection of reference genes for qRT-PCR studies of exercise induced stress in horses. Our results, based on different algorithms and analytical procedures, clearly indicate SDHA and HPRT as the most stable reference genes of our pool.
Knowledge of the molecular mechanisms underlying the stress response in athlete horse is a fundamental prerequisite for planning an appropriate training schedule to obtain better performances, preserve animal welfare and avoid overtraining-syndrome [1,2].
It is universally accepted that moderate physical activity may have beneficial effects in terms of general health conditions and could favour the functioning of the immune system. Conversely, strenuous exercise, like exhaustive endurance races, may have detrimental effects on the immune system, determine changes in the cellular composition of peripheral blood and induce the expression of genes that appear to be related to the overtraining-syndrome [3-5]. The list of candidate genes is nevertheless far from being complete, as the athlete's reaction to exercise is a coordinated response of multiple organ systems, and likely involves multiple and complex regulatory changes: induction of heat shock proteins, inflammatory response modulation (pro and anti-inflammatory cytokines) and generation of reactive oxygen and nitrogen species (ROS and RSN) that, besides their damaging potential, play a crucial role in cellular signalling [6,7].
Since exercise has been shown to be an important factor in regulating immune cells and their functions, and considering that stress evokes inflammatory reactions, lymphocytes are considered the best candidate cell type to study physiological changes associated with exhaustive exercise .
Quantitative real-time PCR (qRT-PCR) is the technique of choice when trying to detect modifications in transcription levels in a reliable and reproducible manner. Nevertheless, there are some technical issues that must be taken into account, such as quality and quantification of the starting material, enzyme efficiency, and primer design. Different approaches have been proposed to normalize measurements of expression levels , but this is generally done using an internal control gene, known as a reference gene or as housekeeping gene (HKG), under the assumption that this has a constant level of expression in the chosen tissue, is not affected by the treatment and has no inter-individual variability. In addition, the reference gene and the target gene should have similar ranges of expression to avoid analytical problems.
Widely expressed genes like ACTB, GAPDH or R18S are generally preferred, without preliminary analysis of their expression profiles under the specific study conditions [5,6,9]. Nevertheless, a number of studies report how commonly accepted HKGs do not always constitute reliable controls [10-14], because of unexpected variation in their expression profiles.
More appropriately, multiple HKGs should be evaluated before their employment, and their stability should be measured in the context of the relevant experimental conditions.
The aim of this paper is to identify the best reference genes for qRT-PCR experiments investigating horse lymphocyte gene expression in exercise induced stress. Statistical algorithms implemented in geNorm , BestKeeper , and NormFinder  were used.
To assess which are the most stable genes during strenuous exercise, nine potential HKGs were tested in ten endurance horses with a time course sampling strategy. A qRT-PCR assay, based on SYBR® Green detection, was designed for the transcription profiling of the nine genes (ACTB, B2M, GAPDH, HPRT, R18S, RPL32 SDHA, TFRC and UBB, Table 1). The specificity of the amplifications was confirmed by melting curve analyses (Additional files 1, 2, 3, 4, 5, 6, 7, 8, 9). For each assay, a standard curve was generated by using 4-fold serial dilutions of pooled cDNAs. As shown in Table 2, linear correlation coefficients (R2) varied from 0.998 to 1.000 and PCR efficiencies (E) ranged between 93.9 and 102.6%.
Table 1. Details of the nine genes evaluated.
Table 2. Assay conditions for each of the nine genes evaluated.
Additional file 1. Melting curve ACTB gene. Melting curve analyses image (jpg format) collected using using the MxPro software ver. 3.20 (Stratagene) during calibration experiments of the selected primer pair for the ACTB gene.
Format: JPEG Size: 131KB Download file
Additional file 2. Melting curve B2M gene. Melting curve analyses image (jpg format) collected using using the MxPro software ver. 3.20 (Stratagene) during calibration experiments of the selected primer pair for the B2M gene.
Format: JPEG Size: 131KB Download file
Additional file 3. Melting curve GAPDH gene. Melting curve analyses image (jpg format) collected using using the MxPro software ver. 3.20 (Stratagene) during calibration experiments of the selected primer pair for the GAPDH gene.
Format: JPEG Size: 129KB Download file
Additional file 4. Melting curve HPRT gene. Melting curve analyses image (jpg format) collected using using the MxPro software ver. 3.20 (Stratagene) during calibration experiments of the selected primer pair for the HPRT gene.
Format: JPEG Size: 138KB Download file
Additional file 5. Melting curve R18S gene. Melting curve analyses image (jpg format) collected using using the MxPro software ver. 3.20 (Stratagene) during calibration experiments of the selected primer pair for the R18S gene.
Format: JPEG Size: 140KB Download file
Additional file 6. Melting curve RPL32 gene. Melting curve analyses image (jpg format) collected using using the MxPro software ver. 3.20 (Stratagene) during calibration experiments of the selected primer pair for the RPL32 gene.
Format: JPEG Size: 140KB Download file
Additional file 7. Melting curve SDHA gene. Melting curve analyses image (jpg format) collected using using the MxPro software ver. 3.20 (Stratagene) during calibration experiments of the selected primer pair for the SDHA gene.
Format: JPEG Size: 134KB Download file
Additional file 8. Melting curve TFRC gene. Melting curve analyses image (jpg format) collected using using the MxPro software ver. 3.20 (Stratagene) during calibration experiments of the selected primer pair for the TFRC gene.
Format: JPEG Size: 135KB Download file
Additional file 9. Melting curve UBB gene. Melting curve analyses image (jpg format) collected using using the MxPro software ver. 3.20 (Stratagene) during calibration experiments of the selected primer pair for the UBB gene.
Format: JPEG Size: 136KB Download file
Expression levels of candidate reference genes
Cycle threshold values (Cts) for the nine HKGs tested ranged between 17.9 (ACTB) and 26.6 (TFRC). The gene encoding 18S rRNA is largely over expressed (Ct 9.1) compared to the protein coding genes. Each single control gene appeared to be equally expressed in the tested cDNA samples, and the variations of the Ct values (calculated for each single gene in the ten horse individuals subtracting the Min Ct from the Max Ct values) was always smaller than one (Figure 1).
Figure 1. Average Ct of candidate HKGs. Expression levels of candidate control genes in the ten endurance horses. Values are given as qRT-PCR cycle threshold numbers (Ct values). Circles represent mean Ct values, bars indicate the standard deviation.
Profiles obtained for each horse and HKG were analysed using three different methods, implemented in the software geNorm, NormFinder and BestKeeper.
GeNorm provides a ranking of the tested genes, based on their expression stability, determining the two most stable HKGs for normalization purposes. Selected HKGs were ranked according to the stability measure M (average pair-wise variation of each gene against all others), from the most stable (lowest M value) to the least stable (highest M value): SDHA/HPRT, R18S, B2M, UBB, ACTB, RPL32, TFRC, GAPDH (Table 3). All genes displayed a relatively high stability over the three time course samplings, with M values (M < 0.8) far below the accepted limit of 1.5 . The two most stably expressed genes of our pool (SDHA and HPRT) allow an optimal normalization of qRT-PCR data, and the addition of a third HKG (R18S) would not significantly increase the statistical reliability of this calculation (V2/3 = 0.090, abundantly below the default cut-off value of 0.15 ).
Table 3. Candidate reference genes ranking according to geNorm.
The NormFinder algorithm uses a model-based approach for the estimation of modifications among the HKG expressions, also taking into account variation across sub-groups and avoiding artificial selection of co-regulated genes . The results of the NormFinder analysis are shown in Table 4. This ranking appeared to be slightly different from what obtained using geNorm. GAPDH, TFRC and RPL32 still occupy the lowest positions, while SDHA remains the most stable gene. ACTB gained the second position stepping over HPRT and R18S defined as the least reliable controls.
Table 4. Candidate reference genes ranking according to NormFinder.
BestKeeper measures HKG stability by using a pair-wise correlation analysis of all pairs of candidate genes and calculating the geometric mean of the best candidates [15,19]. A preliminary analysis, based on the inspection of raw Ct values, estimated the variation of all HKGs to be compatible with an overall stability in gene expression (Table 5), with SD values lower than 1. All genes were retained for the calculation of the BestKeeper index, which similarly exhibited a moderate SD variation (0.58). BestKeeper allows a comparative analysis across HKGs, by estimating correlations in the expression levels between all the possible candidates. Highly correlated control genes are combined into an index. Afterwards, the pair-wise correlation between genes and the correlation between each gene and the index are calculated, describing the consistency between the index and each HKG .
Table 5. Statistical output from BestKeeper analysis.
The nine control genes tested in our analysis correlated well one with one another and with the BestKeeper index (Table 5). The best correlation between one HKG and the BestKeeper index was obtained for SDHA (r = 0.983), followed by HPRT, R18S and UBB. The statistically significant correlation shown by SDHA with the BestKeeper index appeared to be consistent with the good performance of this gene as assessed by geNorm and NormFinder. It is also remarkable how GAPDH, RPL32 and TFRC, are again classified as the least reliable HKGs, showing the worst correlations with the determined BestKeeper index (Table 5).
A number of authors have studied gene expression profiles in exercise induced stress using forefront technologies, like gene chips and qRT-PCR. This study is the first solid contribution in assessing which reference genes have to be used to validate and normalize qRT-PCR outcomes.
Several methods have been proposed to allow accurate normalization of gene expression using qRT-PCR [9,19-22] but at present there is no consensus on which algorithm should be used to measure reference gene stability. A comparison of different methods of reference gene selection allows a better identification of the most reliable controls and reduces the risk of artificial selection of co-regulated transcripts .
We compared three different statistical approaches (geNorm, NormFinder and BestKeeper) to evaluate nine potential HKGs, in order to select the best reference gene to be used in studying exercise-induced stress in horses.
The uniformity in gene ranking between the three software packages was generally high: SDHA is the most stable HKG according to all the three methods. HPRT similarly displays a constant significant stability. R18S always ranks third, B2M fourth and UBB fifth, and can be therefore considered plausible HKGs, even if the addition of supplementary reference genes would not significantly enhance the reliability of the normalization according to the geNorm analysis (V value, Table 3).
Regarding ACTB, it is difficult to formulate a final judgement because, as already reported in a previous study , its classification is not consistent between the three software packages (6th in geNorm and BestKeeper, 2nd in NormFinder). Nevertheless, this gene shows an overall reduced variability, as attested by the M value calculated by geNorm and by its good correlation with the BestKeeper index (r = 0.890).
GAPDH, TFRC and RPL32 were classified as the least stable genes and they are not likely to be useful in this given experimental system. Notably, the expression of GAPDH, that has been used as HKG in a previous exercise induced stress study , appears to be the least stable.
In contrast with what reported elsewhere [23-26], R18S appears to be a good potential reference gene. Despite its good performance, the usefulness of this gene as a control is often doubted: some authors  tend to consider it unsuitable for normalization because its transcription is carried out by RNA polymerase I and because of its well known over-expression in comparison with mRNAs (as confirmed even in our experiments, Figure 1). Considering that HKGs that have expression levels comparable to the gene of interest are generally preferred , its usage should be carefully considered if used for the normalization of genes that exhibit low level of expression.
Our results indicate SDHA and HPRT as the most stable reference genes with a very good statistical reliability according to all the three software employed. Moreover, the use of only two genes (SDHA and HPRT) appears to be sufficient for a reliable normalization of the genes of interest; this result is of special interest for future high throughput applications of the technique.
Blood collection, RNA extraction and cDNA synthesis
Ten horses were chosen among participants to national endurance races (90–120 km). Blood samples were taken from the jugular vein and collected at three different time points: before, at the end of the race, and 24 hours after the race. Immediately after collection, peripheral blood mononuclear cells (PBMCs) were isolated by the Ficoll-Hypaque method (GE Healthcare, Pollards Wood, United Kingdom) from 8 ml of blood. Total RNA was extracted from approximately 1 × 107 PBMCs, using the Aurum Total RNA Fatty and Fibrous Tissue kit (Bio-Rad, Herculers CA, USA) according to the manufacturer's instructions. Genomic DNA was eliminated by a DNase treatment supplied with the kit. Extracted RNA was quantified using the Quant-It RNA assay (Invitrogen, Dorset, United Kingdom) in a VersaFluor fluorometer (Bio-Rad) and checked for integrity in a denaturing agarose gel electrophoresis with ethidium bromide staining. Successful removal of DNA contaminants was tested by absence of PCR amplification in the MC1R gene (GenBank accession number X98012, primers from ). 1.5 μg of total RNA were retro-transcribed using random hexamers and Superscript III Reverse Transcriptase (Invitrogen) according to the manufacturer's specifications. A PCR with ACTB primers (Table 2) was performed on each cDNA to check for successful retro-transcription.
Reference genes selection and primer design
Nine widely used reference genes were evaluated: β-actin (ACTB), glyceraldehyde-3P-dehydrogenase (GAPDH), hypoxanthine ribosyltransferase (HPRT1), β-2-microglobin (B2M), succinate dehydrogenase complex subunit A (SDHA), transferrin receptor (TFRC), ubiquitin B (UBB), ribosomal protein L32 (RPL32) and 18S ribosomal rRNA (R18S). In order to minimize the possibility of co-regulation, genes were selected from different functional classes.
Primers were designed based on available sequences using the Primer3 software.
Mfold  was used to check the chosen sequences to avoid designing primers in the region of template secondary structure; amplicon lengths were optimized to 68/138 bp to ensure optimal polymerization efficiency. Specificity of amplification was confirmed by sequencing.
For each primer pair, a preliminary real-time assay was performed to evaluate the amplification of non-specific products or primer dimer artefacts (Additional files 1, 2, 3, 4, 5, 6, 7, 8, 9). Efficiency of RT-PCR (E), slope values, and correlation coefficients (R2) were determined (Table 2) using serial 1:4 dilutions of a template cDNA (pooled from the ten individuals studied, previously diluted 1:5).
PCR products were subsequently resolved on 2% agarose gel to check for size specificity of the amplicon.
Real-time quantitative PCR
Five microliters of cDNA template (previously diluted 1:10) were added to the master mix FastStart SYBR Green Master (Roche Applied Science, Penzberg, Germany) with the ROX fluorochrome internal check. PCR reactions, in a volume of 25 μl were performed on a MX3000P machine (Stratagene, La Jolla CA, USA). PCR conditions were the same for all primer pairs: initial denaturation at 95°C for 10' followed by 40 cycles of denaturation at 95°C for 30", annealing at 58°C for 30" and extension at 72°C for 30". Fluorescence data were collected at the end of the extension step. Following cycling, the melting curve was determined in the range 58°–95°C, with a temperature slope of 0.01°C/sec. Each reaction was run in triplicate with appropriate negative controls.
Baseline and threshold values were automatically determined for all plates and genes using the MxPro software ver. 3.20 (Stratagene). In order to ensure comparability between data obtained from different experimental plates, threshold values for each gene were manually set to the arithmetic mean between the thresholds as automatically determined following each run. Corrected Ct values were transformed to quantities based on the comparative Ct method. Following appropriate formatting, values were imported into geNorm (version 3.4), NormFinder (version 0.953) and BestKeeper (version 1) VBA applets.
KC, MF, SC performed all experiments and data analysis, and drafted the manuscript. GS supervised the study design and the data analysis. MS conceived the project, AVS supervised and coordinated the project and participated in writing the manuscript. All authors read and approved the final manuscript.
Work sponsored by MIUR – PRIN 2006 (coordinator Prof. Maurizio Silvestrelli)
The authors thank Mr Gianluca Alunni for his valuable technical support and Dr. Francesco Nardi for his suggestions in writing this manuscript. We are also grateful with the horse breeders: Mr Giuseppe Sereni, Mr. Carlo Formica, Miss Chiara Rosi for allowing collection of samples and Dr. Marcello Conte for his kind help.
Journal of Endocrinological Investigation 2004, 27:603-612. PubMed Abstract
Zieker D, Fehrenbach E, Dietzsch J, Fliegner J, Waidmann M, Nieselt K, Gebicke-Haerter P, Spanagel R, Simon P, Niess AM, Northoff H: cDNA microarray analysis reveals novel candidate genes expressed in human peripheral blood following exhaustive exercise.
Genes and Immunity 2005, 6:279-284. Publisher Full Text
Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Speleman F: Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes.
Molecular and Cellular Probes 2001, 15:307-311. Publisher Full Text
Bas A, Forsberg G, Hammarstrom S, Hammarstrom ML: Utility of the housekeeping genes 18S rRNA, beta-actin and glyceraldehyde-3-phosphate-dehydrogenase for normalization in real-time quantitative reverse transcriptase-polymerase chain reaction analysis of gene expression in human T lymphocytes.
Biotechniques 2004, 37:112-119. PubMed Abstract
Banda M, Bommineni A, Thomas RA, Luckinbill LS, Tucker JD: Evaluation and validation of housekeeping genes in response to ionizing radiation and chemical exposure for normalizing RNA expression in real-time PCR.
Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP: Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper - Excel-based tool using pair-wise correlations.
Andersen CL, Jensen JL, Orntoft TF: Normalization of real-time quantitative reverse transcription-PCR data: A model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets.
Haller F, Kulle B, Schwager S, Gunawan B, Heydebreck AV, ltmann H, zesi L: Equivalence test in quantitative reverse transcription polymerase chain reaction: Confirmation of reference genes suitable for normalization.
Journal of Biochemical and Biophysical Methods 2000, 46:69-81. Publisher Full Text
Bond BC, Virley DJ, Cairns NJ, Hunter AJ, Moore GBT, Moss SJ, Mudge AW, Walsh FS, Jazin E, Preece P: The quantification of gene expression in an animal model of brain ischaemia using TaqMan real-time RT-PCR.
Biochemical and Biophysical Research Communications 2004, 313:856-862. Publisher Full Text
Rieder S, Taourit S, Mariat D, Langlois B, Guerin G: Mutations in the agouti (ASIP), the extension (MC1R), and the brown (TYRP1) loci and their association to coat color phenotypes in horses (Equus caballus).