Selection of Reference Genes for qRT-PCR Analysis of Gene Expression in Stipa grandis during Environmental Stresses

Stipa grandis P. Smirn. is a dominant plant species in the typical steppe of the Xilingole Plateau of Inner Mongolia. Selection of suitable reference genes for the quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) is important for gene expression analysis and research into the molecular mechanisms underlying the stress responses of S. grandis. In the present study, 15 candidate reference genes (EF1 beta, ACT, GAPDH, SamDC, CUL4, CAP, SNF2, SKIP1, SKIP5, SKIP11, UBC2, UBC15, UBC17, UCH, and HERC2) were evaluated for their stability as potential reference genes for qRT-PCR under different stresses. Four algorithms were used: GeNorm, NormFinder, BestKeeper, and RefFinder. The results showed that the most stable reference genes were different under different stress conditions: EF1beta and UBC15 during drought and salt stresses; ACT and GAPDH under heat stress; SKIP5 and UBC17 under cold stress; UBC15 and HERC2 under high pH stress; UBC2 and UBC15 under wounding stress; EF1beta and UBC17 under jasmonic acid treatment; UBC15 and CUL4 under abscisic acid treatment; and HERC2 and UBC17 under salicylic acid treatment. EF1beta and HERC2 were the most suitable genes for the global analysis of all samples. Furthermore, six target genes, SgPOD, SgPAL, SgLEA, SgLOX, SgHSP90 and SgPR1, were selected to validate the most and least stable reference genes under different treatments. Our results provide guidelines for reference gene selection for more accurate qRT-PCR quantification and will promote studies of gene expression in S. grandis subjected to environmental stress.


Introduction
Stipa grandis P. Smirn. is a C3 perennial bunchgrass [1] that is distributed widely in eastern Eurasian steppes and the middle Eurasian steppe zone [2]. This species is dominant in the typical steppe of the Xilingole Plateau of Inner Mongolia [1]. The S. grandis steppe acts as a natural a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 green barrier and plays an important role in sandstorm prevention [3]. To date, studies of S. grandis have mostly focused on ecological aspects, such as the relationship between the photosynthetic capacity of S. grandis and the soil moisture gradient [4], the effects of grazing intensity and duration on the biomass of S. grandis [5], and the effects of global warming and elevated CO 2 on plant growth [6]. The physiological responses of S. grandis to high temperature and drought have also been reported [2]. Several genetic diversity studies have been carried out on S. grandis related to grazing and location [1,3,7]. Transcriptome analysis was carried out to reveal the gene expression profiles of overgrazed and non-grazed S. grandis in a typical grassland; the results showed that many stress-related genes had differential expression under long-term grazing conditions [8]. However, the mechanism underlying S. grandis's responses to grazing and related stresses remain unclear. To further understand the stress response mechanism of S. grandis, details of the expression patterns of more stress-related genes are necessary.
Quantitative real-time reverse transcription PCR (qRT-PCR) is an efficient tool for expression analysis because of its high-throughput, accuracy, sensitivity, and ease of use compared with traditional methods, such as northern blotting and semi-quantitative reverse transcription-PCR [9][10][11]. The inevitable errors or experimental deviations caused by sample preparation and data analysis can affect the reliability of the results; therefore, it is essential to normalize the data using appropriate internal reference genes [12]. The stability of the reference gene's expression is crucial for the accurate normalization of target gene expression following qRT-PCR quantification [13]. Extensive studies on reference gene stability and selection have been reported in various plants for different tissues and development stages, and under diverse stresses [14][15][16][17][18][19][20][21][22]. Traditionally, reference genes such as 18S rRNA, GAPDH (glyceraldehyde-3-phosphate dehydrogenase), ACT (β or γ actin), TUB (α or β tubulin), EF1α (elongation factor 1α), and UBQ (poly-ubiquitin) have been used most frequently as internal control genes [13,17]. Recently, taking advantage of a variety of databases comprising data from microarrays, expressed sequence tags, and transcriptomes, other stably expressed genes have been used as novel reference genes for qRT-PCR, such as CUL (Cullin), UBCP (ubiquitin carrier protein) [17], CAP (adenylyl cyclase-associated protein), SKIP1 (ASK-interacting protein 1) [22], F-BOX (F-box/kelch-repeat protein) [23], PP2A (protein phosphatase 2A) [24], SAMDC (s-adenosyl methionine decarboxylase) [12,25], and SNF (sucrose non fermenting-1 protein kinase) [26]. However, the transcript levels of these reference genes vary among different conditions to some extent [20,27]. Thus, the selection of ideal reference genes to normalize target gene expression is important to improve the reliability of qRT-PCR results.
In this study, 15 candidate reference genes that could be used potentially as internal control genes to normalize gene expression were selected. The expression stabilities of these reference genes under various abiotic stresses and three hormone treatments were assessed using four algorithms (geNorm, NormFinder, Bestkeeper, and RefFinder). According to our analysis, we were able to recommend the most suitable reference genes to quantify gene expression in S. grandis under various experimental conditions.

