Reference gene selection for qRT-PCR assays in Stellera chamaejasme subjected to abiotic stresses and hormone treatments based on transcriptome datasets

Background Stellera chamaejasme Linn, an important poisonous plant of the China grassland, is toxic to humans and livestock. The rapid expansion of S. chamaejasme has greatly damaged the grassland ecology and, consequently, seriously endangered the development of animal husbandry. To draft efficient prevention and control measures, it has become more urgent to carry out research on its adaptive and expansion mechanisms in different unfavorable habitats at the genetic level. Quantitative real-time polymerase chain reaction (qRT-PCR) is a widely used technique for studying gene expression at the transcript level; however, qRT-PCR requires reference genes (RGs) as endogenous controls for data normalization and only through appropriate RG selection and qRT-PCR can we guarantee the reliability and robustness of expression studies and RNA-seq data analysis. Unfortunately, little research on the selection of RGs for gene expression data normalization in S. chamaejasme has been reported. Method In this study, 10 candidate RGs namely, 18S, 60S, CYP, GAPCP1, GAPDH2, EF1B, MDH, SAND, TUA1, and TUA6, were singled out from the transcriptome database of S. chamaejasme, and their expression stability under three abiotic stresses (drought, cold, and salt) and three hormone treatments (abscisic acid, ABA; gibberellin, GA; ethephon, ETH) were estimated with the programs geNorm, NormFinder, and BestKeeper. Result Our results showed that GAPCP1 and EF1B were the best combination for the three abiotic stresses, whereas TUA6 and SAND, TUA1 and CYP, GAPDH2 and 60S were the best choices for ABA, GA, and ETH treatment, respectively. Moreover, GAPCP1 and 60S were assessed to be the best combination for all samples, and 18S was the least stable RG for use as an internal control in all of the experimental subsets. The expression patterns of two target genes (P5CS2 and GI) further verified that the RGs that we selected were suitable for gene expression normalization. Discussion This work is the first attempt to comprehensively estimate the stability of RGs in S. chamaejasme. Our results provide suitable RGs for high-precision normalization in qRT-PCR analysis, thereby making it more convenient to analyze gene expression under these experimental conditions.


INTRODUCTION
Stellera chamaejasme Linn (Thymelaeaceae), a perennial herb and dominant plant of grassland desertification, is native to the northern and southwestern regions in China (Tseng, 1999). The whole plant is toxic, and its main toxic component is chamaejasmin, which can poison and kill cattle, sheep, and other livestock . The rapid spread of S. chamaejasme speeds up the process of grassland desertification and also poisons a large number of livestock in pasturing areas, causing great damage and loss to the local grassland ecology and livestock husbandry . Thus, it is of fundamental importance to elucidate the mechanisms of the rapid spread and stress adaptation of S. chamaejasme. However, limited genome sequence information is available, which greatly hinders the study of stress functional genes, ultimately resulting in a slow advancement of prevention and control measures. For the above reasons, our group established local transcriptome data for S. chamaejasme seedlings at five different time points (300 mM NaCl treatment for 0, 3, 12, 24, and 72 h; three biological replicates) using the Illumina HiSeq 4000 sequencing platform. After transcriptome sequencing and data analysis, fragments per kilobase of exons per million fragments mapped (FPKM) converted from RSEM (RNA-Seq by expectation maximization) were used to estimate unigene expression, which in some cases led to a few false-positive results.
Quantitative real-time polymerase chain reaction (qRT-PCR) is one of the most widely applied technologies to detect the expression levels of selected genes in many different samples  because of its relatively accurate quantification, simplicity, specificity, high sensitivity, and high throughput capacity (Qi et al., 2016;Wang et al., 2016a). In the relative quantitative method of qRT-PCR data processing, the choice of internal genes is particularly important, and small changes in reference gene (RG) stability will significantly influence the accuracy of the relative expression of target genes . Generally speaking, an ideal RG should be an endogenous gene that does change in any of the tested tissues or under any of the experimental conditions (Derveaux, Vandesompele & Hellemans, 2010;Li et al., 2016a;Wang et al., 2016b). In cells, some endogenous housekeeper genes with consistent relative expression are often used as RGs (Taylor et al., 2016).
As of now, there is no available internal control gene for qRT-PCR data normalization in S. chamaejasme, so we were unable to verify transcriptome sequencing results, analyze the expression patterns of salt or stress-related genes, or further clarify its spread mechanism. To solve this problem, in our study, we selected 10 candidate RGs based on the local salt S. chamaejasme transcriptome database and then determined their expression profiles in five different stages under various abiotic stresses (drought, salt, and cold) and with three hormone treatments (abscisic acid, ABA; gibberellin, GA; ethephon, ETH) by qRT-PCR and further evaluated their expression stabilities using three popular software packages: geNorm (Vandesompele et al., 2002), NormFinder (Andersen, Jensen & Orntoft, 2004), and BestKeeper (Pfaffl et al., 2004). The 10 candidate genes were 18S, 60S, CYP, EF1B, GAPCP1, GAPDH2, MDH, SAND, TUA1, and TUA6. Two target genes, Delta 1-pyrroline-5-carboxylate synthetase 2 (P5CS2), which encodes a crucial enzyme in the proline synthesis pathway under stress conditions by activating glutamate 5-kinase and glutamate-5-semialdehyde dehydrogenase (Strizhov et al., 1997), and GIGANTEA (GI), a circadian regulated gene whose protein product has not only been shown to regulate photoperiodic flowering and various developmental processes but has also been implicated in mediating cold stress and salinity stress responses (Cao, Ye & Jiang, 2005;Penfield & Hall, 2009;Li et al., 2017), were used to verify the selected RGs.

