Selection of Suitable Reference Genes for Normalization of Quantitative Real-Time Polymerase Chain Reaction in Human Cartilage Endplate of the Lumbar Spine

The quantitative real-time polymerase chain reaction (qRT-PCR) is one of the most widely used methods to study gene expression profiles, and it requires appropriate normalization for accurate and reliable results. Although several genes are commonly used as reference genes (such as GAPDH, ACTB, and 18S rRNA), they are also regulated and can be expressed at varying levels. In this study, we evaluated twelve well-known reference genes to identify the most suitable housekeeping gene for normalization of qRT-PCR in human lumbar vertebral endplate with Modic changes, by using the geNorm, NormFinder, and BestKeeper algorithms. Our results showed that the rarely-used SDHA was the most stable single reference gene, and a combination of three, SDHA, B2M, and LDHA, was the most suitable gene set for normalization in all samples. In addition, the commonly-used genes, GAPDH, ACTB and 18S rRNA, were all inappropriate as internal standards. The rankings of reference genes for the three types of Modic change differed, although SDHA and RPL13A uniformly ranked in the first and last position, respectively. Further simulated expression analysis validated that the arbitrary use of a reference gene could lead to the misinterpretation of data. Our study confirmed the necessity of exploring the expression stability of potential reference genes in each specific tissue and experimental situation before quantitative evaluation of gene expression by qRT-PCR.


Introduction
Signal intensity changes in the vertebral endplate and subchondral bone marrow on magnetic resonance imaging, also known as Modic changes (MCs), are often observed in patients with degenerative spinal diseases [1]. Three different types have been described. Type I changes are hypointense on T1-weighted imaging (T1WI) and hyperintense on T2-weighted imaging (T2WI) and indicate edema and hypervascularity in the lesions as confirmed histologically. Type II changes are hyperintense on T1WI and isointense or hyperintense on T2WI and reflect fatty replacement of the red bone marrow. Type III changes are hypointense on both T1WI and T2WI and represent subchondral bone sclerosis. MCs are strongly associated with intervertebral disc degeneration and low back pain [2][3][4]. Identification of the mechanisms and factors involved in the progression of MCs is of great importance for clinical interventions to repair or retard the development of MCs. To exploit the mechanisms of MCs, it is necessary to understand the pathophysiological changes of the vertebral cartilage endplate with MCs at the molecular level.
Comprehensive evaluation of gene expression patterns is important for understanding the biological processes occurring in vertebral cartilage endplate with MCs.
Gene expression analysis is widely used in many fields of biological research [5,6]. Gene expression analysis in cartilage has been satisfactorily performed in recent years since the major problem of poor ribonucleic acid (RNA) content in human cartilage is compensated by improving the RNA yield and quality and via cDNA amplification by in vitro transcription [7][8][9][10]. Quantitative real-time polymerase chain reaction (qRT-PCR) is currently one of the most precise and frequently-used methods to study the expression profiles of genes; it can quantify both the absolute and relative amounts of a gene's RNA [11,12]. Relative quantification is more common in qRT-PCR, and its accuracy, reliability, and reproducibility are highly dependent on the choice of suitable internal controls within each sample to normalize experimental variations [13]. The use of reference genes can correct biases caused by variations in RNA input, or reverse transcription efficiency, or amplification efficiency.
An ideal reference gene is presumed to be expressed at a constant level in all tissues and cells and under different experimental conditions. Although several genes are commonly used as controls (such as GAPDH, ACTB, and 18S rRNA), they are also regulated and can be expressed at varying levels [11]. The use of inappropriate reference genes may weaken the detection sensitivity of the target genes and even lead to wrong results [14,15]. Therefore, the expression stability of potential reference genes should be explored in each specific tissue and experimental situation. And the use of more than one reference gene is required for high-quality data, as suggested by many investigators [13,16].
The suitability of reference genes in human osteoarthritic articular cartilage (hip and knee) has been evaluated by Pombo-Suarez et al. [17]. Their analysis showed that the expression levels of GAPDH, ACTB, and 18S rRNA varied between samples. On the contrary, the rarely-used TBP, RPL13A, and B2M were the most stable and it was necessary to use several of these together to obtain the best results. So far, no appropriate reference genes have been identified in human lumbar cartilage endplate with MCs. In this study, we investigated the expression levels of 12 well-known reference genes to identify those most suitable for normalization of qRT-PCR in this tissue.

