Selection of reference genes for normalization of quantitative real-time PCR in organ culture of the rat and rabbit intervertebral disc

Background The accuracy of quantitative real-time RT-PCR (qRT-PCR) is often influenced by experimental artifacts, resulting in erroneous expression profiles of target genes. The practice of employing normalization using a reference gene significantly improves reliability and its applicability to molecular biology. However, selection of an ideal reference gene(s) is of critical importance to discern meaningful results. The aim of this study was to evaluate the stability of seven potential reference genes (Actb, GAPDH, 18S rRNA, CycA, Hprt1, Ywhaz, and Pgk1) and identify most stable gene(s) for application in tissue culture research using the rat and rabbit intervertebral disc (IVD). Findings In vitro, four genes (Hprt1, CycA, GAPDH, and 18S rRNA) in rat IVD tissue and five genes (CycA, Hprt1, Actb, Pgk1, and Ywhaz) in rabbit IVD tissue were determined as most stable for up to 14 days in culture. Pair-wise variation analysis indicated that combination of Hprt1 and CycA in rat and the combination of Hprt1, CycA, and Actb in rabbit may most stable reference gene candidates for IVD tissue culture. Conclusions Our results indicate that Hprt1 and CycA are the most stable reference gene candidates for rat and rabbit IVD culture studies. In rabbit IVD, Actb could be an additional gene employed in conjunction with Hprt1 and CycA. Selection of optimal reference gene candidate(s) should be a pertinent exercise before employment of PCR outcome measures for biomedical research.


Background
Quantitative real-time RT-PCR (qRT-PCR) is a powerful tool for detection and quantification of gene expression owing to its high sensitivity, specificity and reproducibility [1]. However, the relevance and magnitude of absolute measures obtained from qRT-PCR are subject to inherent sample variations that usually lead to statistical uncertainty. Such outcomes may also lead to inferences that may be biologically obscure [2]. To this end, an appropriate normalization strategy is employed for reliable data interpretation and the most common method is the use of an internal reference or housekeeping gene [3]. A reference gene is weakly regulated in experimental conditions of interest along with comparable expression characteristics to target genes. Thus, selection process of an ideal reference gene is critical for applicability of PCR in research. However, commonly used reference genes for tissue and cell-based PCR normalization such as glyceraldehydes-3-phosphate dehydrogenase (GAPDH), β-actin (Actb) and 18S ribosomal RNA (18S rRNA) are frequently applied without appropriate validation for their gene expression stability [4][5][6].
Gene expression analyses are used as one of the major contemporary research tools in understanding the pathology of intervertebral disc (IVD) degeneration. Clinically termed as "Degenerative Disc Disease (DDD)", the condition is believed to be a significant source of low back pain [7,8]. The clinical significance of understanding the onset and progress of DDD is well documented and there is increasing need to establish relevant experimental models to study this disease. Numerous studies on molecular level changes in disc biology have been reported by normalization using GAPDH [9][10][11][12] and Actb [13][14][15] without validation for their stability. Although the assumption of certain genes being constitutively expressed may be valid in certain cases, this assumption cannot be taken for granted under rapidly changing conditions such as growth, remodeling and disease.
The aim of this research was to evaluate the stability of seven potential reference genes (Actb, GAPDH, 18S rRNA, CycA, Hprt1, Ywhaz, and Pgk1) and select the most stable genes or a combination of stable genes for the purpose of normalization in IVD gene expression studying of rat and rabbit organ culture. The stability of selected reference genes under different experimental culture periods and species was analyzed using geNorm [6], NormFinder [5] and BestKeeper [16].

Results
Quantitative real-time RT-PCR Seven candidate reference genes were selected from commonly used housekeeping genes which have different biologic function (Table 1). Their primer information was summarized in Table 2.
All RNA samples were examined for their purity. The absorbance ratio at A260/A280 nm of all samples was ranged from 1.86 to 2.09, indicating all the samples were pure during the RNA extraction procedure. For all the candidate reference genes, the melt curve analyses of PCR reactions were performed. The specificity and integrity of the products were confirmed by the presence of a single peak in dissociation curve (additional file 1). The standard curves were made by serial dilutions to determine PCR efficiencies. The qRT-PCR efficiency (E) of each primer pair was ranged from 1.901 to 2.141 (90.1% to 114.1%) with linear correlation coefficient (R 2 ), making all assays suitable for quantitative analysis (Table 3). Figure 1 represents C T values of candidate reference genes from IVD organ culture of SD rat and NZW rabbit. C T values of reference genes obtained by the rat IVD were more variable with higher standard deviation than those of the rabbit IVD. To quantify the variation, C T range was calculated by maximum and minimum values through four harvesting time points. Hprt1 was shown the most invariable expression between all samples for both rat (1.09) and rabbit (0.67). In rat, CycA (1.61) and GAPDH (1.88) were ranked as second and third gene, respectively. In rabbit, on the other hand, Pgk1 (1.02) and CycA (1.04) were ranked as second and third genes, respectively. This implies that Hprt1, CycA, and GAPDH in rat and Hprt1, Pgk1, and CycA are more potential reference genes compared to the others.