Plant materials and stress treatments
Stellera chamaejasme seeds were collected from Qilian, Qinghai province. After peeling, the seeds were treated with 98% H 2 SO 4 for 9-11 min and were then rinsed for 30 min with running water and planted in individual pots (14.5 Â 14.5 Â 6.5 cm) filled with nutrition soil, vermiculite and perlite (6:1:1). Germinated seeds were grown seven weeks and were then transferred to nurseries potted with double-layered filter paper for three days of adaptation cultivation. All of the nursery pots were placed in an artificial climate chamber at a temperature of 25 ± 2 C during the day and 15 ± 2 C at night, with a relative humidity of 50-55% and illumination intensity of 300 mmol m -2 s -1 (14/10 h, day/night). Three pots of seven-week-old seedlings (three biological replicates) with a consistent growth status for each group were chosen and treated with abiotic stresses and hormone treatments.
For drought and salt treatments, 20% PEG-6000 (w/v, Sangon, Shanghai, China) (Zhuang et al., 2015) and 300 mM NaCl (Sangon, Shanghai, China)  were applied to irrigate the seedlings, respectively. For cold stress, the seedlings in the nursery pots were shifted to another artificial climate chamber at 4 C. For hormone treatments, the leaves were sprayed with 0.1 mM ABA (Reddy et al., 2016;Wan et al., 2017), 0.1 mM GA (Li et al., 2016b), or 1.5 mM ETH (Wu et al., 2016). Seedlings were irrigated or sprayed every 12 h during the course of the experiment. Complete seedlings were carefully collected at 0, 3, 12, 24, and 48 h after treatments; immediately frozen in liquid nitrogen; and stored at -80 C refrigerator until total RNA isolation.
Total RNA isolation and first strand cDNA synthesis Five random individual plants, approximately 100 mg of seedlings in each sample, were used for total RNA isolation with a TRNzol reagent kit (TIANGEN, Beijing, China). The concentration and 260/280 and 260/230 ratios of the RNA samples were detected with a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), and the integrity of all of the RNA samples was verified by 1.0% (w/v) agarose gel electrophoresis (AGE). Subsequently, for reverse transcription PCR (RT-PCR) and qRT-PCR, a total of 3.0 mg of RNA was DNase I (Ambion, Waltham, MA, USA) treated and purified and then used to synthesize first strand cDNA by reverse transcription (Roche, Basel, Switzerland) in a 20 ml reaction system. Finally, cDNA diluted 50-fold with ddH 2 O, was used as the template for PCR amplification.
Candidate RG selection and primer design Ten candidate RGs from the local S. chamaejasme transcriptome database were selected by using local NCBI-blast (version 2.4.0+). The sequences of these genes were used to design the qRT-PCR primers using Primer 5.0, Oligo 7.60, and Beacon Designer 8.20 software with the following criteria: melting temperature TM of 50-65 C, primer lengths of 17-25 bp, GC contents of 45-55%, and product lengths of 90-300 bp. The specificity of all of the selected primer pairs was observed via RT-PCR using the cDNA of control groups at 0 h as the template, and each gene fragmentation was underpinned by 2.0% (w/v) AGE and sequenced to ensure its reliability.

