Reference genes selection for transcript normalization in kenaf (Hibiscus cannabinus L.) under salinity and drought stress

Plant materials and treatments

Kenaf (H. cannabinus L.) variety Fuhong 992 was used for all experiments. To achieve disease-free materials, seeds were rinsed under running water for 10 min and sterilized with 5% sodium hypochlorite for 10 min, washed three times with sterile water, and then germinated on filter paper that had been saturated with water in complete darkness at 28 °C. After 3 days, seedlings were grown in the greenhouse in 25% Hoagland solution under a 16/8-h light/dark photocycle at 28/26 °C (day/night). The most consistent seedlings at the 3–5 leaf stage were used for excess salinity and drought treatments or for the harvesting of different treated leaf samples. For the excess salinity treatment, seedlings were subjected to 200 mM NaCl and harvested at 2, 4, 6, 8, 12, and 24 h. For the drought treatment, seedlings were treated with 20% (w/v) PEG6000, and samples were collected at several time points (2, 4, 6, 8, 12, and 24 h). Untreated seedlings were harvested as a control at 0 h. Three biological replicates were obtained from each group. All treated-leaf samples were employed in experiments that evaluated the stabilities of candidate reference genes under excess salinity and drought stresses.

Total RNA extraction and cDNA synthesis

All samples were snap-frozen in liquid nitrogen and stored at −80 °C before RNA extraction. The OMEGA isolation kit (R6827-01, USA) was used for RNA extraction, following the elimination of genomic DNA by RNase-free DNase I (TaKaRa, Japan); RNA integrity was then assessed by 2% agarose gel electrophoresis, and RNA sample quality was determined using a NanoDrop 2000 spectrophotometer (NanoDrop, Thermo Scientific). Finally, RNA samples with an A260/A280 ratio of 1.9–2.1, and an A260/A230 ratio greater than 2.0 were used for further analyses. Subsequently, first-strand cDNA was synthesized in a 20-µl reaction using the PrimeScript^® RT reagent kit (TaKaRa, Japan) following the manufacturer’s protocol. The quality and integrity of the cDNA were checked by NanoDrop 2000 spectrophotometry and agarose gel electrophoresis, respectively, and stored at −20 °C until use.

Primer design, verification of PCR products, and qRT-PCR