Hormone and abiotic stress treatments
Four-week-old plants were used in the following experiments, and the different treatment conditions were mainly achieved according to methods described previously [12,22,28,29]. Before treatment with hormones, salinity, high pH, and drought, S. grandis plants were removed from the soil gently, washed clean with tap water, and then immediately subjected to the different treatments. For plant hormone, high pH and salinity treatments, the plants were dipped separately in solutions of 100 μM salicylic acid (SA), 100 μM jasmonic acid (JA), 100 μM abscisic acid (ABA), high pH (pH 10), and 300 mM NaCl. For the drought treatment, the plants were placed on filter paper and incubated in a growth chamber. For the wounding treatment, plants growing in pots were wounded using tweezers and placed in a growth chamber. For cold and heat treatments, plants growing in pots were maintained at 4˚C or 42˚C, respectively. Untreated plants were used as the control. Shoot samples (from four to five plants) were taken at 1, 3, 6, 12, 24, and 48 h time points for all treatments, and an additional 0.5-h time point for drought and wounding treatments, and were frozen immediately in liquid nitrogen and stored at −80˚C until RNA extraction. Three independent biological experiments were performed.

RNA isolation and cDNA synthesis
Total RNA was extracted from the samples using the Trizol regent. The purity and integrity of the RNA samples were estimated by calculating their A260/280 absorbance ratios and by agarose gel electrophoresis analysis. Total RNA (0.5 μg), pretreated with RNase-free DNase I (Takara, Dalian, China), was reverse transcribed into cDNA using an M-MLV reverse transcriptase kit (Takara) and oligo (dT) 18 primers (Takara) following the manufacturer's instructions. The generated cDNAs were diluted 20-fold for qRT-PCR analysis.

Reference gene selection and primer design
The sequences of the candidate reference genes used in this study were selected based on our previous transcriptome analysis [8], and included EF1 beta, ACT, GAPDH, SamDC, CUL4, CAP, SNF2, SKIP1, SKIP5, SKIP11, UBC2, UBC15, UBC17, UCH, and HERC2. According to these sequences (S1 Text), specific primers for qRT-PCR were designed using the Premier5 software, with a melting temperature (T) of 59-61˚C, a length of 18-25 bp, an amplicon product size of 75-160 bp, and a GC content of 40-60%; the primers are listed in Table 1.

qRT-PCR analysis
The qRT-PCR was performed on a LightCycler 480 Real Time PCR system (Roche, Basel, Switzerland) with a final volume of 20 μL per reaction, as previously described [30]. Each 20-μL reaction mixture contained 5 μL of diluted cDNA, 10 μL of SYBR Premix Ex Taq II (Takara), 0.8 μL of each primer, and 3.4 μL of sterile distilled water. The thermal cycling program was 95˚C for 30 s, followed by 40 cycles of 95˚C for 5 s, 60˚C for 30 s, and 72˚C for 15 s. Meltcurve analyses were performed using a program with constant heating from 70˚C to 95˚C, followed by 95˚C for 15 s. Each reaction was carried out in three technical replicates.