Human Lumbar Cartilage Endplate Sample Collection
Lumbar cartilage endplate specimens were obtained from patients who underwent lumbar interbody fusion for various spinal diseases in the Orthopaedic Surgery Department. All patients received a magnetic resonance imaging examination to determine the presence or absence of MCs in the lumbar vertebral endplate. Patients with MCs at the operated segment were recruited into the MC group, and those without MCs formed the control group. Finally, 12 MC and 12 non-MC samples, age-and sex-matched, were used. The Medical Ethics Committee of Sir Run Run Shaw Hospital, School of Medicine, Zhejiang University approved this study. All participants gave their written informed consent prior to participation. Patient characteristics are summarized in Table S1.
The cartilage endplate specimens were separated from disc and bone tissues immediately after they were harvested from the intervertebral space. They were then snap-frozen in liquid nitrogen and stored at 280uC until the extraction of total RNA.

RNA Extraction Quality Control and cDNA Synthesis
RNA extraction from lumbar cartilage endplate was performed following the method of Untergasser [18]. Samples were first cut into small pieces under sterile conditions and ground in liquid nitrogen. Total RNA was then extracted using Trizol (Invitrogen) and further purified with an RNeasy Mini Kit (Qiagen) according to the manufacturer's protocol with minor modifications. RNasefree DNase (Qiagen) was used twice to remove any trace of genomic DNA.
The concentration and purity of the isolated RNA were estimated in triplicate using a NanoDrop 2000 spectrophotometer (Thermo Scientific). Samples with concentrations $50 ng/ml and optical density absorption A 260 /A 280 between 1.8 and 2.1 were taken for cDNA synthesis. The integrity of RNA samples was confirmed by electrophoresis on 2% Sybr Green agarose gel (Invitrogen) as indicated in the MIQE (Minimum Information for Publication of Quantitative Real-Time PCR Experiments) guidelines [19].
Total RNA (200 ng) was reversed-transcribed to first-strand cDNA using a Reverse Transcription System (Promega) in a total volume of 20 ml, according to the manufacturer's instructions.

