Exercise induced stress in horses: Selection of the most stable reference genes for quantitative RT-PCR normalization

Background 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. Results 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. Conclusion 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.


Background
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][4][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 [5].
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 [8], but this is generally done using an internal control gene, known as a ref-erence 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][11][12][13][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.
A number of statistical methods have been proposed to evaluate stability of gene expression and select the best HKGs in a given experimental setting [9,[15][16][17][18].
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 [9], BestKeeper [19], and NormFinder [20] were used.

Results
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 (R 2 ) varied from 0.998 to 1.000 and PCR efficiencies (E) ranged between 93.9 and 102.6%.

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).

Data analysis
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 pairwise 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 [9]. 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 (V 2/3 = 0.090, abundantly below the default cut-off value of 0.15 [9]).
The NormFinder algorithm uses a model-based approach for the estimation of modifications among the HKG expressions, also taking into account variation across subgroups and avoiding artificial selection of co-regulated genes [20]. 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.
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). Best-Keeper 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 [15].
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   worst correlations with the determined BestKeeper index (Table 5).

Discussion
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][20][21][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 coregulated transcripts [16].
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 [15], its classification is not consistent between the three software packages (6 th in geNorm and BestKeeper, 2 nd 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 [5], appears to be the least stable.
In contrast with what reported elsewhere [23][24][25][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 [27] 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 [15], its usage should be carefully considered if used for the normalization of genes that exhibit low level of expression.

Conclusion
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
Primers were designed based on available sequences using the Primer3 software.
Mfold [29] 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.
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.