Identification of stable normalization genes for quantitative real-time PCR in porcine articular cartilage

Background Expression levels for genes of interest must be normalized with an appropriate reference, or housekeeping gene, to make accurate comparisons of quantitative real-time PCR results. The purpose of this study was to identify the most stable housekeeping genes in porcine articular cartilage subjected to a mechanical injury from a panel of 10 candidate genes. Results Ten candidate housekeeping genes were evaluated in three different treatment groups of mechanically impacted porcine articular cartilage. The genes evaluated were: beta actin, beta-2-microglobulin, glyceraldehyde-3-phosphate dehydrogenase, hydroxymethylbilane synthase, hypoxanthine phosphoribosyl transferase, peptidylprolyl isomerase A (cyclophilin A), ribosomal protein L4, succinate dehydrogenase flavoprotein subunit A, TATA box binding protein, and tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein—zeta polypeptide. The stability of the genes was measured using geNorm, BestKeeper, and NormFinder software. The four most stable genes measured via geNorm were (most to least stable) succinate dehydrogenase flavoprotein, subunit A, peptidylprolyl isomerase A, glyceraldehyde-3-phosphate dehydrogenase, beta actin; the four most stable genes measured via BestKeeper were glyceraldehyde-3-phosphate dehydrogenase, peptidylprolyl isomerase A, beta actin, succinate dehydrogenase flavoprotein, subunit A; and the four most stable genes measured via NormFinder were peptidylprolyl isomerase A, succinate dehydrogenase flavoprotein, subunit A, glyceraldehyde-3-phosphate dehydrogenase, beta actin. Conclusions BestKeeper, geNorm, and NormFinder all generated similar results for the most stable genes in porcine articular cartilage. The use of these appropriate reference genes will facilitate accurate gene expression studies of porcine articular cartilage and suggest appropriate housekeeping genes for articular cartilage studies in other species.


