Selection of reliable reference genes for quantitative real-time RT-PCR in alfalfa

Real-time quantitative RT-PCR (qRT-PCR) is the most commonly used method for accurately detecting gene expression patterns. of internal control genes, β , , β , dehydrogenase GAPDH and of unknown function, in a of alfalfa ( Medicago ) samples representing different tissues and abiotic stress challenges, using geNorm BestKeeper software. The results that the eight candidate genes are inconsistently different was analyzed using selected reference genes. These results provide an experimental guideline for future research on gene expression in alfalfa using qRT-PCR.


MAINBODY
Gene expression analysis is a vital tool for developing deeper insight into the biological processes of plants. At present, real-time quantitative RT-PCR (qRT-PCR) is the most reliable method for the detection of target genes (Bustin, 2000(Bustin, , 2002, particularly for genes expressed at low abundance. Compared with classical molecular techniques, the major advantages of qRT-PCR are higher sensitivity, better reproducibility and specificity, and higher throughput (Heid et al., 1996;Gachon et al., 2004;Udvardi et al., 2008). qRT-PCR analysis requires internal control gene(s) to normalize the raw data in order to eliminate errors due to differences between samples in RNA production, RNA quality and reverse transcription efficiency (Andersen et al., 2004;Bustin and Nolan, 2004). In plants, reference gene evaluation studies have mainly been conducted in model plants and crops (Nicot et al., 2005;Czechowski et al., 2005;Expósito-Rodríguez et al., 2008;Jian et al., 2008;Libault et al., 2008;Wan et al., 2010;Yan et al., 2012). Ideally, the expression of reference genes should be stable across different organs and different stress conditions. However, the most frequently used reference genes are not always stably expressed. Numerous reports have demonstrated the lack of a "universal reference gene", and reference gene expression can vary across different tissues, different developmental stages and different environmental conditions (Lee et al., 2002;Gutierrez et al., 2008). The use of an inappropriate reference gene in qRT-PCR analysis can lead to significant deviation in data, resulting in altered findings (Suzuki et al., 2000;Bustin and Nolan, 2004;Huggett et al., 2005). Therefore, it is necessary to Edited by Shuhei Nasuda * Corresponding author. E-mail: gaohongwen@263.net select appropriate control genes according to the experimental conditions and samples in order to obtain the most accurate gene expression analysis results.
High salinity is one of the most severe abiotic stresses: it reduces agricultural production by affecting every aspect of plant physiology and metabolism. In recent decades, great progress has been made in elucidating salt stress signaling (Hasegawa et al., 2000). Abscisic acid (ABA) is a powerful signaling molecule that accumulates in response to abiotic stress and is involved in salt stress tolerance (Zhu, 2002;Zhang et al., 2006). Under stress, ABA accumulates to a high level to help plant survival by inhibiting processes such as stomatal opening and plant size expansion. On the other hand, ABA functions as a cellular signal mediating the expression of genes responsive to drought and salt stresses (Jia et al., 2002;Zhang et al., 2006).
Alfalfa (Medicago sativa) is the largest forage crop worldwide (Wang et al., 2015). Considerable progress has been made to date by genetic and molecular studies on alfalfa, and high-throughput gene expression studies have become necessary for helping us to better understand and utilize this legume forage crop. In the present study, we evaluated eight housekeeping genes as internal controls for normalizing qRT-PCR data from alfalfa tissues and plant samples challenged by abiotic stress in order to find reliable control genes for expression analysis in alfalfa.
Surface-sterilized alfalfa seeds were germinated in petri dishes at 24°C under 12-h light/12-h dark conditions. After 2-3 days, seedlings were transplanted into pots with 3:1 (v/v) vermiculite:perlite and grown in a growth chamber. The plants were irrigated with Hoagland solution (Epstein, 1972)  Total RNA was extracted using TRIzol reagent (Invitrogen, USA) according to the manufacturer's instructions. The integrity of RNA samples was assessed by agarose gel electrophoresis. Single-stranded cDNA was synthesized with Oligo(dT)18 primer (Bogdanović et al., 2013) using a Thermo Scientific RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, USA) according to the supplier's standard protocol.
A SYBR Premix Ex Taq kit (Takara, Japan) was used for PCR experiments. The PCR mixture contained 10 μl of 2×SYBR Premix Ex Taq II, 0.4 μl of ROX Reference Dye II (50×), 1 μl (10 pmol) each of gene-specific primers, and 1 μl of diluted cDNA in a final volume of 20 μl. The realtime PCR reactions were performed with an ABI Prism 7500 Sequence Detection System (ABI, USA). Thermal conditions were 95°C for 30 s, followed by 40 cycles at 95°C for 5 s and annealing and extension at 60°C for 34 s. Quantification analysis was performed by the comparative cycle threshold (Ct) method. All reactions were performed in triplicate for each of three biological replicates represented by pools of tissue collected in different treatments. Relative expression level data of the eight candidate reference genes were imported into the geNorm (version 3.5) (Vandesompele et al., 2002; http://medgen. ugent.be/~jvdesomp/genorm/) and BestKeeper (Pfaffl et al., 2004; http://www.gene-quantification.de/bestkeeper.html) analysis tools, which were used as described in the product manuals.
The specificity of the PCR primers and the expression levels of these eight candidate internal reference genes were tested using a set of 16 samples from three different groups: the tissue series, which included samples from four different organs; the salt treatment series, which included six samples representing different NaCl stress challenge times; and the ABA treatment series, which included six samples representing different ABA stress challenge times. All eight primer pairs produced a single PCR product of the expected size, as determined by a dissociation curve analysis performed by the quantification PCR system after 40 cycles of amplification and agarose gel electrophoresis ( Supplementary Fig. S1). In addition, the PCR products were cloned into a T/A cloning vector and sequenced. The orthologous sequences showed high similarity with the original sequences. Amplification curves for these eight genes were generated and grouped across all the samples, and Ct was determined at a fluorescence value of 0.2. The highest Ct value was for β -Actin (29.3), representing the lowest level of transcription, and the lowest value was for 18S rRNA (13.0), representing the highest transcript level. The other six genes had average Ct values ranging from 15.7 to 26.7. Each reference gene showed variation in Ct values across different sample types ( Supplementary Fig.  S2). According to these Ct data, it is clear that the expression levels of all these candidate genes, except 18S rRNA, are inconsistent under different experimental conditions.
geNorm 3.5 software was used to rank the tested reference genes based on their expression stability value, M. Genes with lower expression stability have higher M values, whereas the most stable gene has the lowest M value (Fig. 1). The tissue samples were analyzed as four different data sets: the first included all the samples, the second the NaCl treatment series, the third the ABA hormone treatment series, and the fourth the tissue type series. The average M value for the eight candidate reference genes across all the samples was lowest for Msc27 and 18S rRNA, followed by UBC2 and EF-1α, whereas that for TUB was the highest (Fig. 1A). Msc27and 18S rRNA also showed the most stable expression levels under salt stress conditions, followed by EF-1α and ACT2, with β -Actin the least stable gene (Fig. 1B). In the ABA treatment series, ACT2 and 18S rRNA had the lowest M value and GAPDH had the highest (Fig. 1C). For the tissue-type series, Msc27 and 18S rRNA were again the most stable genes with the lowest M value (Fig. 1D). The M values for the eight candidate reference genes in different tissues and under different treatment conditions were all below the default limit of 1.5 in geNorm; as noted above, a lower value of M indicates more stable expression (Fig. 1). Remarkably, the M values of Msc27 and 18S rRNA were the lowest in most of the series (all samples, NaCl stress conditions and different tissue types), indicating that Msc27 and 18S rRNA are stably expressed.
To obtain accurate results, a single reference gene is usually not enough (de Jonge et al., 2007). geNorm provides a method that can calculate pairwise variation (V n /V n+1 ), to determine the optimal number of internal reference genes to be used in an expression analysis. Generally, 0.15 is used as a cutoff value, below which the inclusion of an additional reference gene is not required. A value higher than the cutoff value means that the added reference gene (n+1) has a significant effect and should be included for calculation of a reliable normalization factor (Vandesompele et al., 2002;Czechowski et al., 2005). The pairwise variation for the salt stress series showed a V 3/4 value of 0.123 (< 0.15), indicating that the gene expression normalization needs three reference genes, EF-1α, Msc27 and 18S rRNA (Figs. 1B and 2). For the ABA treatment series, the

β-Actin TUB GAPDH UBC2 ACT2 EF-1α Msc27 18SrRNA
Average expression stability M B V 3/4 value was 0.130, suggesting that the three most stable candidate reference genes, ACT2, 18S rRNA and Msc27, would be sufficient for normalization of the mRNA transcription level analysis data. For different tissue types, Msc27 and 18S rRNA were the most appropriate reference genes since V 2/3 was 0.143, below the cutoff value of 0.15 (Fig. 2). When all samples were taken together, a decrease in the pairwise variation appeared with the inclusion of a fifth gene. The V 4/5 value was 0.118 (Fig.  2), indicating that at least four reference genes should be included for gene expression studies. The qRT-PCR results were further analyzed using the program BestKeeper. Variation was calculated for all the reference genes in all the different samples based on the variables of standard deviation (SD [±Ct]) and percentage covariance (CV [%Ct]). Reference genes with SD values greater than 1 were considered as inconsistent and were excluded (Wan et al., 2010). Based on this analysis, 18S rRNA had the lowest SD and CV values among all eight genes, indicating that it is the most stable reference gene (Table 1). ACT2 and Msc27 were ranked just after 18S rRNA and also had low SD and CV values.
TUB, GAPDH and β -Actin had SD values higher than 1, and were therefore unsuitable for use as reference genes in alfalfa under the conditions in this study.
WRKY is an important transcription factor gene which is involved in both biotic and abiotic stress responses (Luo et al., 2005;Jiang and Deyholos, 2009). Using RACE technology, we isolated a WRKY I type gene (MsWRKY33, accession no. KP642130) from alfalfa and found that it was up-regulated by salt stress. Therefore, the relative expression of MsWRKY33 was analyzed using the three most stably expressed reference genes (Msc27, EF-1α, 18S rRNA) to test the reference gene selection result and to study this gene's expression profile in response to NaCl stress. Similar expression patterns were observed when these three reference genes were used for normalization of the MsWRKY33 gene, indicating their ability to yield reliable results in relative gene expression analysis. As shown in Fig. 3, salt stress dramatically induced the expression of MsWRKY33: expression was rapidly induced after 4 h of stress treatment, and rose continuously to a peak at 12 h. These results indicate that the MsWRKY33 gene is involved in the salt stress response of alfalfa.
GAPDH has been reported to be a suitable normalization reference for measuring gene expression in different tissues/organs of sugarcane (Iskandar et al., 2004) and in hop (Nagel et al., 2008). However, in our study, the M values of GAPDH and β -Actin were higher than those of most of the other candidate genes in all four test series, suggesting that GAPDH and β -Actin are not suitable ref- Fig. 2. Pairwise variation analysis to determine the optimal number of reference genes for normalization of qRT-PCR data for all samples, NaCl stress series, ABA treatment series and different tissue type series. Arrows indicate the optimal number of genes for normalization.   (Bustin, 2002). Therefore, a reference gene expressed stably in one organism must be validated before its use in qRT-PCR on samples from another organism under certain conditions. Some researchers believe that total RNA can be used to normalize qPCR data (Huggett et al., 2005). However, this method cannot normalize errors that arise during reverse transcription, PCR amplification and other experimental procedures. Furthermore, it relies heavily on the accuracy and reliability of RNA quantification. Generally speaking, spectrophotometer readings have several shortcomings, such as low sensitivity, and can easily be affected by residual DNA in the cuvette. A fully automated biochemical analyzer, such as the Beckman Au680, is superior to a spectrophotometer, but because its cost is often prohibitive, the use of reference genes remains the first choice for most researchers. In this study, we used two different statistical tools, geNorm and BestKeeper, to evaluate the eight candidate reference genes. Using these software tools, any number of housekeeping genes can be evaluated in any samples from different tissue types or conditions, and two or more reference genes can be selected rather than the traditional single gene to normalize RT-PCR results. This process will reduce system deviation and help to achieve reliable results. Therefore, this approach is of great significance in obtaining accurate quantitation of gene expression, particularly when examining subtle expression differences between genes. Use of this approach requires that the following conditions be met: 1) extraction of high-quality total RNA; 2) use of an equal amount of total RNA when performing reverse transcription under uniform experimental conditions; and 3) minimization of experimental errors that occur during the reverse transcription and PCR process.
Controlling these factors will make the results of analysis more accurately reflect the regulation of housekeeping genes and reduce the influence of errors when selecting reference genes.
According to the BestKeeper analysis results, the four reference genes that had the least gene expression variation [SD (± x-fold)] were 18S rRNA (1.20), ACT2 (1.29), Msc27 (1.40) and EF-1α (1.41). Results from geNorm showed that the same four genes were needed for normalization across all the samples (Fig. 2). These results indicate consistency between the geNorm and BestKeeper software. However, in the geNorm analysis results, Msc27 was ranked second to 18S rRNA in expression stability, whereas in the BestKeeper analysis, ACT2 was ranked second. This difference may be due to the distinct statistical algorithms used by these two methods.
From the qRT-PCR and stability analysis results, we found an interesting phenomenon: 18S rRNA and Msc27 have the lowest M values, and their expression levels were the highest among all candidate genes, whereas TUB and β -Actin showed higher M values than other genes and their expression levels were also lower than others. Therefore, it can be supposed that, at least in this study, highly expressed genes have more stable expression than weakly expressed genes, and vice versa.
In summary, we evaluated eight candidate reference genes to select suitable control genes for normalization of real-time PCR data in alfalfa. Using two analysis methods, our results indicate that EF-1α, ACT2, Msc27 and 18S rRNA are suitable reference genes under all of the conditions tested in this study. Msc27 and 18S rRNA are the most stably expressed and should be considered suitable reference genes in different tissues of alfalfa. Under ABA treatment, ACT2, Msc27 and 18S rRNA are the most appropriate control genes. Finally, under different NaCl conditions, reliable reference genes are EF-1α, Msc27 and 18S rRNA. Of note, other frequently used reference genes such as GAPDH or β -Actin were found to be less reliable and are therefore unsuitable for use as reference genes. Such results are critical for more accurate and reliable normalization of gene expression research on this important forage crop plant in future.
We thank Dr. Ian Smith and the anonymous reviewers for their comments on the manuscript. This study was supported financially by the National Natural Science Foundation of China