RT-PCR and qRT-PCR analysis
To confirm the specificity of each primer that we designed, we performed RT-PCR in a 25 ml system using the Bio-Rad C1000 PCR system (Bio-Rad, Hercules, CA, USA). The reaction system was as follows: 2.5 ml of Ex Taq buffer, 2 ml of dNTPs, 0.125 ml of TaKaRa Ex Taq (TaKaRa, Beijing, China), 60 ng of cDNA template, 0.2 mM reverse primer, 0.2 mM forward primer, and ddH 2 O to 25 ml. The RT-PCR reaction parameters were: 95 C for 3 min, 40 cycles at 95 C for 30 s, 58 C for 30 s, 72 C for 20 s, and 72 C for 5 min. The amplification products were evaluated by 2.0% (w/v) AGE. To further confirm that the amplicon corresponded to the target sequence, PCR products contained in the agarose gel were extracted using a TIANgel Midi Purification Kit (TIANGEN, Beijing, China) and then sequenced using the dideoxy chain-termination method by Sangon Biotech (Shanghai, China) Co., Ltd.
qRT-PCR reactions were carried out with the Fast Start Universal SYBR GreenMaster (Roche, Basel, Switzerland) on a Bio-Rad CFX96 Real-Time PCR system (Bio-Rad, Hercules, CA, USA) in accordance with the manufacturer's instructions. Reactions were conducted at 95 C for 3 min as an initial denaturation, followed by 40 cycles at 95 C for 10 s, 58 C for 10 s, and 72 C for 20 s. The melting curves, ranging from 58 C to 95 C, were determined to check the specificity of the amplicons. In the negative control group, qRT-PCR was performed using water instead of cDNA as the template. Three technical replicates were analyzed for each biological sample, and the final Ct values for each set of samples were the average of three biological replicates. A total of 45 cDNA samples from five time points in the control groups were used to determine the mean amplification efficiency (E) of each primer pair with the LinRegPCR program (Ruijter et al., 2009;Zhuang et al., 2015;Vavrinova, Behuliak & Zicha, 2016;Wu et al., 2016).

Data analysis of gene expression stability
Three different types of statistical tools: geNorm (version 3.5), NormFinder (version 0.953), and BestKeeper (version 1.0), were applied to rank the expression stability of the RGs across all of the experimental sets. For geNorm and NormFinder, the raw Ct values calculated by the CFX equipment (Tm) software were converted into the relative quantities using the formula 2 -ÁCt (ÁCt = each corresponding Ct valuelowest Ct value) for gene expression profiling. For BestKeeper, the raw Ct values and amplification efficiencies estimated by the LinRegPCR program were used to calculate the coefficient of variation (CV) and standard deviation (SD). The RG with the lowest CV ± SD value was identified as the most stable gene, and the RG with SD value greater than 1.0 was judged to be unstable and should be avoided for gene expression normalization (Guenin et al., 2009). geNorm software was also used to determine the proper RG numbers with pair wise variation (V n /V n+1 , n refers to the RGs number) between two sequential normalization factors.

Validation of RGs
To test the accuracy of the results, the geometric mean from the sort results of geNorm, NormFinder, and BestKeeper in each subset were used to calculate the comprehensive ranking of the candidate genes. The smaller the comprehensive ranking results, the better the gene expression stability. Then, the combination of the top two best RGs, best ranked RG and worst ranked RG were used to standardize the expression of two target genes, i.e., P5CS2 and GI, under different experimental conditions. Furthermore, the expression levels of P5CS2 and GI under salt stress calculated by the combination of the top two best RGs were also compared with the FPKM values in the S. chamaejasme transcriptome database.

Selection of candidate RGs and target genes
After comparing the reported RGs in other species with the local transcriptome database of S. chamaejasme using the local Blast program, 10 RGs and two target genes were chosen to perform the gene normalization studies. The results showed that the E value of each blast gene indicated high homology. The untranslated region of these full-length unigene sequences were used to design the specific primers for RT-PCR and qRT-PCR. The unigene ID, NCBI accession number, gene symbol, gene name, homolog locus of 10 candidate RGs and two target genes, and E value compared with those of the homologous genes are listed in Table 1.