Background
With relative quantitative real-time reverse transcriptase PCR (qPCR), multiple genes across many specimens may be evaluated to measure changes in expression. However, to accurately determine the relative expression levels, and the corresponding fold changes, a reference gene is necessary. Reference genes, frequently termed "housekeeping genes," are used to normalize the expression results for differences in cDNA quantity between different specimens and thus enable comparisons between genes of interest across treatments. In order to act as a reference, a housekeeping gene's expression should remain unchanged regardless of treatment. Genes whose expression is generally unchanged with treatment conditions are most often associated with basic cellular processes such as metabolism. Our goal was to identify the most appropriate reference genes for analyses of porcine articular cartilage.
Regardless of the tissue being examined, housekeeping genes have usually been selected based on genes used in previous studies in various human tissues, and typically include beta actin (actb), beta-2-microglobulin (b2m), glyceraldehyde-3-phophate dehydrogenase (gapdh), hydroxymethylbilane synthase (hmbs), hypoxanthine guanine phosphoribosyl transferase (hprt), ribosomal protein L13a (rpl13a), ribosomal protein S18 (s18), succinate dehydrogenase flavoprotein subunit A (sdha), TATA box binding protein (tbp), and tyrosine 3-monooxygenase/ tryptophan 5-monooxygenase activation protein-zeta polypeptide (ywhaz). A variety of genes have been used in the past as housekeeping genes in cartilage studies in various species. Gapdh has been used as a housekeeping gene in studies of human, bovine, porcine, and caprine articular cartilage, including both normal and osteoarthritic (OA) samples [1,2]. Swingler et al. [3] used sdha as a reference gene in their study of human OA cartilage. These genes appeared to be selected based on literature, not selected based on evaluation of a panel of genes to identify the most stable gene. Pombo-Suarez et al. [4] evaluated nine of these same reference genes in addition to ubiquitin C in human cartilage with advanced OA and found the rarely used housekeeping genes TATA box binding protein (tbp), ribosomal protein L13a (rpl13a) and beta-2-microglobulin (b2m) to be the most stably expressed genes while they found the most commonly used genes (gapdh, actb and 18s) to be the least stable. Pombo-Suarez et al. [5] therefore recommended that tbp, rpl13a and b2m be used as housekeeping genes for human cartilage research. In a canine study of normal and OA cartilage, rpl13a and sdha were identified as the most stable reference genes [6]. The pig has been used as a model of human OA disease, cartilage repair, xenotransplantation, and gene transfer research, but no one has yet conducted a study to determine the ideal reference gene(s) for gene expression studies in porcine articular cartilage.
The reason for using a reference gene is to control for differences in the amount of starting material, efficiency of amplification, and differences in expression from cells and the overall level of transcription [7]. Therefore, selecting a stable housekeeping gene presents a circular problem: determining a stable gene when that gene is expressed differently across samples/tissues. Several methods have been developed to identify the best housekeeping gene(s) from an initial panel of potential reference genes. Three of the most commonly used methods are geNorm, BestKeeper, and NormFinder. All of these programs attempt to provide a relative measure of the stability of a panel of genes by comparing their individual stability in relation to that of the entire panel.
In geNorm [8], the average pairwise gene expression variation of each potential housekeeping gene is compared to all other evaluated reference genes. The genes that demonstrate the least variance in comparison with all other genes are ranked as the most stable genes and are therefore likely to be the best reference genes. The authors developed a Visual Basic Application for Microsoft Excel (geNorm; [8]) to carry out the analysis.
BestKeeper, developed by Pfaffl et al. [9], uses an Excel based application to determine the most stable gene from a panel of up to ten candidate genes. The geometric mean of the cycle threshold values (Ct values) for each sample across all housekeeping genes are combined together to form the BestKeeper index. Subsequently, each individual gene is compared in a pair-wise fashion via Pearson correlation coefficients to the BestKeeper index. The outcome is a ranked order of genes in terms of their stability. The highest ranked gene is the most stable. Rather than using only one housekeeping gene or the impractical method of using all potential housekeeping genes, the authors recommended the use of the best 3 or 4 genes as that provides a realistic number of housekeeping genes while still providing adequate normalization of results.
NormFinder was developed by Ohl et al. [10] and also uses an Excel based application to determine the most stable genes from a panel. This program uses a modelbased approach, where all expression values are compared via analysis of variance, and all genes and specimen results are utilized for estimation of the expected expression values. A stability measure is calculated to identify the genes that deviate the least from the calculated values [11].
Nygard et al. [12] evaluated a panel of nine genes using the geNorm approach to determine the best housekeeping genes across 17 different porcine tissues. That study included tissues such as muscle, adipose, heart, bladder, kidney, liver, skin, intestine, pancreas, bone marrow, and different portions of the brain, but no cartilage. They identified actb, ribosomal protein L4 (rpl4), tpb, and hprt as the most stably expressed housekeeping genes across the 17 tested tissues. Though cartilage was not included in the set of tissue they evaluated, their set of potential housekeepers included all of the genes previously discussed as commonly used in cartilage in other species with the exception of rpl13a and 18s. Rpl13a, like rpl4, encodes a protein of the 60S subunit of ribosomes and is still a good candidate while 18s has been shown to vary in proportion to total RNA and is therefore no longer considered a good candidate for normalization [9]. In this study we propose to determine the best housekeeping genes for use in porcine articular cartilage and to evaluate three software packages, geNorm, BestKeeper, and NormFinder for determining overall gene stability. We used the nine genes identified by Nygard and co-workers [12] as potential housekeeping genes as a starting point with the addition of peptidylprolyl isomerase A (ppia). Ppia was added because it has been used as a normalizing gene in cartilage for both OA-related [13,14] and non-OA related studies [15,16] and it exhibited no differential expression in impacted and control cartilage specimens in our previous work [17].

Methods
RNA was extracted from the articular cartilage of 40 porcine patellae obtained from an in vitro study examining gene expression changes following an applied impact injury. Patellae were subjected to one of three treatments-axial impaction, shear impaction or no impaction (non-impacted control)-and were maintained in culture for 0 (no culture), 3, 7 or 14 d. The expression of ten potential housekeeping genes: actb, b2m, gapdh, hmbs, hprt, ppia, rpl4, sdha, tbp, and ywhaz were evaluated using quantitative real-time PCR (qPCR). The relative stability of the genes was evaluated using Best-Keeper [9], NormFinder [10] and geNorm [8].