geNorm analysis
The gene expression stability of seven candidate reference genes over the organ culture of the rat and rabbit IVD was analyzed using the geNorm software applications. The analysis by geNorm showed that only four genes (GAPDH, CycA, Hprt1, and 18S rRNA) in the rat IVD reached a high expression stability with low M values, below the default limit of M = 1.5 [6] ( Figure  2a). GAPDH and CycA were identified as the best pair of reference genes. On the other hand, all reference genes showed stable M value (M < 1.5), and Actb and Ywhaz were ranked as the most stable reference genes in rabbit (Figure 2b). CycA and Hprt1 were ranked as third and fourth genes, respectively. Similarly, CycA and Hprt1 were identified to be suitable for normalization in the combined analysis obtained from rat and rabbit ( Figure 2c). To determine the optimal number of reference genes necessary for accurate normalization, calculation of the pair-wise variation (V n/n+1 ) was evaluated ( Figure 2d). The combination of two reference genes was suitable for normalizing gene expression data in rat and combined species, V 2/3 (0.238) and V 2/3 (0.316) respectively, which is close to the cutoff value 0.15. On the other hand, in rabbit samples, all pair-wise variations were close to cutoff value and three reference genes would be sufficient to normalize the target genes.

NormFinder and BestKeeper analysis
Analysis on the stability of the reference genes using NormFinder and BestKeeper showed different ranking order compared to geNorm analysis (Table 4). Norm-Finder indicated Hprt1 as the best stable gene in rat and combined species, followed by 18S rRNA, CycA, and GAPDH. All reference genes in rabbit showed better stability value compared to those of rat, and CycA ranked as best control gene with a stability value of 0.227. The BestKeeper calculated variations for all the reference genes based on geometric mean of the C T values. Genes with standard deviation (SD) higher than 1 were defined as unstable. Similar to NormFinder, Hprt1 showed the lowest C T value variation in rat and combined species, and GAPDH and CycA with lower than 1 (SD) were indicated as good candidate genes. In rabbit, all reference genes except GAPDH showed low SD values and Hprt1 which was second ranking order in NormFinder was calculated to be the most stable gene with a SD value of 0.23.

Overall ranking order and selection of the best genes
All results obtained from three programs showed different ranking order according to their different calculated algorithm. To find the combination of stable reference genes, we summarized the output results which were M value (geNorm), stability (NormFinder), and standard deviation (BestKeeper) and calculated the average and ranking order in Table 5. Based on the pair-wise variation analysis by geNorm, the number of stable reference genes was applied as two for rat and combined species and three for rabbit. The overall ranking order of SD rat was Hprt1, CycA, GAPDH, 18S rRNA, Actb, Ywhaz, and Pgk1 and the combination of Hprt1 and CycA was finally selected as best control genes. Similarly, the combination of Hprt1 and CycA were determined as the most stable housekeeping genes in analyzing of combined speices. In experiments with rabbit intervertebral disc, however, the most preferred reference genes were CycA, Hprt1, and Actb.