Selection of Candidate Reference Genes and Primer Design
Twelve candidates were chosen based on their common use as reference genes and previous screening from microarray expression data of osteoarthritic cartilage [20,21]. These genes represented several distinct functional classes so as to reduce the chances of co-regulation and false-positive reference gene selection.
For all genes, primer pairs were designed using Primer 3 ver. 0.4.0 (http://frodo.wi.mit.edu/primer3/) and then checked for the absence of stable hairpins and dimers using Oligo 5.0 (Molecular Biology Insights, Cascade, CO). All the primers were designed to be close to the 39 end of the RNA sequence, and to be located on different exons to avoid genomic DNA contamination, except the primer for 18S rRNA ( Table 1). The primers were synthesized by Sangon Biotech Co., Ltd (Shanghai, China).

Quantitative Real-time PCR
qRT-PCR was performed using the GoTaqH qPCR Master Mix kit (Promega). Reactions were run in triplicate on a 7500 Real Time PCR System with 96-well plates (ABI). SYBR Green was used to detect dsDNA synthesis. Each reaction was performed in 20 ml containing 2 ml cDNA, 10 ml Gotaq Master Mix, 0.4 ml upstream and 0.4 ml downstream PCR primers (0.2 mM), and 7.2 ml nuclease-free water. Amplification was performed with an initial holding period at 95uC for 2 min, followed by a two-step PCR program consisting of 95uC for 5 s and 60uC for 34 s for 40 cycles. After 40 cycles, a melting analysis was performed by heating the amplicon from 60uC to 95uC. A reverse transcriptase negative control was included to ensure the absence of genomic DNA contamination, and the no-template control was also run to exclude contamination or dimer formation for each primer pair. The amplification specificity was confirmed by melting curve analysis and agarose gel electrophoresis of the products. For each primer pair, a series of 10-fold of three dilutions of cDNA (10-to 1,000-fold dilution) were made to generate a standard curve. The PCR amplification efficiency (E) was determined by the slope of the standard curve: E(%) = (10 [21/slope] 21)6100%.

Statistical Analysis
The Ct values for each sample were compiled and run through the Microsoft Excel-based software programs, geNorm (ver. 3.5) [16], NormFinder (Ver 0.953), [13] and BestKeeper [22]. These statistical algorithms were used to evaluate the stability of candidate reference genes, and then the overall ranking of the 12 candidate reference genes was determined according to the method described by Chen et al. [23].
For geNorm and NormFinder, data were analyzed by transforming raw Ct values into relative quantities using the DCt method. The lowest Ct value was subtracted from the raw Ct values of qRT-PCR for each gene to give the DCt value. The equation E 2DCt was applied to each data point. Therefore, all data were expressed relative to the expression of the most highlyexpressed gene. BestKeeper analysis was based on the raw Ct values.

Amplification Specificity and Efficiency and Expression Levels of 12 Candidate Reference Genes
Agarose gel electrophoresis and melting-curve analysis gave a single product of the expected length for each candidate gene. No non-specific amplicons or primer dimers were detected in the notemplate control, and the absence of signals in the reverse transcriptase negative control suggested no genomic DNA contamination. All PCR assays showed efficiency values between 95.1% and 104.7% (Table 1).  The expression levels of these 12 reference genes varied widely with Ct values ranging from 17.6 (18S rRNA) to 33.5 cycles (ACTB), and most of the Ct values were between 24 and 33 cycles. 18S rRNA was the most abundantly transcribed with a mean Ct value of 20.2 cycles. SDHA, GAPDH, and HRPT1 were moderately expressed with most of the Ct values between 26 and 31 cycles. ACTB showed the lowest level of expression with a mean Ct value as high as 32.1 cycles (Figure 1).

Expression Stability of 12 Candidate Reference Genes
GeNorm analysis. GeNorm is a program designed to analyze the expression stability of candidate reference genes on the assumption that the ratio of the expression level of two ideal reference genes is constant in all samples. And the average expression-stability M value for each investigated gene is calculated with the average of pairwise variations, according to which the expression stability of all reference genes is ranked. Genes with the lowest M value are the most stable, and a value of 1.5 is recommended as the cut-off for the selection of qRT-PCR reference genes by Vandesompele et al. [16]. In our analysis, when the results from all 24 samples of cartilage endplate were combined, the genes with the smallest M value were SDHA and B2M (0.45), therefore these were the most stable genes. The sequence from most to least stable was SDHA, B2M, LDHA, TBP, GUSB, GAPDH, IPO8, HMBS, ACTB, HPRT1, 18S rRNA, and RPL13A, with the largest M value being 0.97 for RPL13A (Figure 2A). In subgroup analysis, SDHA and LDHA, SDHA and TBP, and SDHA and LDHA were the most stable genes for MC types I, II, and III, respectively, and RPL13A was uniformly the most unstable ( Figure 2B-2D). The M values of the 12 genes in all samples and subsets were less than 1, below the default limit of M #1.5, indicating relatively high expression stability (Figure 2).
The geNorm algorithm also calculates the pairwise variation (V n /V n+1 ) between two sequential normalization factors NF n and NF n+1 to determine the optimum number of reference genes. As a general rule, the stepwise inclusion of reference genes is performed until V n /V n+1 drops below the theoretical threshold of 0.15, when the benefit of adding an extra gene (n+1) is limited for accuracy normalization [17,24]. And the use of at least the three most stable reference genes is recommended [16]. In this study, the pairwise variations V 3/4 for total and subgroup analyses were all below 0.15, therefore the addition of the fourth-best gene to the gene set composed of the three best ones was not needed (Figure 3).
NormFinder analysis. NormFinder is a model-based algorithm, which identifies the most stable reference genes based on combing samples into groups. The main goal of this approach is to calculate the inter-and intra-group variation of the candidate reference genes and then combine both results into a stability value M. Genes with the lowest value are considered to be the most stable [13].
The calculated stability values of the 12 genes are listed in Table 2. According to NormFinder, the most stable reference gene in all samples was LDHA, with a M value of 0.102, followed by SDHA, B2M, GUSB, GAPDH, IPO8, HMBS, ACTB, TBP, 18S rRNA, HPRT1, and RPL13A. LDHA was also top-ranked for those with MC types I and III, while B2M was the most stable for those with MC type II. SDHA was among the three most stable genes and RPL13A was the most unstable in the total sample set as well as in the three subsets, as for the geNorm algorithm.
BestKeeper analysis. BestKeeper analyzes the expression stability of reference genes using raw Ct values. Gene expression variation is determined by the calculated standard deviation (SD) and coefficient of variance (CV) for all candidate reference genes based on the whole data set of their Ct values [22]. Those with the lowest CV6SD were identified as the most stable genes. Those with SD values .1 were considered to be inconsistent and were excluded. In the current analysis, SDHA and ACTB had the smallest CV 6 SD of 2.61460.776 and 2.75760.884, respectively, thus being the most stable reference genes in all samples. This was different from the results produced by geNorm and NormFinder, in which ACTB ranked poorly (8/12 and 9/12, respectively). In the three subsets, SDHA was the most stably expressed and 18r RNA or RPL13A the least stably expressed (Table 3).
Final ranking of candidate reference genes. Since the discrepancies in expression stability of candidate reference genes among the algorithms, a method taking into account all the three sets of results was applied to calculate the final ranking. Specifically, the geometric means of the three ranking numbers produced by geNorm, NormFinder, and BestKeeper were calculated for each candidate reference gene; those with the smallest geometric means were considered to be the most stable [23]. As a result, SDHA was the most stable single gene in all samples, and SDHA, B2M, and LDHA comprised the optimal reference gene set. Although the rankings of reference genes for the three types of MC were different, SDHA and RPL13A uniformly ranked in the first and last position, respectively ( Table 4).

Effect of Choice of Reference Genes
To validate the importance of selecting the appropriate reference genes for normalization, a simulated expression analysis was performed using data from samples with MC types I and III. In these two types of specimen, SDHA and LDHA consistently turned out to be the most stable genes, and RPL13A the most unstable. The simulation was conducted taking LDHA as target and SDHA and RPL13A as reference genes (Figure 4). The relative fold expression of LDHA was calculated and normalized to the lowest value: 2.5460.57 in samples with MC type I and 1.0060.21 in those with MC type III; the different between these two subgroups was significant (P = 0.012). When the expression level of LDHA was normalized to SDHA, it was higher in the subgroup with MC type I than that with type III (1.2760.08 versus 1.0060.06, P = 0.009). However, LDHA was non-significantly less expressed in the former than the latter group when using RPL13A as the reference (1.0060.89 versus 7.0968.20, P = 0.219) ( Figure 4). Therefore, the difference of gene expression levels between the subgroups could be masked by using unsuitable reference genes.

Discussion
The purpose of this study was to choose appropriate internal controls to ensure credible evaluation of gene expression levels in human lumbar cartilage endplate with MCs. All the 12 candidates were selected from previous investigations, and they were reported to have relatively stable expression levels [13,16,17,25]. Our studies showed that the most stable gene in all samples was SDHA, which is among the less commonly used. And a combination of three reference genes was recommended, i.e., SDHA, B2M, and LDHA, based on a comprehensive consideration of the results of all algorithms. The application of multiple references is beneficial for normalization [13,17,26].
We compared our results with those of a trial which selected reference genes for normalization in human osteoarthritic hip and knee cartilage [17]. B2M was identified as one of the best reference genes in both studies. The surprising result in this analysis was that GAPDH as well as ACTB and 18S rRNA performed poorly, and this has also been reported for osteoarthritic articular cartilage [17].
These genes are often given preference in many studies [16,27]. Our results, however, further confirmed the necessity of assessing the reference genes in each tissue and experimental condition. Yet there were discrepancies between our results and the published data for hip and knee cartilage. They concluded that RPL13A is among the most stable reference genes, while SDHA is not. On the contrary, we found that SDHA was the most stable while RPL13A the most unstable. Thus, in tissues from different anatomical regions, hip and knee joint cartilage versus lumbar vertebral endplate cartilage, the stability of reference genes is different.
Interestingly, the results for the three types of MC differed. For instance, IPO8, which was one of the most stably-expressed genes in samples with MC type I, was not ranked so high in MC types II and III. Similarly, TBP was relatively stable in samples from MC types II and III, however, it was ranked lower in MC type I. Therefore, even for a tissue with different subtypes of pathologic change, previous testing of the stability of reference genes is required.
The geNorm, NormFinder, and BestKeeper are now widely used for the selection of stable reference genes with invariable expression from a set of candidate genes. GeNorm analyzes the expression stability of the tested genes in all samples, and ranks them according to a stability measure [16]. In contrast, NormFinder evaluates the expression stability of each single reference gene independently, and takes into account intra-and inter-group variations for normalization [13]. While BestKeeper analyzes the stability of candidate reference genes depending on the standard deviation and coefficient of variance of their Ct values. Together they seek genes with stable expression levels, by either the relative stability with reference to other candidates, or to clinically-relevant groups, or to their intrinsic degree of variation [13,16,22]. The use of these programs provides complementary information. Because of the distinct statistical algorithms used by these three programs, it was not surprising that they gave somewhat different results. Nevertheless, there was general agreement; SDHA was in the top two positions and RPL13A in the last two in all samples and subsets for these three programs.
There were some limitations in our study. First, we only included a limited number of candidate genes, and it is likely that some other genes may be better internal references for human lumbar endplate cartilage with MCs. Second, our results only apply directly to vertebral cartilage endplate with MCs in the lumbar region. It is unclear how well our results could be extended to other regions of vertebral cartilage endplate, i.e., the cervical and thoracic regions.