Selection and Validation of Reference Genes for the qRT-PCR Assays of Populus ussuriensis Gene Expression under Abiotic Stresses and Related ABA Treatment

Populus ussuriensis Kom. is one of the most important tree species for forest renewal in the eastern mountainous areas of Northeast China due to its fast growth, high yield, and significant commercial and ecological value. The selection of optimal reference genes for the normalization of qRT-PCR data is essential for the analysis of relative gene expression. In this study, fourteen genes were selected and assessed for their expression stability during abiotic stress (drought, high salinity, and cold stress) and after the treatment with the drought-related hormone ABA. Three algorithms were used, geNorm, NormFinder, and BestKeeper, and a comprehensive ranking of candidate reference genes was produced based on their output. The most appropriate reference genes were UBQ10 and RPL24 for drought and ABA treatment, UBQ10 and TUB3 for cold stress, and UBQ10 and 60S rRNA for high salinity. Overall, UBQ10 was the most stable reference gene for use as an internal control, whereas PP2A was the least stable. The expression of two target genes (P5CS2 and GI) was used to further verify that the selected reference genes were suitable for gene expression normalization. This work comprehensively assesses the stability of reference genes in Populus ussuriensis and identifies suitable reference genes for normalization during qRT-PCR analysis.


Introduction
Gene expression analysis is the basis of modern molecular biology and is widely used in genetic and developmental studies. Gene expression can be analyzed by several RNA quantification methods, including Northern blots, RNase protection assays, semi-quantitative reverse-transcription PCR, and quantitative real-time PCR (qRT-PCR), all of which rely on normalization procedures to quantitatively compare multiple samples [1]. qRT-PCR has become a common technique for gene expression analysis because of its rapidity, high sensitivity, and quantitative accuracy [2,3]. Moreover, qRT-PCR is regarded as the only effective way to measure the expression of genes with low mRNA copy number [4]. However, due to differences in RNA quality and quantity, in the efficiency of DNA synthesis, in PCR amplification, in pipetted volumes, and especially, in the activity of tissues and cells, qRT-PCR results can be severely biased [5]. A set of MIQE (Minimum Information for Publication of Quantitative Real-Time PCR Experiments) guidelines aimed at improving the reliability and repeatability of research using qRT-PCR, has been published, which refers to a suitable reference gene (RG) for normalization as being essential for obtaining reliable results [6]. RGs whose expression

Plant Materials and Stress Treatments
Populus ussuriensis clone Donglin plants were grown in vitro on half-strength Murashige and Skoog (1/2 MS) medium (Phytotech) with 0.6% (w/v) agar and 2% (w/v) sucrose [15]. Plants were cultured in 250 mL plastic containers that held 100 mL of medium and were housed in a growth chamber at 25 • C with a 16 h light/8 h dark photoperiod and a light intensity of 46 µmol m −2 s −1 . For stress treatments, 5-leaf-stage in vitro plants were transferred to 1/2 MS medium supplemented with 7% PEG 6000 (drought), 150 mM NaCl (high salinity), or 100 µM abscisic acid (ABA) medium that had been precooled at 4 • C (cold stress). All media were adjusted to pH 5.8 and sterilized by autoclaving at 121 • C for 20 min. Leaf and root tissues from each treatment were independently collected at six time points (0 h, 6 h, 12 h, 24 h, 48 h, and 72 h). Non-treated plants were used as controls. All samples were immediately frozen in liquid nitrogen and stored at −80 • C for analysis. There were three biological replicates of each experimental condition.

Total RNA Extraction and cDNA Synthesis
RNA extraction was carried out according to the MIQE guidelines to ensure the reliability of the results [6]. Total RNA from each sample (approximately 100 mg of leaf or root tissue) was isolated using the CTAB (cetyltrimethylammonium bromide) method [24] and treated with DNaseI to remove genomic DNA contamination. The concentration and quality of each RNA sample were measured spectrophotometrically using the OD260/OD280 and OD260/OD230 absorption ratios, and all RNA samples were verified by 1% (w/v) agarose gel electrophoresis. cDNA was synthesized using the HiScript II Q RT SuperMix for qPCR (+gDNA wiper) (Vazyme, Nanjing, China) in a 20 µL reaction system according to the manufacturer's instructions. The cDNA was diluted 50-fold with ddH 2 O and used as the template for PCR amplification.