Stability analysis of reference genes
First, three analytical tools, geNorm (version 3.5) [11], NormFinder (version 0.953) [31], and BestKeeper (version 1) [32], were used to calculate the expression stabilities of the 15 reference genes in all samples from the different treatments. geNorm determines the most stable reference genes from a set of tested genes in a given sample set under a certain experimental condition and calculates the gene expression normalization factor via geometric averaging of all the reference genes [11]. geNorm evaluates gene expression stability by the M value, and ranks the tested genes according to their expression stability by stepwise exclusion of the gene with the highest M value. The lower the M value, the more stable the gene's expression. The optimal number of reference genes was determined by pairwise variation (V) with a threshold of 0.15 [11]. If the value of Vn/n+1 was below 0.15, additional genes were not required for normalization. NormFinder is an algorithm based on a mathematical model of gene expression and uses a solid statistical framework to estimate the expression variation of candidate normalization genes, as well as the variation in sample subgroups, to identify the optimal normalization gene among a set of candidates [31]. A stability value, which is a direct measure of the estimated expression variation for each gene, is provided by NormFinder; thus, the systematic errors introduced during the process of normalization can be evaluated. In geNorm and NormFinder, the quantities obtained from Ct (cycle threshold) values via the delta-Ct method were used for data input, whereas Bestkeeper used raw Ct values. BestKeeper uses the standard deviation (SD) and coefficient of variance (CV) of the Ct values to calculate the stability of a candidate reference gene; the highest SD and CV values equate to the least stable reference genes [32]. Finally, RefFinder (http://fulxie.0fees.us/?type=reference), which is a comprehensive webbased tool that integrates geNorm, NormFinder, BestKeeper, and the comparative delta Ct method, was applied to determine the most suitable reference genes for the overall final ranking [33].
The transcript levels of SgPOD and SgPAL under all treatments, SgLEA under drought, cold, salt, high pH and ABA treatments, SgLOX under wounding and JA treatments, SgHSP90 under heat treatment, and SgPR1 under SA treatment were quantified using the most and least stable reference genes, respectively. The primers were as follows: [34] was used to calculate the relative expression levels. Each biological sample was analyzed using three technical replicates.

Results
Performance of primers and expression profiles of candidate reference genes Fifteen fragments of candidate reference genes (Table 1), obtained from previous S. grandis transcriptome data [8], that exhibited stable expression after grazing (within a 2-fold change in expression level) were selected. Agarose gel electrophoresis and melting curve analysis showed that a specific fragment of the expected size and a single peak were observed (S1 Fig), respectively, after reverse transcriptase-PCR amplification using cDNA templates with each primer pair. The PCR amplification efficiencies of each primer pair across all samples ranged from 99.3% to 107.3% (Table 1).
The expression levels of the candidate reference genes were determined by qRT-PCR across all samples (including normal and different stress conditions) (Fig 1), and a wide range of Ct values, from 15.25 to 30.03, were found, with the values for most of the genes lying between 20 and 26. SamDC exhibited the highest expression abundance, with a mean Ct of 16.91, while SKIP5 had the lowest expression level, with a mean Ct of 26.08 across all samples.
SamDC (3.40 cycles) and UCH (4.33 cycles) showed the smallest expression variations across all the samples. In contrast, SKIP5 (6.34 cycles) and SKIP11 (6.21 cycles) exhibited the largest relative expression variations among all the samples.

geNorm analysis of reference genes
The geNorm results are shown in Fig 2. For the drought treatment, UBC15 and UBC17 were the most stably expressed genes with the lowest M value of 0.0844. For ABA, salt, and JA treatments, the most stably expressed genes were EF1beta and UBC15, with M values of 0.096, 0.106, and 0.117, respectively. For heat treatment, HERC2 and ACT had the lowest M value of 0.137. Under cold treatment, UBC17 and SKIP5 had the lowest M value, and for high pH treatment, UBC17 and CUL4 were the most stable genes. After SA and wounding treatments, the M values of the combinations EF1beta and HERC2, and UBC2 and UBC15 were the lowest, respectively. Under heat and cold treatments, the least stable genes were SKIP5 and SKIP1, respectively. In the other seven treatments, SKIP11 was the most unstable gene. Across all samples, UBC17 and CUL4 ranked as the most stable genes, while SKIP11 was the least stable gene.
The optimal number of reference genes for normalization was calculated by pairwise variation (V). Under various conditions, the Vn/n+1 values were lower than 0.15 (Fig 3), which suggested that adding a further reference gene would not improve the normalization quality and a combination of the two most stable reference genes was sufficient for normalization.