Verification of the primer specificity and qRT-PCR amplification efficiency
The specificity of each primer was tested by 2.0% AGE, sequencing and melting curves analysis, which provided the expected amplicon length (Fig. S1) and single peak melting curves (Fig. S2). The primer sequences, amplicon size, product Tm, amplification efficiencies, and other relevant information are given in Table 2. The amplification product length of PCR varied from 94 to 267 bp. The Tm for all PCR products spanned from 76.0 C for MDH to 83.5 C for GAPCP1. The E values of these genes were between 1.824 (MDH) and 1.930 (GAPDH2), and the linear correlation coefficients (R 2 ) varied from 0.994 (SAND) to 0.998 (CYP). In conclusion, we had every reason to believe that all of these specificity and efficiency estimates of the amplification were reliable for further analysis.

Expression profiles of candidate RGs
Boxplot analysis of the Ct values of different RGs in all of the experimental samples was performed using origin 2017 software (Fig. 1). The results demonstrated that the mean Ct values of the 10 candidate RGs presented a relatively wide field, from 19.26 to 30.76. 60S showed the least expression variation, while 18S exhibited the highest variation, with the Ct values ranging from 15.58 to 22.59. Since the Ct values are negatively related to the gene expression levels, the smaller the Ct value, the higher the gene expression level. As Fig. 1 shows, 18S was the highest-expressed RG for its lowest mean Ct value (15.58), and GAPCP1 had the lowest expression level on account of its maximum mean Ct value (32.58).
Analysis of gene expression stability geNorm analysis. geNorm calculates the gene expression stability measure M value as the average pairwise variation V for the RG and other tested RGs (Vandesompele et al., 2002). The smaller the M value, the more stable the gene, and vice versa. In our study, the M values of the 10 candidate RGs of S. chamaejasme calculated by geNorm software were below 1.5 in all of the experimental settings (Fig. 2), suggesting that these genes should be considered relatively stable. As described in Figs (Fig. 2D), GA (Fig. 2E), and ETH (Fig. 2F) treatment groups, SAND and TUA6, TUA1, and CYP, GAPDH2 and 60S were considered to be the most stable genes with the lowest M values of 0.21, 0.22, and 0.19, respectively. In addition, for all of the sample sets (Fig. 2G), GAPCP1 and CYP were suggested to be two most stable RGs. On the contrary, 18S was the least stable gene in all of the sets except for ETH treatment, in which TUA6 was the least stable gene. NormFinder analysis. NormFinder provides a stability value for each gene by analyzing expression data obtained through qRT-PCR, which is a direct measurement for estimating expression variation when the gene is used for normalization . The orders based on the stability values calculated by NormFinder (Table 3) were similar to those determined by geNorm. The stability ranking results under cold stress and in the GA treatment subsets were completely consistent with the results determined through geNorm; meanwhile, TUA6 and 18S were the two least stable genes for ETH treatment and the rest of the treatments. For the cold stress group, GAPCP1 and EF1B were the two most stable RGs (also ranked first by geNorm). For the salt stress group, GAPCP1 and TUA1 were the two most stable RGs, which was different from the geNorm results. For all samples, ABA-treated and ETH-treated subsets, NormFinder suggested that GAPCP1 and 60S, TUA1 and SAND, GAPDH2 and 60S were the most stable RGs, respectively, which were not exactly the same as the geNorm analysis results.
BestKeeper analysis. BestKeeper evaluates the RG expression stability by calculating the CV and SD of the average Ct values. A lower CV value indicates more stable RG expression (Guenin et al., 2009). As shown in Table 4, under drought stress and for all of the sample subsets, TUA1 had the lowest CV ± SD values of 0.52 ± 0.16 and 0.53 ± 0.16 and was considered to be the most stable RG. Under the cold stress condition and salt stress and ABA treatment subsets, EF1B, which had the lowest CV ± SD values of 1.16 ± 0.31, 1.35 ± 0.36, and 1.04 ± 0.27, respectively, was identified as the best RG.   In the GA treatment subset, TUA6 had the lowest CV ± SD value of 0.82 ± 0.22 and was the most stable RG. In the ETH treatment subset, BestKeeper suggested that GAPDH2 was the most stable RG with the lowest CV ± SD value of 0.68 ± 0.18. Additionally, only a few genes had a SD value greater than 1.0, indicating that most of the candidate RGs were relatively stable. Except for the ETH treatment subset, the most unstable RG among all of the experimental settings was 18S, which was the same as the results of geNorm and NormFinder.