Candidate Reference Gene Selection and PCR Primer Design
Fourteen candidate genes from a local P. ussuriensis transcriptome database were selected for assessment of their stability and suitability as RGs for the normalization of qRT-PCR data under various experimental treatments. Primers were designed based on the cDNA sequences of the candidate reference genes using Primer Premier 6 software (http://www.premierbiosoft.com/) with the following criteria: primer length 17-25 bp (optimum 20), product size 75-300 bp, melting temperature 50-65 • C (optimum 60 • C), and GC content 40%-60% (optimum 50%). The specificity of all selected primer pairs was assessed by PCR using the cDNA of control groups at 0 h as the template, and gene fragments were visualized following 2.0% (w/v) agarose gel electrophoresis to ensure reliability. All primers mentioned above are listed in Table 1.

qRT-PCR
TransStart ® Top Green qPCR SuperMix (TransGen, Beijing, China) was used to perform qRT-PCR in 96-well optical reaction plates (Applied Biosystems, Foster City, CA, USA). Reactions were prepared in 20 µL volumes containing: 10 µL of 2 × TransStart ® Top Green qPCR SuperMix, 7.2 µL of ddH 2 O, 2 µL of 50-fold diluted cDNA, and 0.4 µL of each specific primer, prepared to a final concentration of 10 µM. The qRT-PCR program consisted of an initial step at 94 • C for 30 s, followed by 40 denaturation cycles at 94 • C for 5 s and primer annealing at 60 • C for 30 s. Next, melting curves ranging from 60 • C to 95 • C with an increment of 0.2 • C were recorded in each reaction to check the specificity of the amplicons. Three technical replicates were analyzed for each biological sample. The threshold cycle (Ct) was automatically measured. A total of 48 cDNA samples from six time points in the control groups were used to determine the mean amplification efficiency (E) of each primer pair using the LinRegPCR program. The relative expression levels of two target genes (PuP5CS2 and PuGI) were calculated with the equation 2 -∆∆Ct under abiotic stress and ABA treatment, respectively [25].

Statistical Analysis to Determine Expression Stability with geNorm, NormFinder, and BestKeeper
To analyze the expression stability of 14 candidate RGs, the cycle thresholds (Cts) from different tissues under different stresses and treatments were analyzed using box plots and three different algorithms of the statistical tools geNorm [21], NormFinder [22], and BestKeeper [23]. The geNorm tool calculates the average expression stability value (M value) of RGs based on the average pairwise variation (V) value. With 1.5 as the threshold, a lower M value corresponds to a more stable gene. The V value calculated by geNorm determines the optimal number of RGs required to accurately normalize the data and is based on a threshold value of Vn+1 < 0.15. The NormFinder tool analyzes the stability of candidate RG expression based on the results of intra-and inter-group variation analysis and automatically ranks the genes according to their stability value (SV) [22]. A lower SV corresponds to higher stability. BestKeeper assesses the stability of RGs based on the standard deviation (SD) and coefficient of variation (CV) of the Ct values. An RG with an SD less than 1 is considered to be stably expressed. The CV decreases with the SD, indicating that the RG is more stable [23,26].

Validation of Candidate Reference Genes
To assess the accuracy of the expression stability of candidate RGs, the geometric means of the ranking results from geNorm, NormFinder, and BestKeeper for each stress or treatment were used to calculate a comprehensive ranking of candidate genes. The lower the comprehensive ranking, the better the stability of gene expression. Finally, the two best RGs for each experimental condition were used to normalize the expression of two target genes, PuP5CS2 and PuGI, under that condition.

Selection, Amplification Efficiency, and Ct Value Range of Candidate Reference Genes
We compared RGs used previously in other plant species, such as Arabidopsis thaliana, with a local P. ussuriensis transcriptome database using a local Blast search in Bioedit. Fourteen candidate RGs were selected for use in the gene normalization studies (Table 1). To evaluate the amplification efficiency of primer pairs, PCR amplification specificities were assessed using 2.0% agarose gel electrophoresis ( Figure S1), melting curves ( Figure S2), and direct sequencing. The results showed that we obtained specific target fragments of the expected lengths and sequences. The qRT-PCR amplification efficiency (E) and correlation coefficient (R 2 ) were obtained from the slope of the calibration curves (Table 1). For the 14 candidate RGs and two target genes, the E values ranged from 1.85 to 2.08 and the R 2 values ranged from 0.9918 to 0.9996 (Table 1). The results showed that all 16 pairs of primers met the requirements for use in qRT-PCR experiments.

Expression Profiles of Candidate Reference Genes under Abiotic Stresses and Hormone Treatment
The Ct value is the cycle number at which the reaction curve intersects the threshold line; it indicates the number of cycles required to detect a real signal from the sample. High abundance cDNA samples reach this threshold at a lower Ct value during PCR amplification, indicating a higher gene expression level [27]. The transcript abundances of 14 RGs were determined based on Ct values from the qRT-PCR analysis of 48 samples, including six time points for roots and leaves under four treatment conditions (drought, high salinity, cold stress, and ABA treatment). Mean Ct values ranged from 15.15 (UBQ10) to 29.86 (SAND) (Figure 1). Under all stresses and treatments, UBQ10 had the highest expression and SAND had the lowest expression. The coefficient of variation (CV) is related to the degree of fluctuation in Ct values. The highest and lowest CVs were 7.25% (TUA2) and 0.85% (UBQ10) under drought, 9.4% (TUA2) and 0.94% (F-box) under ABA treatment, 13.30% (UBC) and 1.80% (UBQ10) under cold stress, and 11.07% (RPL25) and 1.18% (SAND) under high salinity ( Figure 1). The average CVs of the candidate RGs were ranked as UBQ10

Expression Stability of Candidate Reference Genes under Abiotic Stress and Hormone Treatment
To identify the most suitable RGs for gene expression analysis under drought, high salinity, cold stress, and ABA treatment in P. ussuriensis, the expression stability of candidate RGs was evaluated by three different algorithms in geNorm, NormFinder, and BestKeeper.

geNorm Analysis
The geNorm tool calculates an RG stability value (M value) based on the average pairwise variation V value [21]. The more stable the gene, the lower the M value, and vice versa. As shown in Figure 2 and Table 2  Sometimes a single RG does not meet the stability requirements for use in normalization, and two or more RGs are therefore needed to reduce error and obtain more accurate results [28]. geNorm can determine the optimal number of candidate RGs based on the calculation of pairwise variation (Vn/Vn+1). The program suggests a threshold of 0.15 as a standard for selecting the best pairwise variation value (V value) for normalization. There is no need to add more candidate RGs if the V value is lower than 0.15. As shown in Figure 3, the value of V5/V6 among 14 candidate RGs under salt stress was 0.125, which is smaller than 0.15, suggesting that the most suitable RG combination contained four genes, namely 60S rRNA, ACT7, TUA2, and RPL24. Similarly, the most appropriate RG combinations under drought, ABA, and cold stress were ACT7 and UBQ10, ACT7 and RPL24, and TUB3 and EF1-α, respectively. Table 2. Expression stability of candidate reference genes calculated by geNorm.

Rank
Drought ABA Cold High Salinity Sometimes a single RG does not meet the stability requirements for use in normalization, and two or more RGs are therefore needed to reduce error and obtain more accurate results [28]. geNorm can determine the optimal number of candidate RGs based on the calculation of pairwise variation (Vn/Vn+1). The program suggests a threshold of 0.15 as a standard for selecting the best pairwise variation value (V value) for normalization. There is no need to add more candidate RGs if the V value is lower than 0.15. As shown in Figure 3, the value of V5/V6 among 14 candidate RGs under salt stress was 0.125, which is smaller than 0.15, suggesting that the most suitable RG combination contained four genes, namely 60S rRNA, ACT7, TUA2, and RPL24. Similarly, the most appropriate RG combinations under drought, ABA, and cold stress were ACT7 and UBQ10, ACT7 and RPL24, and TUB3 and EF1-α, respectively.

NormFinder Analysis
The NormFinder program analyzes the stability of candidate RG expression based on the results of variance analysis and ranks the genes according to their stability [22]. As shown in Table 3, UBQ10 was the most stable candidate RG under drought, similar to the results of geNorm. But unlike the results of geNorm, UBQ10 was also ranked first in cold stress and ABA treatment. However, TUB3 (ranked first by geNorm) in cold stress and RPL24 and ACT7 (ranked first by geNorm) in ABA treatment had stabilities similar to that of UBQ10, a result that differed from the geNorm analysis. For high salinity, 60S rRNA was the most stable RG, which is again different from the result of the geNorm analysis.

BestKeeper Analysis
BestKeeper assesses the stability of RGs by calculating the standard deviation (SD) and coefficient of variation (CV) of the Ct values. An RG with an SD of less than 1.0 is considered to be stably expressed. The CV decreases with the SD, indicating that the RG is more stable [23,26]. As shown in Table 4, UBQ10 was the most stable RG with the lowest SD and CV values under drought, high salinity, cold stress, and ABA treatment. Additionally, the SD and CV values of twelve candidate genes were less than 1.0 under drought stress, and only a few genes had SD and CV values greater than 1.0 under salinity stress and ABA treatment, indicating that the expression of most candidate RGs was relatively stable.

Comprehensive Stability Analysis of the Reference Genes
The final comprehensive RG ranking is shown in Table 5 and Figure 4, based on the results from geNorm, NormFinder, and BestKeeper. UBQ10 and RPL24 were the most stable RGs under drought stress and ABA treatment, UBQ10 plus TUB3 was the best combination under cold stress, and UBQ10 plus 60S rRNA was the best combination under high salinity stress. UBQ10 was the most stable RG under all experimental conditions.

Reference Gene Validation
To verify the accuracy of the stable expression of the RGs, two marker genes (PuP5CS2 and PuGI) were used as positive controls in roots and leaves under drought, high salinity, cold stress, and ABA

Reference Gene Validation
To verify the accuracy of the stable expression of the RGs, two marker genes (PuP5CS2 and PuGI) were used as positive controls in roots and leaves under drought, high salinity, cold stress, and ABA treatment (Figures 5 and 6). P5CS2 is a regulatory enzyme in proline biosynthesis and plays a crucial role in abiotic stress resistance and tolerance [29][30][31][32]. GI encodes a nuclear-localized protein that regulates both freezing tolerance and various development processes such as photoperiod-mediated flowering [33]. We selected two RGs based on the comprehensive ranking results for each stress and used them to normalize the PuP5CS2 and PuGI expression data. The selected RGs were UBQ10 and RPL24 for drought, ABA, UBQ10, and TUB3 for cold stress, and UBQ10 and 60S rRNA for high salinity. When the two optimal RGs were used for normalization, the expression levels of PuP5CS2 and PuGI differed among treatments (Figures 5 and 6). Under drought stress, PuP5CS2 was induced in roots and reached a maximum expression value at 48 h (Figure 4). PuP5CS2 was also upregulated by ABA treatment and reached a maximum value at 24 h (Figure 4). Under cold stress, PuP5CS2 was upregulated in both roots and leaves and was significantly induced at 48 h and 72 h. PuGI was significantly upregulated (>20 fold) in leaves under cold stress, and its expression also responded to drought and high salinity ( Figure 5). Similar expression patterns were observed when we used the two single appropriate genes as RGs. However, the expression patterns of PuP5CS2 and PuGI showed significant fluctuations when the least stable gene, PP2A, was selected as the RG, and their expression patterns were not consistent with those observed when more appropriate RGs were used in specific experimental conditions (Figures 5 and 6).

Discussion
To improve forest productivity and maintain the forest's economic benefits, it is necessary to identify stably-expressed RGs for use in qRT-PCR normalization when studying transcriptional regulatory pathways and the functions of key genes in P. ussuriensis. Abiotic stresses such as drought,

Discussion
To improve forest productivity and maintain the forest's economic benefits, it is necessary to identify stably-expressed RGs for use in qRT-PCR normalization when studying transcriptional regulatory pathways and the functions of key genes in P. ussuriensis. Abiotic stresses such as drought,

Discussion
To improve forest productivity and maintain the forest's economic benefits, it is necessary to identify stably-expressed RGs for use in qRT-PCR normalization when studying transcriptional regulatory pathways and the functions of key genes in P. ussuriensis. Abiotic stresses such as drought, high salinity, and cold temperatures are the main environmental factors that limit forest growth. ABA is a key signaling molecule during abiotic stress response and has particularly important biological functions in the response to drought [34], high salinity [35], and low temperature [36]. Under these stresses, increased levels of endogenous ABA can upregulate many stress-related genes in an ABA-dependent manner [37]. Numerous studies have identified and validated RGs in rice [38], maize [39], soybean [40,41], and carrot [42] under abiotic stress. However, there are currently no reports of RGs that are universally and stably expressed under abiotic stress in P. ussuriensis. The purpose of our study was therefore to select suitable RGs for qRT-PCR in P. ussuriensis under various abiotic stresses.
Fourteen candidate RGs were selected, including traditional housekeeping genes (ACT7, UBQ10, GAPDH, TUA2, TUB3, and EF1-α) and widely-used RGs (60S rRNA, SAND, CYP2, PP2A, F-box, RPL24, RPL25, and UBC). The suitability of these candidate RGs was then evaluated using three algorithms including geNorm, NormFinder, and BestKeeper. Under drought, all algorithms ranked UBQ10 as having the highest stability. However, in other treatments, differences were observed in gene stability values and relative rankings generated by the different algorithms. For example, ACT7 was ranked second by geNorm and NormFinder, but it was ranked in the middle of the 14 candidate RGs by BestKeeper. Under ABA treatment, RPL24 and ACT7 were ranked among the top three most stable genes by geNorm and NormFinder, but they were ranked in intermediate positions by BestKeeper. These differences may be due to the different algorithms used to calculate stability. Previous studies have reported similar differences in stability ranking [43]. In general, the use of more than two algorithms is recommended to screen for stably-expressed RGs [21].
In our study, we calculated a comprehensive ranking of candidate genes based on the geometric mean of the rank results from three algorithms in each experimental condition. UBQ10 ranked as the most stable RG under all stresses and ABA treatment in P. ussuriensis. Our results were consistent with the stable expression of UBQ reported in Platycladus orientalis under NaCl and ABA treatments [44], white clover leaves and stolons under water-limited and well-watered conditions [45], Brachypodium distachyon in different tissues and under various hormone treatments [46] and A. thaliana green siliques and seeds under hormone treatment [47]. However, UBQ showed less-stable expression in salt-stressed Suaeda aralocaspica seeds [48] and in different developmental stages and photoperiodic treatments in soybean [13]. UBQ, which is primarily involved in proteolysis in the ubiquitin-proteasome system, appears to be regulated differently in different plant species and cell types and in response to different experimental conditions [46,49]. Therefore, although UBQ10 was a universally reliable RG in our study, further verification is needed to determine its suitability for use as an RG under other experimental conditions.
According to the comprehensive ranking, RPL24 emerged as a relatively-stable choice for drought and high salinity stresses, and ABA treatment. RPL24, named as 60S ribosomal protein L24, is one of several large ribosomal subunit proteins widely used as an RG in many studies. For example, 60S rRNA was stably expressed among different developmental stages in Panax ginseng [50]. Others include RPL30, which has been used as an RG in roots and shoots of soybean seedlings under cold stress and ABA treatment [40], and RPL2, which is stably expressed under nitrogen, cold, and light stress in tomato [8]. Ribosomal protein genes are also widely used as RGs in animals [51][52][53]. In our study, RPL25, which encodes a different large ribosomal subunit protein, was also evaluated. However, RPL25 was not an optimal RG under any condition in P. ussuriensis. In our study, the performances of RPL24 and 60S rRNA were within an acceptable range in all stresses and ABA treatments. Actin has been used extensively as a traditional HKG for normalization of gene expression data [54]. In our study, geNorm and NormFinder ranked ACT7 among the top four stable genes in roots and leaves under drought, high salinity, and ABA treatment. As noted previously, the expression of ACT7 is relatively stable under ultraviolet treatment in rice, drought, high salinity, and cold stress in common bean, and high salinity and cold stress in peanut [11,55,56]. However, its expression varies in pollen, seeds, and roots in A. thaliana [47]. Therefore, the expression of traditional housekeeping genes is not always stable under all experimental conditions. SAND and TUB3 were two of the top three stable genes under cold stress. Han et al. [57] and Reid et al. [58] reported that SAND was a suitable RG for gene expression analysis in iron-deficient A. thaliana and in developing grape berry tissues. However, SAND has not been reported previously as a stable RG under low-temperature stress in any species. Similar to our results, previous studies have reported that TUB was stably expressed as an RG in various developmental stages of soybean [13], different tissues of poplar [59], and the barks of Populus yunnanensis cuttings [60]. The expression of TUB was relatively stable in leaves from powdery mildew-infected wheat [61] and in P. orientalis seedlings under cold, heat, salinity, PEG, and ABA treatment [44]. TUB was also reported as a suitable RG in all tissues (root, stem, and leaf) under multiple abiotic stress (cold, heat, salinity, and drought), and hormone treatments (ABA, SA, ET, and GA3) in maize [62]. These results suggest that the expression of the same RG varies among organs, developmental stages, and environmental conditions.
We analyzed the expression levels of two target genes (PuP5CS2 and PuGI) to validate the applicability and stability of suitable, high-ranking RGs. Consistent results were obtained by comparing the top two high-ranking RGs in roots and leaves under drought, high salinity, cold stress, and ABA treatment. The expression of PuP5CS2 was significantly induced in roots under drought, high salinity, and ABA treatment. Similar expression patterns of SbP5CS1 and SbP5CS2 under drought and high salinity were reported in sorghum [29]. In A. thaliana, AtP5CS1 responded to ABA and salt stress through controlling proline accumulation [63]. The expression of PuGI was upregulated significantly in leaves during cold stress but not following ABA treatment, which is consistent with the behavior of GI in A. thaliana. In addition, GI was not induced by drought stress or ABA treatment in A. thaliana, but PuGI was upregulated under drought stress and ABA treatment in our study. This result may be because different species and RGs were used in the experiments.
We compared the results of PuP5CS2 and PuGI expression analysis performed using appropriate and inappropriate RGs. When the appropriate RGs were used for normalization, the expression patterns of PuP5CS2 and PuGI were similar, but some differences still emerged. As previous studies have noted, more RGs are often required to correct for non-specific experimental variation and small differences in qRT-PCR conditions [27]. Nevertheless, when the expression data for PuP5CS2 and PuGI were normalized using inappropriate RGs, there were significant biases. These results indicate that the RGs screened in this study were reliable.
In summary, we evaluated and validated stable RGs in P. ussuriensis under abiotic stresses (drought, high salinity, and cold stress) and following treatment with the drought-related hormone ABA. The appropriate RGs for the qRT-PCR normalization of target gene expression data were UBQ10 and RPL24 for drought and ABA treatment, UBQ10 and TUB3 for cold stress, and UBQ10 and 60S rRNA for salinity stress. Suitable RGs differed among various stresses and treatments. However, UBQ10 appeared to be the most reliable RG across all experimental conditions in this study.

Conclusions
We selected fourteen candidate RGs to validate their expression stability using three algorithms, namely geNorm, NormFinder, and BestKeeper during abiotic stress (drought, high salinity, and cold stress) and after the treatment with the drought-related hormone ABA in P. ussuriensis. The results of the three algorithms were analyzed using comprehensive ranking. We identified the most appropriate reference genes, which were UBQ10 and RPL24 for drought and ABA treatment, UBQ10 and TUB3 for cold stress, and UBQ10 and 60S rRNA for high salinity. UBQ10 was the most stable reference gene for use as an internal control, whereas PP2A expression was the least stable. Furthermore, the expression of two target genes (P5CS2 and GI) confirmed the importance of selecting appropriate RGs in gene expression studies. This study provides a guideline for molecular studies of gene expression in P. ussuriensis under various abiotic stresses and provides the foundation for genomic researches in woody species.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/4/476/s1, Figure S1: Agarose gel electrophoresis of PCR products for each of the fourteen RGs and the two target genes; Figure S2: Dissociation curves of the qRT-PCR amplicons for 14 RGs.

Conflicts of Interest:
The authors declare no conflict of interest.