The sequences of 10 candidate reference genes (18S rRNA, ACT4, EF1 α, MZA, PP2A, PTB, RAN, TUB α, UBC and UBQ) from the transcriptomic data (L Zhang, 2015, unpublished data) of kenaf were subjected to a BLAST search using A. thaliana sequences in GenBank. The probe sequences from A. thaliana were used to search the transcriptome of kenaf and the kenaf expressed sequence tags database (http://www.ncbi.nlm.nih.gov/nucest/?term=kenaf). Target sequences were identified by querying homologous kenaf sequences together with A. thaliana complete CDS and probe sequences. According to these potential reference gene sequences, primers were designed using Primer 3 (http://bioinfo.ut.ee/primer3-0.4.0/primer3/) based on the following criteria: GC content 45–80%, melting temperatures 58–62 °C, primer lengths 18–24 bp, and amplicon lengths 100–250 bp (see Table 1 for detailed primer information). To check the specificity of the amplicon, all primer pairs were initially tested via standard RT-PCR using the Premix Ex Taq™(TaKaRa, Japan); the amplification products of each gene were verified by 2% agarose gel electrophoresis. Real-time amplification reactions were performed with the Applied Biosystems 7500 Real-Time PCR System using SYBR^® Premix Ex Taq. Reactions were prepared in a 20-µl volume containing 2 µl cDNA template, 0.4 µl each amplification primer, 0.4 µl ROX Reference Dye II, 10 µl 2 × SYBR Premix Ex Taq, and 6.8 µl dH2O. Amplifications were performed with an initial step of 95 °C for 30 s, followed by 40 cycles of denaturation at 95 °C for 5 s and primer annealing at 60 °C for 34 s. The melting curve ranged from 60 °C to 95 °C, and the temperature was increased in increments of 0.2 °C per 10 s for all PCR products. ABI Prism Dissociation Curve Analysis Software was used to confirm the occurrence of specific amplification peaks. All reactions were carried out in triplicate with template-free negative controls being performed in parallel.

Statistical analyses

To select a suitable reference gene, geNorm, NormFinder, and BestKeeper were used to analyze the stability of each candidate reference gene. All analyses using these packages were performed according to the manufacturer’s instructions. The raw Ct values from qRT-PCR were transformed into the required data input format. The maximum expression level (the minimum Ct value) of each gene was used as a control and was set to a value of 1. Then relative expression levels were calculated from Ct values using the formula: 2^−ΔCt, where ΔCt = each sample Ct value—minimum Ct value. The geNorm and NormFinder algorithms further calculated the expression stability value (M) for each gene with the obtained data, whereas the BestKeeper analyses were based on untransformed Ct values. Standard curves were generated using Microsoft Excel 2003 by plotting cycles at threshold fluorescence (Ct) against the logarithmic values of standard RNA amounts. Quantities of standard RNA were prepared by diluting cDNA (1, 10⁻¹, 10⁻², 10⁻³, 10⁻⁴, and 10⁻⁵; each gene in triplicate). Only Ct values less than 40 were used to calculate correlation coefficients (R² values) and amplification efficiencies (E) from the slope generated in Microsoft Excel 2003, according to the equation E = [10^−(1/slope) − 1] × 100%. All PCR assays showed efficiency values between 97% and 118% (Table 1). All other multiple comparisons were performed with the statistical analysis software SPSS 22.0 (SPSS Inc., USA). To validate the reference gene(s) selected in the current study, the relative expression level of WRKY28 was normalized using the 2^−ΔΔCt method after collecting the mean Ct value of each biological replicate from the samples treated under conditions of salinity or drought for 0, 4, 6, 8, 12, and 24 h. Finally, the relative increases in expression level of WRKY28 were used to calculate the differences in the normalization of each reference gene.

Selection of candidate reference genes, primer specificity, and amplification efficiency

A total of 10 candidate reference genes, including four commonly used housekeeping genes, 18S rRNA, ACT4, TUB α, and UBQ, and six novel candidate reference genes, EF1 α, MZA, PP2A, PTB, RAN, and UBC, were identified in this study (Table 1). The four novel candidates were validated in A. thaliana, O. sativa, or G. hirsutum for expression stabilities under different abiotic stresses. RAN and UBC have been evaluated as the optimal reference gene in Cucumis melo (Kong et al., 2014) and Platycladus orientalis (Chang et al., 2012), respectively. Additionally, the specificity of the designed primers was verified by a single band with the expected size after agarose gel electrophoresis (Fig. S1). Specificity was further confirmed by a single peak in the melting curve analysis, which was done prior to performing qRT-PCR (Fig. S2). A standard curve was generated using a 10-fold dilution of cDNA in the qRT-PCR assay to determine the amplification efficiency for each primer pair. Both E and R² were calculated using the slope of the standard curve. Results indicated that the average amplification efficiency values of all primers ranged from 97.13% to 118.60% (Table 1).

As shown in Fig. 1 and Table 1, for all tested samples, the mean Ct values of 10 candidate reference genes had a wide range (17.89–34.28), and the standard deviation (SD) varied from 0.65 to 1.94, and the coefficient of variation (CV) ranged from 2.15% to 8.64%. Comparing to Ct values of all the candidates, 18S rRNA had a highest expression level (mean Ct ± SD = 17.89 ± 1.55), following ACT4 (mean Ct ± SD = 22.17 ± 1.24) and EF1 α (mean Ct ± SD = 22.46 ± 1.54), whereas MZA (mean Ct ± SD = 34.28 ± 0.74) accumulated less than all the others. Additionally, PP2A, PTB, RAN, TUB α, UBQ and UBC showed a narrow range of mean Ct values, indicating that these genes were stably expressed. However, it is insufficient for evaluating gene expression stability just by the comparison of raw Ct values. To obtain accurate gene expression data, other approach should be combined with to validate a set of appropriate reference genes under a given conditions.

Stability of reference genes analysis by geNorm

The geNorm-based analysis was conducted to determine which reference gene(s) would be optimal in each tested sample set. As shown in Fig. 2, the candidate genes were ranked accordingly to their M values. Since a lower M values indicates a greater stability of the reference gene, TUB α and 18S rRNA were ranked as the most stable, with M values of 0.28 and 0.65, respectively, whereas UBC was the least stably expressed under conditions of drought and excess salinity (Fig. 2). The results were consistent with the pattern observed in Tables 2 and 3. When the combination of drought and salinity was considered, the same results (TUB α and 18S rRNA) were acquired for normalization (M = 0.51) (Table 4 and Fig. 2C). The geNorm algorithm can also be used to identify the optimal number of reference genes by calculating the pairwise variation (V) between normalization factors (NFn). It is proposed that an additional reference gene makes no sense to the normalization when Vn/n+1 is less than 0.15. In this study, the data showed that a V2/3 of 0.13 was less than 0.15, which indicated that the combination of TUB α and 18S rRNA was sufficient for the normalization of gene expression under drought stress (Fig. 2D). For salinity stress in kenaf, a V3/4 of 0.15 was less than a V2/3 of 0.19, which suggested that three reference genes, TUB α, 18S rRNA, and RAN were the best options for accurate normalization under salinity stress (Fig. 2D). For all drought and excess salinity samples, the same effects were observed. These results revealed that TUB α together with 18S rRNA (V2/3 = 0.14) could provide a reliable reference for the normalization of gene expression.

Stability of reference genes analysis by NormFinder

The NormFinder program analyzes candidate reference genes according to inter- and intragroup variations in expression. As in the geNorm method, the gene with the lowest M value has the most stable expression. As shown in Tables 2–4, the NormFinder analysis also identified that 18S rRNA and TUB α were the most stably expressed genes with values of 0.12 and 0.13 under drought stress, respectively, with slight differences in the ranking order (Table 3). For the salinity samples, PP2A and ACT4 were the top two reference genes followed by 18S rRNA, identified by NormFinder (Table 2). For the drought and excess salinity samples, ACT4 and 18S rRNA were the top two genes calculated by NormFinder (Table 4). Compared with the results of geNorm and NormFinder analyses, 18S rRNA was the best reference gene under salinity stress. Overall, the results obtained from NormFinder analysis were consistent with those from geNorm analysis.

Stability of reference genes analysis by BestKeeper

The BestKeeper program utilizes the unconverted Ct values to accomplish parametric tests on normally distributed expression levels. It evaluates the geometric mean of Ct values, determines the coefficient of variance (CV), and calculates a Pearson’s coefficient of correlation (r) of each gene. The standard deviation (SD) of the average Ct values are performed to generate a weighted index of most suitable normalizing genes across selected biological samples and exclude genes that are not stably expressed. The most stable genes are identified as those that exhibit the lowest CV and SD (CV ± SD), in relation to r and p values. Genes with SD greater than 1 are considered unacceptable and should be excluded. In this study, 18S rRNA had CV ± SD values of 1.96 ± 0.67, showing remarkably stable expression (r = 0.911; p = 0.004). TUB α displayed the lowest CV ± SD values of 1.50 ± 0.46, but with a lower r value (r = − 0.186), and higher p-value (p = 0.692), it is not possible to conclude on stable gene. UBQ, PTB, UBC, EF1 α and ACT4 had an SD > 1.0, thus excluding them as suitable reference genes for NaCl-treated samples (Table 2). These results are consistent with those acquired from the geNorm method. For the PEG-treated samples, MZA displayed the lowest CV ± SD = 0.72 ± 0.25, but with r = 0.246, p = 0.593, it does not act as a suitable reference gene. TUB α (CV ± SD = 1.54 ± 0.46, r = 0.837, and p = 0.019) and 18S rRNA (CV ± SD = 4.47 ± 0.85, r = 0.940, and p = 0.002) showed the potential for being normalization factors. Whereas UBQ and RAN had an SD value greater than 1.0, indicating that the two reference genes were not suitable for normalization (Table 3). When the two datasets from conditions of excess salinity and drought were analyzed together, although MZA and TUB α had a lower CV ± SD, a lower r value and higher p value was observed. It is notable that 18S rRNA with SD >1.0 was identified as the least stable by BestKeeper compared to the other candidate reference genes (Tables 3 and 4).

Reference gene validation

To further validate the reliabilities of the selected reference genes in the current study, the relative expression level of two transcription factors, HcWRKY28 and HcWRKY32, were investigated in two kenaf samples under conditions of drought and excess salinity, using one or two of the most stable reference genes, as well as the least stable gene for normalization, which had been determined by geNorm, NormFinder, and BestKeeper as described above (Tables 2–4 and Fig. 2). In A. thaliana, AtWRKY28 (At4g18170) plays a crucial role in protecting plants against pathogen infections and oxidative stress. In our previous study, we found that kenaf WRKY28, a homolog of AtWRKY28, was differentially expressed under abiotic stress conditions according to the gene differential expression profiling analysis based on RNA-seq. We also found that WRKY32 was actively upregulated in kenaf in response to drought and salinity stresses (X Niu, 2015, unpublished data). Therefore, the two genes were selected to further validate the usefulness of the selected reference genes by qRT-PCR. In the current study, we selected the two reference genes with lowest M values, 18S rRNA and TUB α and their combination (18S rRNA + TUB α), and the gene (UBC) with highest M value, as the normalization factors to determine the expression patterns. As shown in Fig. 3, the expression abundance of WRKY28 increased significantly after 6 h of salinity treatment, peaked at 8 h, and thereafter decreased (P < 0.01) (Fig. 3A), while WRKY32 peaked at 12 h (P <0.01) (Fig. 3C) in the NaCl-treated samples. For PEG-treated samples, qRT-PCR analysis showed that the expression level of WRKY28 decreased from 0 h to 6 h, increased at 8 h, peaked at 12 h, and then decreased (P < 0.01) (Fig. 3B). In contrast, the expression level of WRKY32 increased at 12 h, and peaked at 24 h with approximately 7.5-fold accumulation (P < 0.01) (Fig. 3D). Salinity and drought stresses significantly increased the expression of WRKY28 by 4.5- and 3.5-fold at 8 h and 12 h, WRKY32 by 6.5- and 7.5-fold at 12 h and 24 h, compared with that at 0 h, respectively (P < 0.01). These results are in accordance with the patterns of WRKY28 and WRKY32 in previously analysis. Moreover, the expression patterns of WRKY28 and WRKY32 showed similar trends across different experimental sets to those of 18S rRNA and TUB α as reference genes either singly or in combination (Fig. 3), although some slight differences exist. However, when the expression of WRKY28 and WRKY32 was normalized using the reference gene with the highest M value (UBC) as an internal control, the expression pattern of WRKY28 and WRKY32 showed a significantly different changes, and the expression level was 5-fold lower than that of the genes with the lowest M values (P < 0.01) (Fig. 3). These results indicate that 18S rRNA and TUB α are suitable reference genes for gene expression normalization under conditions of drought and excess salinity in kenaf.