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

Plant materials and stress treatments

Stellera chamaejasme seeds were collected from Qilian, Qinghai province. After peeling, the seeds were treated with 98% H₂SO₄ 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 μmol m⁻²s⁻¹ (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) (Wang et al., 2015) 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 μg 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 μl reaction system. Finally, cDNA diluted 50-fold with ddH₂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 ™ 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 μl system using the Bio-Rad C1000 PCR system (Bio-Rad, Hercules, CA, USA). The reaction system was as follows: 2.5 μl of Ex Taq buffer, 2 μl of dNTPs, 0.125 μl of TaKaRa Ex Taq (TaKaRa, Beijing, China), 60 ng of cDNA template, 0.2 μM reverse primer, 0.2 μM forward primer, and ddH₂O to 25 μl. 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 value − lowest 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²) 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. 2A–2C, GAPCP1 and EF1B under drought stress, GAPCP1 and 60S under cold stress, and EF1B and 60S under salt stress were the most stable RGs with the lowest M values of 0.07, 0.05, and 0.19, respectively. At the same time, in the ABA (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 (Dheda et al., 2005). 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 downregulated at 3 h after ETH treatment and then began to continuous increase, reaching the maximum at 48 h (Fig. 5F). The maximal expression levels of GI under drought (Fig. 6A), cold (Fig. 6B), salt (Fig. 6C), ABA (Fig. 6D), and ETH (Fig. 6F) 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.