Tissue collection
Porcine knee joints were obtained from a local slaughterhouse. Patellae were sterilely removed from the joint and assigned to one of three treatment groups: control, axial impacted, or shear impacted. Patellae to be impacted were positioned in a custom holder in a hydraulic load frame (MTS MiniBionix, MTS, Minneapolis, MN). Impactions were carried out with a stainless steel impactor measuring 10mm long by 10 mm in diameter. For the axial impactions, a 2,000 Newton load was rapidly applied (loading rate of 25 mm/s) normal to the patella surface in the center of each facet. For the shear impaction a 500 Newton axial load was slowly applied (loading rate of 0.5 mm/s) followed by a rapid (200 mm/s) 10mm horizontal displacement to induce larger shear forces. Intact patellae were then placed in culture at 37°C with 5% CO 2 . Media (Delbecco's/MEM, 10% fetal bovine serum, ascorbic 2-Phosphate (25 μg/mL), penn 100 units/mlstrep 100 μg/mlamphotericin B 25 μg/mL; Gibco, Grand Island, NY) was changed daily. Patellae were maintained in culture for 0, 3, 7, or 14 d at which point 5 mm × 10 mm full thickness sections of cartilage were removed and immediately flash frozen in liquid N 2 and stored at -80°C until RNA extraction was performed. Zero day tissue was harvested on the day of impaction approximately 2 hours after impaction.

RNA extraction
Total RNA was extracted by first grinding the cartilage specimens in a mortar and pestle cooled by liquid nitrogen. The resulting powder was dissolved in Tri Reagent (Molecular Research Center Inc., Cincinnati, OH). The tissue was then homogenized in a BeadBeater W (Biospec Products, Bartlesville, OK) for 10 s at 4,800 oscillations per minute. The manufacturer's protocol was followed except the RNA was first precipitated in the presence of acetic acid and then in the presence of ammonium acetate [18][19][20]. Finally, on-column DNAse digestion was accomplished with an RNeasy kit (Qiagen, Valencia, CA). Purity of the RNA was measured on a Nanodrop spectrophotometer (Thermo Scientific, Wilmington, DE) and a sampling of samples were run on a 1% agarose gel to ensure little to no degraded RNA.
qPCR A High Capacity cDNA reverse transcription kit (Applied Biosystems Inc., Foster City, CA) was used to reverse transcribe 250 ng of total RNA per the manufacturer's protocol. Subsequently, reactions were diluted 1:10 to provide enough template for all genes to be evaluated. PCR primer sequences for the evaluated genes were obtained from Nygard et al. [12] with the exception of ppia (NM_214353.1). The ppia primers were designed with Beacon Designer software (Premier Biosoft Intl., Palo Alto, CA) from porcine gene sequences as previously described (F: 5'-GCAGACAAAGTTCCAAAGACAG-3' , R: 5'-AGATGCCAGGACCCGTATG-3') [17] spanning an intron to detect genomic contamination.
qPCR was performed in a volume of 20 μL, consisting of 1 μL of diluted cDNA, 400 nmol/L of forward and reverse primers, 10 nmol/L fluorescein, and 1X Power SYBR Green Master Mix. A three-step amplification protocol was performed in an iCycler IQ (Bio-Rad, Hercules, CA); an initial denaturation was performed with one cycle at 95°C for 7 min. Subsequently, target amplification involved 40 cycles of 30 s at 95°C, 30 s at 56°C to 62°C for annealing, then extension for 30 s at 72°C. After 40 amplification cycles, PCR products were evaluated for quality using melt curve analysis, which entailed 5 min at 72°C, 1 min at 95°C, and 1 min at 55°C. Reactions were performed in duplicate and Ct values were averaged for the replicates and negative controls were included to detect contamination.
Standard curves were evaluated for each primer by combining equal amounts of cDNA from each specimen into a pool. The pool was then diluted in serial dilutions of 1:3, 1:9, 1:27, 1:81, and 1:243. The dilutions were evaluated in triplicate by iCycler iQ Real-Time PCR Detection System Software v3.1 (Bio-Rad, Hercules CA) to calculate amplification efficiency.

Data analysis
BestKeeper, geNorm, and NormFinder were used to select the most stable genes. For the BestKeeper program, raw Ct values were entered and a BestKeeper index, which is the geometric mean of all housekeeping gene Ct values, was calculated. Pearson correlations between each individual gene and the BestKeeper index were calculated and reported as the BestKeeper correlation coefficient. Genes with the highest BestKeeper correlation coefficient were considered the most stably expressed. While there is no specific threshold for the BestKeeper correlation coefficient, Pfaffl et al. [9] recommended the use of multiple genes geometrically averaged to control for outliers. They suggested three genes were a realistic number to use for most studies while still ensuring accurate normalization [9].
The geNorm program uses normalized Ct values, where Ct values for a particular gene are normalized to the specimen with the highest expression (minimum Ct value) for that gene. The normalized Ct values (Q) are calculated via the delta-Ct formula (Equation 1).
where: Q = normalized Ct value for a given gene in the current specimen, E = PCR amplification efficiency (ranging from 1 to 2 with 100% = 2) calculated from standard curve, minCt = minimum Ct value for the gene among all specimens, and sampleCt = the Ct value of the gene for the current specimen.
In geNorm, pairwise comparisons of each gene with every other gene are performed to determine their relative stability in gene expression. Vandesompele et al. defined the stability measure M j of a given gene (j), as the mean of all pairwise variations V jk , between gene j and all other examined genes [7] (Equation 2).
where: M j = gene stability measure, V jk = pairwise variation of gene j relative to gene k, and n = total number of number of examined genes.
Lower M values represent genes with more stable expression across specimens being compared.
NormFinder also relies on Q values (Equation 1) as input, calculated from the Ct values. The program then log transforms the data and a model based approach is used with analysis of variance to calculate the expected value for each sample. The deviation of the measured value from the expected is used to calculate a stability value that ranks the genes, with the lowest value indicating the most stable [10,11].
Vandesompele et al. [7] suggested using 3 or 4 of the most stable genes for accurate normalization, using the geometric mean of the Ct values of the chosen housekeeping genes (Equation 3).
Geometric mean ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 1 a 2 . . . a n n p ð3Þ where: a = individual Ct values for the specimen's housekeeping genes and n = total number of housekeeping genes employed. The geometric mean is better able to control for outliers and abundance differences than the arithmetic mean (the sum of the individual Ct values divided by the n -the total number of values). Thus the most accurate normalization strategy is to use the geometric mean of the 3 or 4 most stable genes for normalization [7,9].

Results
Ten potential housekeeping genes for articular cartilage were evaluated in an in vitro porcine patella organ culture model that included non-impacted control, axially impacted and shear impacted tissue subjected to culture times of 0 (no culture), 3, 7, or 14 d. Two or more specimens at each impaction treatment/time point combination were evaluated using a total of 40 patellae (Table 1).
PCR amplification products were obtained for all genes but hprt, which was excluded from analysis as it displayed consistently high Ct values (greater than 35) and failed to amplify in five samples in which all other genes amplified, suggesting it is not expressed in sufficient quantity to be used as an effective housekeeping gene in these specimens. The Ct values from each of the candidate genes were input directly in the BestKeeper software [9] and were used to calculate the input values (Q values) for geNorm [8], and NormFinder [10].

Discussion
To make accurate comparisons of changes in gene expression when studying a tissue, it is important to select the best reference gene(s) for normalizing Ct values. A perfect reference gene would be expressed stably in all cells under all conditions, however, a perfect reference gene has not been found. BestKeeper, geNorm, and NormFinder provide three approaches for examining potential genes to select the most stable housekeeping genes for a given set of conditions.
We found all programs easy to use, each providing an easily accessible measure of gene expression stability in a tissue. The three programs agreed on the 4 most stable genes. The geNorm program provided a simpler more user-friendly and structured interface as it was programmed in Microsoft Visual Basic Language (VBL). This made geNorm a simpler program to use, however the equations were hidden from the user and missing values for a specimen for a particular gene could not be accepted, necessitating the removal of the entire specimen from the analysis. In addition, the user is required to calculate a Q value in geNorm, which may be an additional calculation for the user, depending upon the software associated with the qPCR instrument used.
BestKeeper is based upon an Excel spreadsheet but does not employ VBL and therefore did not have as simple a user interface. However it allowed the user to clearly see the equations used and the various steps involved in the calculation. Additionally, BestKeeper allowed for direct input of Ct values, and accepted missing Ct values for a particular gene. NormFinder was an add-in for Excel and relied upon Q values as input like geNorm. Similarly to geNorm, NormFinder could not accept missing data in the input, and the equations used in the calculation of the stability value are hidden from the user. While all programs were relatively easy to use, BestKeeper was our preferred method because of its ability to handle missing data, the ease of entering Ct values directly, and the transparency of calculations in each step. The candidate housekeeping genes that were evaluated in this study were selected from various studies examining gene expression in cartilage. The majority of these genes were also examined by Nygard et al. [12] in a study evaluating reference genes in 17 porcine tissues which showed that the ideal reference genes are tissue specific. Therefore, it is important to evaluate potential housekeeping genes for the particular tissue being utilized in a study. Because our research involves cartilage, the intent of this study was to build on the work of Nygard and colleagues to determine the most appropriate housekeeping genes specifically for porcine articular cartilage. Ppia was added because it has been used as a normalizing gene for other studies examining cartilage [13][14][15] and it exhibited no differential expression in impacted and control specimens in our previous work [17]. We found ppia, sdha, gapdh, and actb to be the most stable reference genes for porcine articular cartilage across our treatments and timepoints.
In addition to the Nygard et al. [12] study, four previous studies have evaluated reference genes for various porcine tissue. Erkens et al. [21] evaluated ten potential reference genes and found that actb, tbp and topoisomerase 2-beta were stable and that sdha appeared to be unstable in porcine backfat and muscle. Four genes were analyzed by Svobodova et al. [22] in seven porcine tissues, including heart, liver, lungs, spleen, kidney and muscle. In contrast to the results of our study, gapdh was found to be relatively unstable while hprt was found to be stable. Kuijk et al. [23] studied seven reference genes in different stages of porcine embryonic development. Of the panel of genes, gapdh, pgk1, s18 and ubc showed high stability. Nygard et al. [12], Svobodova et al. [22], Piorkowska et al. [24], and Erkens et al. [21] found tissue specific regulation of potential reference genes. Therefore our study was critical for identifying the best reference genes specifically for articular cartilage.
Previous gene expression studies in pig, cattle, goat, sheep, dog and human cartilage [1,2,25,26] have used a variety of housekeeping genes, including gapdh, sdha, s18 and actb, but these genes appeared to be selected based on what others had used in similar studies and not because they had been selected as the most stably expressed gene, as we have done here. Evaluation of appropriate housekeeping genes in human cartilage with advanced OA has been reported by Pombo-Suarez and co-workers [5] who suggest that tbp, rpl13a, and b2m be used in such studies. Tbp, rpl4 (which is similar in function to rpl13a), and b2m did not perform as well in our panel of genes with b2m generally being the least stable of the genes we examined. While we found ppia, sdha, gapdh, and actb to be the most stable, Pombo-Suarez et al. [4] found that gapdh in particular was one of the least stable genes while actb and sdha were in the middle of the pack. Differences in the selection of housekeeping genes for cartilage between Pombo-Suarez et al. and this study could represent expression differences due to species (human vs. porcine), tissue condition (in vivo normal and OA tissue vs. in vitro impacted and control), or even age (sample averages from human tissue ranged from 72 to 81 yr).
Using the pig as a model for cartilage repair studies and osteoarthritis research is quite common today. Therefore we think our results will be useful to researchers evaluating gene expression in porcine articular cartilage and as a starting point for choosing appropriate housekeeping genes in other species.

Conclusions
BestKeeper, geNorm, and NormFinder all generated similar results for the most stable genes in porcine articular cartilage. Peptidylprolyl isomerase A, succinate dehydrogenase flavoprotein, subunit A, glyceraldehyde-3phosphate dehydrogenase and beta actin should be used together by taking the geometric mean of the expression to effectively normalize expression levels for the gene of interest. The use of these appropriate reference genes will facilitate accurate gene expression studies of porcine articular cartilage and will facilitate the choice of appropriate housekeeping genes for articular cartilage studies in other species.