Determination of the optimal number of RGs
At the suggestion of the geNorm Service tool, the critical value V n /V n+1 to determine the optimal RG number for qRT-PCR normalization is 0.15, below which the inclusion of an additional RG is not required (Vandesompele et al., 2002). As Fig. 3 shows, the V2/3 values of all of the experimental groups were less than 0.15, which indicated that a two RG combination would be sufficient to use for normalization.
Comprehensive stability analysis of RGs Table 5 and Fig. 4 summarize and rank the determination results obtained from the geNorm, NormFinder, and BestKeeper programs. Based on the analysis, GAPCP1 and EF1B were the most stable RGs under three abiotic stresses; thus, TUA6 and SAND, TUA1 and CYP, GAPDH2 and 60S were the best RG combinations under the ABA, GA, and ETH treatments, respectively. Still, 18S was the most unstable RG under all of the experimental conditions.

Reference gene validation
As shown in Figs. 5 and 6, when the best RG combinations were used for performing normalization, the expression levels of P5CS2 and GI were affected by different treatments. A sustained increase in expression level of P5CS2 was observed after drought stress, and a peak point was observed at 48 h (Fig. 5A). A tendency of first increase, after downward, and then upward in the transcript level of P5CS2 appeared after cold and GA treatments (Figs. 5B and 5E). Additionally, upregulated expression of P5CS2 was observed after salt and ABA treatments, and reached the maximum value at 12 and 24 h following a decrease (Figs. 5C and 5D). Whereas, P5CS2 expression was first treatments also appeared prominent changes, which were 4.66-fold, 29.22-fold, 2.10-fold, 6.45-fold, and 2.45-fold higher than those of the control group, while there was no significant difference under GA treatment (Fig. 6E). Compared with the best RG combinations for normalization of P5CS2 and GI, similar expression patterns were obtained when the most stable single genes, GAPCP1 (drought and cold), EF1B (salt), TUA6 (ABA), TUA1 (GA), and GAPDH2 (ETH), were used for normalization under the above treatments. However, different expression patterns were generated and the expression levels of P5CS2 and GI were overestimated when the least stable gene, 18S, was selected as the RG for normalization.
In particular, as shown in Fig. 5C, under salt treatment, when the RG combination (GAPCP1 and EF1B) was selected for normalization, gene expression of P5CS2 gradually increased from 0 h, reached the maximum at 24 h, and then began to slightly decline at 48 h. In the same way, the expression levels of GI increased at first, then decreased and maintained a lower level until 48 h (Fig. 6C). The expression trends of P5CS2 and GI over the first 24 h were generally consistent with those of RNA-seq (Fig. S3), which further validated the accuracy and reliability of our experimental results. DISCUSSION qRT-PCR is currently viewed as a powerful technique that can be used to quantify target gene expression. The accuracy of qRT-PCR directly depends on the stability of the internal genes used. The use of inappropriate RGs for normalizing qRT-PCR data will lead to deviations in the results (Shivhare & Lata, 2016). In this study, three programs, geNorm, NormFinder, and BestKeeper, were used to select optimum RGs for six different experimental conditions. The 10 potential RGs exhibited differential stability in response to different stresses. Taking ABA treatment as an example, in the experimental subset, geNorm software ranked SAND as the best RG, and NormFinder regarded TUA1 as the most stable RG. However, BestKeeper identified EF1B as the best RG according to its lowest CV value. This means that the three types of software generated different results, and a solution was not found. Our study carried out a comprehensive analysis and provided ultimate stability ordering results by ranking the geometric means of the three software analysis results, which is a common strategy for evaluating the expression stability of RGs reported in previous scientific papers.
EF1B catalyzes the exchange of guanosine diphosphate bound to the G-protein, elongation factor 1-alpha (EF1A), for guanosine triphosphate, an important step in the elongation cycle of protein biosynthesis. It has been considered to be one of the most stable RGs during drought and salt stresses (Wan et al., 2017) and other stress conditions (Ma et al., 2013). In our study, EF1B was ranked as one of the two best RGs under drought and salt stresses according to the comprehensive ranking results, which was the same as in Stipa grandis (Wan et al., 2017). In addition, EF1B performed a better expression stability under three abiotic stresses than those in the three hormone treatments. The results showed that there are no universal RGs that are stably expressed in all biological materials and/or under all trial conditions. The expression stability of two homologous RGs, TUA1 and TUA6, were estimated in our study. According to the results, the stability ranking of TUA1 was always better than that of TUA6 under all conditions except the ABA treatment, under which TUA6 exhibited better expression stability. Nevertheless, it is notable that the homologous RGs showed different rank orders in each subset, and in most cases, TUA1 showed better expression stability than TUA6. Cordoba et al. (2011) found that TUA1 was one of the most suitable RGs under NaCl stress and 2,4-dichlorophenoxyacetic acid treatment in Hedysarum coronarium. However, Gimeno et al. (2014) suggested that TUA6 should be discard for normalization under drought stress, salt stress, cold and heat shock treatment, and flooding treatment in switchgrass. Although these reports indicated the expression stability of TUA1 and TUA6 in different species, respectively, but at present we have not found any reports of the simultaneous use of TUA1 and TUA6 in stability analysis in other species. Therefore, this does not mean that the expression stability of TUA1 in other species must be better than that of TUA6.
18S is a frequently used HKG and is widely used for normalization in qRT-PCR analysis. Wang et al. (2017) reported that 18S rRNA was one of the most stably expressed gene under diverse heavy metals stresses in tea plants; Huang et al. (2017) also found that 18S rRNA was the most stable gene under UV irradiation and hormonal stimuli in Baphicacanthus cusia. However, our analysis results suggested that 18S was the most unstable RG in all of the experiment groups because of its excessively high expression level. In comparison with the best RG combination and the most stable RG, when 18S was selected as a RG to validate the expression of the two target genes P5CS2 and GI, their expression patterns were significantly overestimated, which was consistent with the findings in Oxytropis ochrocephala (Zhuang et al., 2015) and rice (Bevitori et al., 2014).
Two target genes, P5CS2 and GI, were used to verify the stability of the selected RGs for gene expression normalization. Strizhov et al. (1997) stated that expression of Arabidopsis thaliana P5CS (AtP5CS) is root and leaf specific and can be regulated by salinity, drought, and ABA. The same experimental results were reproduced in our experiments. We also found that P5CS could be efficiently expressed during the later period of cold stress, which may be a supplement to previous findings. The induction mechanism remains to be further studied. Cao, Ye & Jiang (2005) reported that GI is induced by cold stress, but not by salt, mannitol, and ABA. By contrast, Park,  and Kim et al. (2013) claimed that GI, as a negative regulator, participated in the regulation of salt stress in Arabidopsis by interacting with salt overly sensitive 2. Moreover, Riboni et al. (2016) revealed that ABA affects flowering through two independent regulatory mechanisms: activation of GI and constant (CO) functions upstream of the florigen genes and down-regulation of the suppressor of overexpression of CO1 (SOC1) signaling. Our findings indicated that the gene expression of GI not only changed under salt stress and cold stress but also underwent a significant change under drought, ABA, and ETH treatments. We have reason to believe that these mechanisms will be revealed with future in-depth experiments.
There is no doubt that it is necessary to select suitable RGs and/or RG combinations for gene normalization studies to obtain more accurate and reliable results. Combined with all of the validation results above, we observed that, in most cases, P5CS2 and GI showed similar response patterns when normalized by the most stable RGs combinations, but some differences still emerged. Unfortunately, we could not tell which choice was better for normalization. However, to eliminate the small variations caused by technical protocols in qRT-PCR, two or more RGs are often required to correct for non-specific experimental variation (Thellin et al., 1999;Bustin et al., 2009). In this study, two RG combinations, whose V2/3 values were less than 0.15 across all of the experimental subsets and thus GAPCP1 and EF1B for drought stress, cold stress and salt stress, TUA6 and SAND for ABA treatment, TUA1 and CYP for GA treatment, and GAPDH2 and 60S for ETH treatment, were suggested for the accurate normalization of target gene expression.

CONCLUSION
This study represents the first attempt to comprehensively analyze the stability of RGs for use as internal controls in qRT-PCR analysis of target gene expression in S. chamaejasme under three abiotic stresses and three hormone treatments by combining results from three different methods. The results indicated that the stability of an identical gene was not exactly the same under different treatments, and the stability ranking of the RGs calculated by three parameters was not identical under the same treatment. As a result, it makes sense to carry out a comprehensive analysis of the results of the three procedures. Moreover, it may be a better choice to select a combination of two or more RGs as an effective internal control to further improve the accuracy and reliability of gene expression normalization under different stresses. In conclusion, this study provides a guideline to select a valid RG combination that can ensure more accurate qRT-PCR-based gene expression quantification and basic data to facilitate future molecular studies on gene expression in S. chamaejasme and other Thymelaeaceae species (Che et al., 2016).