NormFinder analysis of reference gene expression stability
Based on the results calculated by NormFinder (Table 2), the three highest ranked genes under ABA treatment were CUL4, UBC15, and UBC17; under cold treatment were SKIP5, HERC2, and UBC17; under heat treatment were ACT, SKIP1, and GAPDH; under JA treatment were UBC17, EF1beta, and UBC15; under salt treatment were UBC15, EF1beta, and HERC2; under high pH treatment were UBC15, EF1beta, and SKIP5; under SA treatment were HERC2, UBC17, and EF1beta; under wounding treatment were UBC2, UBC15, and UBC17; and under drought treatment were EF1beta, UBC15, and GAPDH. According to all the samples in the study, UBC15, HERC2, and EF1beta were the most stable genes. In addition, the least stable genes ranked by NormFinder and geNorm were consistent: for cold and heat stresses they were SKIP1 and SKIP5, and for the other treatments it was SKIP11.
Bestkeeper analysis of reference gene expression stability As shown in Table 3, Bestkeeper ranked UBC17 as the most stably expressed gene following cold, wounding, and ABA treatments, while UCH had the lowest SD value under drought and salt stresses. GAPDH, SamDC, UBC2, and SKIP1 were ranked as the most stably expressed genes under heat, high pH, JA, and SA treatments, respectively. When all the samples were analyzed together, UCH was rated as the most stable gene with the lowest SD value. RefFinder ranking of the most stable genes RefFinder, which combines the geNorm, NormFinder, and BestKeeper outputs, was applied to generate a comprehensive final ranking of the candidate references. As shown in Table 4, EF1beta (drought and JA), ACT (heat), SKIP5 (cold), UBC15 (salt, high pH, and ABA), UBC2 (wounding), and HERC2 (SA) were ranked as the top stable genes under the different treatments. Across all samples, EF1beta, HERC2, and UCH were suggested as the three best reference genes for normalizing the expression of target genes.