Discussion
It has become increasingly clear in recent years that the accuracy of qRT-PCT analysis strongly depends on the choice of the normalization approach. Among current normalization approaches available, the use of housekeeping gene as an internal control is by far the most convenient to compute, provided the biological assumption of  To further improve sensitivity, normalization with multiple reference genes is also proposed.
Rat and rabbit are preferred models for studying the pathogenesis of disc degeneration. Current knowledge suggests that mechanical stimulation may be a significant contributor to the biology of the IVD. To this end, our working hypotheses states that abnormal mechanical stimulation can induce pertinent biological cascades that may degenerate a healthy normal disc in  To determine the most stable genes, data were analyzed using three different applications (geNorm, Norm-Finder, and BestKeeper) with distinct algorithms for assessing and identifying gene stability. Briefly, GeNorm software identifies stable genes based on the principle that two ideal genes have same expression ratios in all samples. The expression ratio is reported as the stability value 'M' that is essentially an average pair-wise variation of each gene in comparison to other candidate genes. NormFinder utilizes expression variation between inter and intra groups for identifying the ideal reference genes. In other words, the program can calculate not only the overall variation of candidate reference genes but also the variation between subgroups of given sample set. BestKeeper calculates reference gene's standard deviation (SD) based on raw C T values regardless of sample's efficiency. Gene with the lowest SD is considered the most stable gene. We applied these three programs for complementary analyses. As expected, ranking data reported by the programs differed considerably in terms of magnitude of the ranking parameter and order of ranking. Surprisingly (although with varying rank order), four genes (Hprt1, CycA, GAPDH, and 18S rRNA) in rat IVD tissue and five genes (CycA, Hprt1, Actb, Pgk1, and Ywhaz) in rabbit IVD tissue were presented as available reference genes. To determine a consolidated rank ordering of reference genes using data from all the output, we decided to average the outcome measures of the three programs for each reference gene since the three algorithms determined the gene ranking  based on the function of decreasing value of their outcome measures. Previously, Axtner et al. [17] identified the best combination of genes by mean rank derived from the three programs. We employed a different strategy by using the mean of three outcome measures and ranked the genes of interest as a function of their means with the lower means ranking better than their counterparts. We found this strategy to be more appropriate for the purpose of mathematical clarity.
In conjunction with GAPDH, Actb and 18S rRNA, our selection of other potential reference genes were based on demonstrated results of stability in previous investigations Xing et al. [18] demonstrated Hprt1 to be the best reference gene in rat partial hepatectomy model, which is experimental model for the study of liver regeneration. On the other hand, Ywhaz and CycA were the most stable genes in the brain tissue and asphyxia cardiac arrest model, respectively [19,20]. In the current study, combination of Hprt1 and CycA were identified to be suitable for normalization for rat and rabbit IVD culture studies. Although expression stability of reference genes may be influenced by multiple factors, our findings suggest that experimental conditions have a significant effect in this in vitro model.
GAPDH has been successfully employed as a reference gene in IVD organ culture studies [21,22] whereas our study indicates otherwise. The instability of GAPDH in our culture model could partly be the result of oxygen tension gradients in the IVD [23] that may result in regulation of HIF-1α [24]. Actb and 18S rRNA were also observed to be regulated in our culture model despite successful application in previous studies [25][26][27]. We analyzed matrix metalloproteinase-3 (MMP-3) expression which plays a major role in the disc degeneration process using Actb and combination of Hprt1 and CycA. MMP-3 expression normalized by Hprt1 and CycA showed 3.8-fold up-regulation at day 3, while Actb induced 73.7-fold up-regulation in the rat IVD tissue (additional file 2). Our comparative analysis emphasizes the critical requirement of reliable reference gene (s) to avoid erroneous and misrepresented results.
The delta-delta C T algorithm is a convenient and standard method to analyze the relative changes in gene expression. This method requires the C T values for a reference gene(s) to be reliably lower of that of a target gene. However, C T values of rabbit Hprt1 were relatively higher than those of other optimal reference genes (CycA and Actb) in this study. Although the C T values can be reduced by using increased amount of PCR templates, we first need to confirm the range of C T values of target genes for accurate analysis of relative quantification.
Accurate normalization determines the sensitivity and reproducibility of a PCR measure making the selection process of a reference gene a very crucial step in validating the gene expression tool. It is also suggested by previous investigations that single control normalization may still lead to erroneous results, urging the need to use two or more reference genes to improve sensitivity and also maintain low expression variation [6]. It is also generally advised to maintain Ct values < 30 so that the initial abundance of the target gene is considered biologically consequential. We realized the significance of multiple reference genes and have reported two most stable genes for each individual small animal species. We believe the recommendations of this study are applicable for future investigations in IVD biology that use rat and rabbit disc tissues. Also, our inference is solely concerned with the experimental conditions stated in this investigation.

Conclusions
We evaluated the stability of seven potential reference genes (Actb, GAPDH, 18S rRNA, CycA, Hprt1, Ywhaz, and Pgk1) and attempted to identify the most stable control gene (s) for normalization of qRT-PCR data from rat and rabbit organ culture. Using geNorm, NormFinder, and BestKeeper, we determined that Hprt1 and CycA are ideal reference gene candidates in both the species of interest. Actb was also found suitable for normalization in the rabbit IVD whereas its stability is questionable in the rat model. GAPDH was found to be unsuitable for normalization in both rat and rabbit IVDs under our experimental constraints. Data presented in this work is the first of its kind focusing on the intervertebral disc and may facilitate improvement in reliability and sensitivity of qRT-PCR for IVD organ culture studies.