Validation of candidate reference genes
To validate the utility of the candidate reference genes selected using RefFinder, the expression patterns of six target genes identified from RNA-Seq data, including SgPOD, SgPAL, SgLEA, SgLOX, SgHSP90, and SgPR1, in response to different treatments were examined. When normalization was performed using the two most stable genes under each treatment, EF1beta and UBC15 (drought), SKIP5 and UBC17 (cold), UBC15 and EF1beta (salt), ACT and GAPDH (heat), UBC15 and HERC2 (high pH), UBC2 and UBC15 (wounding), EF1beta and UBC17 (JA), UBC15 and CUL4 (ABA), and HERC2 and UBC17 (SA), the expression of SgPOD was affected by different treatments (Fig 4), and peak points were observed at 3 h under drought, at 6 h under cold, at 3 h under salt, at 1 h under high pH, at 48 h under JA, and at 3 h under SA, whereas minimums were found at 24 h under heat and at 48 h under ABA and wounding treatments. It is worth noting that the transcripts of SgLEA were intensively induced by drought, cold, salt, high pH, and ABA treatments (Fig 5A-5E). Wounding and JA treatments caused transcript accumulation of SgLOX (Fig 5G and 5H), and heat stress increased the expression level of SgHSP90 (Fig 5F). An increased transcript level of SgPR1 was observed at 3 h after SA treatment compared with the control, and a peak point was observed at 12 h following a decrease in expression at 6 h ( Fig 5I). Additionally, upregulated expression of SgPAL was observed after JA, high pH, and SA treatments, whereas SgPAL expression was downregulated by the other treatments (S2 Fig). Similar expression patterns were obtained using the most stable single genes and combinations of two genes for normalization under the above treatments, while different expression patterns were observed when using the least stable reference genes compared with the most stable genes for target gene normalization. Discussion qRT-PCR combines reverse transcription (RT) with the quantitative real-time polymerase chain reaction (qPCR) [35]. The accuracy of qRT-PCR is influenced strongly by the stability of Reference Genes Selection in Stipa grandis the reference gene chosen for normalization of gene expression data [22]; therefore, it is crucial to select stably expressed genes as internal references [13]. Base on previous transcriptome data [8], we selected 15 genes that exhibited relatively low levels of transcript variation after grazing stress and evaluated their expression stability under different stresses in S. grandis. Four statistical algorithms, geNorm, NormFinder, BestKeeper, and RefFinder, were used to evaluate the stability of these genes. All the genes displayed different expression patterns following the various treatments. Moreover, high stability levels were observed using geNorm analysis; all the M values were below the accepted threshold of 1.5 (Fig  2), which indicated that all the genes could be used as internal control genes for normalization. However, according to different algorithms and analytical procedures, several differences were observed in the ranking of the most stable genes, which was similar to the results for other species [19][20][21]. For wounding treatment in S. grandis, UBC2 and UBC15 were selected by both geNorm and NormFinder as the two most stable reference genes (Fig 2 and Table 2), but UBC17 and CUL4 were ranked top by Bestkeeper (Table 3); UBC2, UBC15, and UBC17 were determined as the most stable by RefFinder's comprehensive analysis (Table 4). Across all treatments under study, UBC17, CUL4, and HERC2 for geNorm (Fig 2), UBC15, HERC2, and EF1beta for NormFinder (Table 2), and UCH, SamDC, and SNF2 for Bestkeeper were considered the most suitable reference genes (Table 3). Overall, RefFinder recommended using EF1beta, HERC2, and UCH (Table 4).
Traditional reference genes, such as ACT, GAPDH, and EF1α, which are involved in primary metabolism or other cellular processes [19], as well as new candidate reference genes derived from public databases [12,17,22,25,26], have been validated to allow precise interpretation of the qRT-PCR results for specific experimental conditions. However, fluctuating expression of reference genes has been observed under certain experimental conditions. For example, EF1α was identified as the best reference gene under salinity and drought stresses in Glycine max [36]. Nevertheless, EF1α exhibited poor performance in leaves and roots during abiotic stress treatments [15]. In the present study, we also found that the most suitable reference genes varied under different environmental stimuli. The comprehensive combined analysis by RefFinder showed that the top three ranked genes included the traditional genes EF1beta, ACT, UBC15, UBC2, GAPDH, and UBC17, as well as the novel genes UCH, SKIP5, HERC2, and CUL4 (Table 4).
Previously, EF1beta was identified as a stable reference gene in the developing seeds of the Tung Tree [37]. EF1beta has also been determined to be a good reference in soybean under various stresses [36]. Similarly, we found that EF1beta was in the three top ranked reference genes for drought, salt, high pH, JA, and SA treatments, as well as across all samples (Table 4). ACT, a traditional reference gene, is used commonly as an internal control to normalize target gene expression results from qRT-PCR in various species under different conditions. For instance, ACT was determined to be stably expressed during different floral developmental stages and different temperature treatments in Sweet Osmanthus (Osmanthus fragrans Lour.) [14]. In our study, ACT was identified as the best reference gene by the geNorm and RefFinder programs under heat stress (Fig 2 and Table 4).
Ubiquitination is related to post-translational modification of proteins, and three enzymes are involved in this cellular process: E1 (ubiquitin-activating enzymes), E2 (ubiquitin-conjugating enzymes), and E3 (ubiquitin ligases) [38]. Accumulated data indicate that ubiquitinconjugating enzymes (UBC) exhibit stable expression in many experimental conditions; thus, they have been used widely as reference genes in different species [12,14,39]. Three UBC genes (UBC15, UBC2, and UBC17) were analyzed in this study. Notably, these UBC genes performed well under different treatments (excluding heat stress), especially when analyzed by geNorm and RefFinder (Fig 2 and Table 4). UBC15 and UBC17 exhibited better expression stability than UBC2. HERC2, a HECT domain E3 (ubiquitin ligases) encoding gene [40], also showed stable transcript levels under most experimental conditions. Additionally, the de-ubiquitinating protease UCH (ubiquitin carboxyl-terminal hydrolase) was the top ranked gene using Bestkeeper analysis under drought and salt treatments, as well as in the overall assessment (Table 3). Taken together, at least one of the ubiquitin-related genes was ranked in the top three stable genes (Fig 2, Tables 3 and 4). These results indicated the UBC genes could be used more often as an internal control for normalizing gene expression under various stress treatments. SKIP, a spliceosome component and transcriptional activator, is conserved in eukaryotes and has diverse functions in plant development associated with responses to environmental changes [41]. Among the three SKIP encoding genes examined in this study, SKIP5 was ranked as the most stable gene under cold stress using geNorm (Fig 2), NormFinder (Table 2), and RefFinder (Table 4), and SKIP1 was recommended as the best gene under JA treatment by Bestkeeper analysis (Table 3). Conversely, SKIP11 was ranked the least stable reference gene by the different programs (Fig 2 and Tables 2-4), except for heat and cold stresses by geNorm and NormFinder, heat, cold, ABA, and JA treatments by Bestkeeper, and heat, cold, and JA treatments by RefFinder. These results indicate that SKIP11 is the least suitable gene for use as an internal control in S. grandis under various environmental stimuli.
The effectiveness of the most stable reference genes selected using RefFinder was evaluated by normalization of six target genes in samples under different stress conditions. Peroxidases [POD/POX, EC:1.11.1.7] are heme-containing glycoproteins that oxidize a wide range of hydrogen donors, accompanied by the accumulation of H 2 O 2 [42]. Many peroxidase isoenzymes are found in higher plants and they play roles in plant defense, growth, and development [42]. Guaiacol peroxidase (GPX, EC 1.11.1.7) [43], which is generally considered a stress-related enzyme, is involved in plant defense against various stresses, and the activity of GPX is thought to be induced by diverse environmental stimuli [44], such as drought, salinity, and chilling [43,45,46]. In the present work, the expression of a peroxidase encoding gene was detected under the different treatments (Fig 4), and the changes in its transcript levels suggested that peroxidase plays an important role in S. grandis defense against stressful conditions. The late embryogenesis abundant protein encoding gene LEA has been widely reported to participate in responses to drought, cold, and salt stresses, and ABA treatment, and plays an important role in plant tolerance to abiotic stress [47][48][49][50][51]. The significantly increased expression of SgLEA under drought, cold, salt, high pH and ABA treatments (Fig 5A-5E) implies that SgLEA participates in S. grandis resistance responses against stresses. The expression level of SgLOX was elevated by wounding and JA treatments (Fig 5G and 5H), and SgHSP90 expression was enhanced by heat treatment (Fig 5F); the induction of these genes was in accordance with previously reported expression patterns of LOX2 [52][53][54] and HSP90 [55][56][57]. Additionally, a pathogenesis-related protein 1 encoding gene, SgPR1, which is widely accepted as an SA marker gene [58] and is induced by SA accumulation [59], was selected to assess the utility of the identified reference genes under SA treatment. As a critical defense signaling molecule, SA accumulation is associated with enhanced resistance to pathogen infection [60]. The changes in the transcript level of SgPR1 after SA treatment indicated that SgPR1 could be used as a molecular marker for disease resistance responses in S. grandis (Fig 5I). However, when the target gene normalization was performed with the least stable genes, the expression patterns varied in comparison with most stable reference genes. Taken together, these results suggest that appropriate selection of reference genes is necessary for the normalization of target gene expression to ensure the accuracy of qRT-PCR data.

Conclusions
This study comprehensively evaluated the utility of 15 potential reference genes for quantifying transcript levels in S. grandis exposed to diverse experimental conditions. The stable reference genes identified in this work will help to improve the accuracy of gene expression analysis and will aid the investigation of stress-response genes and the molecular mechanism of stress tolerance in this species.