Methods
Intervertebral disc (IVD) organ culture Young adult male Sprague-Dawley (SD) rats that were 9-weeks (280 g) old and New Zealand White (NZW) rabbits that were 18-weeks (3.6 kg) old were obtained from Harlan Sprague Dawley, Inc. (Indianapolis, IN, USA). The animals were used in accordance with a protocol approved by the University of Iowa Animal Care and Use Facilities. Under sterile condition, animals were sacrificed and lumbar IVD motion segments were dissected from consecutive levels (L1-L6). Posterior elements and soft tissues were removed and cultured in Dulbecco's modified Eagle medium (DMEM) supplemented with 14% fetal bovine serum (FBS), 50 ㎍/㎖ L-ascorbate, 100 U/㎖ penicillin, 100 ㎍/㎖ streptomycin, and 2.5 ㎍/㎖ Fungizone. After 0, 3, 7, and 14 days under standard culture conditions (37°C, 5% CO 2 ), randomly harvested IVDs were isolated from adjacent vertebral bodies and immediately frozen in liquid nitrogen. Three IVDs were pooled together and used for RNA isolation.

RNA isolation
Samples were homogenized with TRIzol ® reagent (Invi-trogen™ Life Technologies, Carlsbad, CA, USA) and total RNA was extracted by the homogenized tissues using the RNeasy Mini Kit (Qiagen, Valencia, CA, USA) according to the manufacture's instructions. Total RNA was quantified using a NanoDrop ® ND-1000 UV-Vis Spectrophotometer (Thermo Scientific, Waltham, MA, USA) at a 260 nm wavelength.

Candidate reference genes and primers for qRT-PCT
Candidate reference genes selected were classical reference genes which are most commonly used as internal control for gene expression studies (Actb, GAPDH, and 18S rRNA) and the others (CycA, Hprt1, Ywhaz, and Pgk1) based on previous reports [19,20,28]. Most primer information was obtained from previously published primer sequences. One rat (Hprt1) and five rabbit primers (Actb, GAPDH, 18S rRNA, CycA, and Hprt1) were designed using the Primer Express ® 3.0 software (Applied Biosystems, Foster City, CA, USA) based on the sequences in the database [29].

Quantitative real-time RT-PCR (qRT-PCR)
qRT-PCR was performed with the SuperScriptTM III Platinum ® SYBR ® Green One-Step qRT-PCR kit (Invitro-gen™ Life Technologies) following the instructions with slight modification. For each sample, 50 ng total RNA was used in the assay for all reference genes except rat GAPDH (25 ng), rat 18S (1 ng), and rabbit Actb (25 ng) and all samples were run in triplicate on a 96-well optical reaction plate with the ABI PRISM 7700 Sequence Detection System (Applied Biosystems) with a Sequence Detection System (SDS) software version 2.3. The PCR Reactions were prepared in a total volume of 25 ㎕ containing 1 ㎕ diluted RNA, 0.5 ㎕ forward and reverse primer (10 μM), 12.5 ㎕ 2X SYBR ® Green Reaction Mix, and 10 ㎕ RNase-free water. The conditions for the PCR were as follows: reverse transcription at 50°C for 3 min, DNA polymerase activation and RT enzyme inactivation at 95°C for 5 min, followed by 40 cycles of denaturation at 95°C for 15 sec, primer annealing at 60°C for 30 sec, elongation at 40°C for 1 min. The quantification values were obtained from the threshold cycle (C T ) number at which the increase in signal associated with an exponential growth of PCR products using SDS software. At the end of the PCR reactions, amplification specificity was confirmed by analyzing dissociation curve. To calculate PCR efficiency for each gene, six points of 2-fold serial dilution were used to build standard curve.

Data analysis
To calculate the C T values, the fluorescence threshold was manually set to 1 in SDS software and the results were directly imported into Microsoft Excel for Best-Keeper (version 1.0) [30] data input. C T values above 30 were excluded in this study. PCR efficiencies (E) were calculated from the slope of each standard curve with the equation, E = 10 (−1/slope) .
Relative quantities (Q) was then calculated from the C T values and efficiencies for geNorm (version 3.5) [31] and NormFinder (version 19.0) [32] data input with the equation, Q = E (Minimum C T −Sample C T ) , (E = 2 for 100% efficiency).
The gene expression stability was ranked from each program and the most stable reference gene or combination genes were calculated.