Three kenaf varieties, Fuhong992, Fuhong952, and GV42, were used for stress treatments. Seeds were germinated on filter paper saturated with water, and seedlings were grown in nutritional soil in an artificial climate chamber (28°C, 16/8 h light/darkness; light intensity 300 μmol·m⁻²·s⁻¹; relative humidity 50 ± 10%). For each variety, three 5-leaf-stage seedlings from the same experimental treatments were collected for 3 replicates. For the salinity and drought treatments, the details refer to the method of Niu et al. (2016). Cold treatment was performed by exposing seedlings to 4°C in light and harvested at the same time points. For the hormone treatments, leaves were sprayed with 0.3 mM abscisic acid (ABA treatment), 5.0 mM ethylene (ET treatment), 0.1 mM jasmonic acid (JA treatment), and 0.2 mM salicylic acid (SA treatment) and then harvested at the same time points above. All samples were harvested from three biological replicates in different pots, giving a total of 132 samples comprised of 57 abiotic and 75 hormone treatment samples for each variety. All samples were frozen in liquid nitrogen and stored at −80°C.

The OMEGA isolation kit (R6827-01, Omega Bio-tech, USA) was used for total RNA isolation. The concentration and purity of each RNA sample was estimated by a NanoDrop ND 2000 spectrophotometer (NanoDrop, Thermo Scientific). RNAs with an OD260/OD280 value between 1.9 and 2.1, and OD260/OD230 > 2.0 were used for further analyses. The first-strand cDNA was synthesized using the PrimeScript® RT reagent kit (TaKaRa, Japan) and stored at –20°C for further use.

Nine candidate reference genes, 18S rRNA, ACT7, UBC, PTB, PP2A, TUB8, MZA, RAN, and EF1α, were screened from the Arabidopsis database. The potential homologs were identified by querying from the kenaf expressed sequence tags (ESTs) database and the kenaf transcriptome database (SRP060459) (Zhang et al., 2015). Based on these potential reference gene sequences, primers were designed for detection and cloning of reference genes using Primer v6. The detailed information of primers, such as GC content, melting temperatures, primer lengths, and amplicon lengths, were presented in Table 1. Cloning information was provided in Table S1. The specificity of PCR products were confirmed by agarose gel electrophoresis and sequencing.

qRT-PCR was performed in 96-well plates by using a Applied Biosystems 7500 Real-Time PCR System (Applied Biosystems, USA). cDNA was amplified by using SYBR® Premix Ex Taq (Tli RNaseH Plus, TaKaRa, Japan), and a final volume of 20 μl reaction contains 10 μl 2 × SYBR Premix Ex Taq, 2 μl cDNA template, 0.4 μl each amplification primer, 0.4 μl ROX Reference Dye II, and 6.8 μl dH₂O. Thermal cycling was performed with an initial step of 95°C for 30 s, 40 cycles at 95°C for 5 s, and 60°C for 34 s. The final melting curve was obtained from 60° to 95°C to verify primer specificity. All assays were carried out in three technical and biological replicates with template-free negative controls being performed in parallel. The amplification efficiency (E) and correlation coefficient (R²) were determined by the standard curves, based on the 10-fold cDNA diluted series.

Three available statistical algorithms, geNorm, NormFinder, and BestKeeper were used for evaluating the expression stability of each candidate gene. All analyses were conducted following the operation instructions. The raw Ct values of all reference genes were transformed into the required data format. The relative quantities were calculated according to the formula: 2^−ΔCt, where ΔCt = corresponding Ct—minimum Ct. A comprehensive ranking of reference genes was generated as suggested by Niu et al. (2015a,b). Statistical analysis of gene expression data were performed using the software SPSS 22.0 (SPSS Inc., USA).

On the basis of the homology analysis with Arabidopsis genes, 9 candidate reference genes (Hc18S rRNA, HcACT7, HcUBC, HcPTB, HcPP2A, HcTUB8, HcMZA, HcRAN, and HcEF1α) were identified from the H. cannabinus transcriptome (SRP060459). The full-length cDNA sequences of these genes were cloned from H. cannabinus variety Fuhong992, and the sequence information was deposited in GenBank, and the accession numbers were showed in Table 1. Subsequently, specific primers were designed, and the specificity of the primers was confirmed on the basis of the amplification efficiency and specificity: (a) all primers were amplified with a single band after 1.5% agarose gel electrophoresis (Figure S1); (b) a single peak in the melting curve analysis; and (c) the primers amplification efficiency (E) ranged from 97.13 to 118.6%, and the correlation coefficients (R²) of the standard curve varied from 0.9949 to 0.9997 (Table 1 and Figure S2).

Transcript abundances of 9 candidate reference genes were assessed by qRT-PCR in 396 samples, and the Ct values were examined under three groups including abiotic stresses, hormonal stimuli, and all samples. Lower Ct values represent higher expression abundance, and higher Ct values represent lower expression abundance. As shown in Figure 1 and Table 1, the mean Ct values ranged from 17.98 to 34.04, with SD varied from 1.69 to 2.46 and CV from 4.9 to 12.5% among these 9 candidate genes. These data means that Hc18S rRNA expressed abundantly (Ct < 18); HcACT7, HcUBC, and HcEF1α were moderately expressed genes (21 < Ct < 24); HcPTB, HcPP2A, and HcTUB8 showed slightly lower expression levels (28 < Ct < 31); and HcRAN and HcMZA showed the lowest expression levels with Ct values as high as 32 and 34 cycles, respectively. Additionally, HcMZA showed the least variable gene (CV = 4.9%), followed by HcPP2A (6.2%), and Hc18S rRNA showed the most variability (12.5%) among the 9 candidate genes. Based on CV values, the stability ranking of all reference genes was as follows: HcMZA > HcPP2A > HcPTB > HcRAN > HcTUB8 > HcUBC > HcEF1α > HcACT7 > Hc18S rRNA (Table 1). In Brief, these data indicate that transcription levels of reference genes are unstable under different experimental conditions.

The geNorm statistical algorithm was used for evaluating expression stability of the 9 candidate references. This algorithm determines the gene expression stability (M) based on the average pairwise variation of all tested genes (Vandesompele et al., 2002; Liu et al., 2016). The lower M means the gene expresses more stable. For kenaf variety Fuhong992, the candidate genes in each experimental group showed much lower M values than the threshold value of 1.5. Therefore, HcPP2A and HcPTB (M = 0.08) were ranked as the most stable, HcUBC (M = 0.19) was the least gene in the abiotic stress group. Under the hormonal stimuli conditions, HcACT7 and Hc18S rRNA (M = 0.10), followed by HcTUB8 (M = 0.11), were identified as the most stable genes, HcUBC (M = 0.26) showed the worst stability. When all samples was considered, Hc18S rRNA and HcACT7 (M = 0.09) featured the most stable genes whereas HcEF1α (M = 0.21) and HcUBC (M = 0.25) were the least genes (Table 2). For kenaf varieties Fuhong952, HcPTB and HcEF1α (M = 0.03) were listed as the best genes, HcUBC and HcMZA were the worst in the abiotic stress group. Under the hormonal stimuli conditions, HcACT7 and HcTUB8 (M = 0.02), followed by HcPP2A (M = 0.04), were identified as the most stable genes, while HcUBC (M = 0.13) gene showed the worst stability. As to total samples, HcPP2A and Hc18S rRNA (M = 0.07) featured the most stable genes whereas HcTUB8 (M = 0.21) was the least (Table S2). For kenaf varieties GV42, HcPP2A and HcACT7 expressed stably, HcUBC and HcEF1α expressed unstably under abiotic stresses, hormonal stimuli, and total samples. (Table S3).

The optimal number of reference genes was also determined by geNorm algorithm based on the pairwise variation (V_n) between normalization factors (NF_n). The threshold value (V_n/V_n+1 = 0.15) indicates that the number of reference genes less than or equal to the value of n is sufficient to use as reference gene. As depicted in Figure 2, the pairwise variation value V₂/V₃ of all experimental samples was less than 0.15, suggesting that two reference genes should be sufficient for normalization under these conditions.

The NormFinder approach determines the stability of candidate reference genes based on inter- and intra-group variations in expression across different sample groups (Andersen et al., 2004). Similar to the geNorm analysis, the lower M value means the gene expressed more stable. We obtained the similar results through the NormFinder analysis, which predicted HcPP2A, HcACT7, and HcPTB to be the three stable genes for abiotic stress group, hormonal stimuli group, and the total samples group in Fuhong992 and GV42 (Table 2; Table S3). For kenaf varieties Fuhong952, HcPP2A, Hc18S rRNA, HcPTB, HcACT7 were the most stable reference genes under different conditions (Table S2). Overall, HcUBC and HcMZA were ranked as the worst reference gene, which was consistently found in the three subsets.

The BestKeeper program was applied for evaluating the candidate references stability based on CV and SD of Ct values (Pfaffl et al., 2004). The lower CV ± SD value, the more stable expressed. The analysis results of BestKeeper were also listed in Table 2. For kenaf variety Fuhong992, HcMZA (CV ± SD = 0.34 ± 0.12) and HcPP2A (CV ± SD = 1.50 ± 0.42) showed remarkably stable expression in the abiotic stress group. HcRAN (1.08 ± 0.37) and HcPP2A (5.10 ± 0.92) were identified as the two best reference genes in the different hormonal stimuli group. When all samples were analyzed together, the two reference genes, HcRAN (3.35 ± 1.19) and HcMZA (3.54 ± 1.21), showed the highest expression stabilities in the total group (Table 2). For kenaf variety Fuhong952, HcRAN (0.71 ± 0.24) and HcEF1α (1.45 ± 0.29) showed stable expression in the abiotic stress group. HcRAN (0.32 ± 0.10) and HcPP2A (2.36 ± 0.65) were identified as the two best reference genes in hormonal stimuli group. HcPTB (2.35 ± 0.64) and HcPP2A (2.37 ± 0.66) showed the highest expression stabilities in the total group (Table S2). For kenaf variety GV42, HcRAN and HcTUB8 were the two best reference genes in abiotic stress group, hormonal stimuli group, and total group (Table S3). All these results were inconsistent with those obtained from the geNorm and NormFinder analyses.

To develop a consistent result, a comprehensive ranking of reference genes was obtained according to the method described by Niu et al. (2015b). The results were shown in Table 2, Tables S2, S3 for Fuhong992, Fuhong952, and GV42, respectively. Based on geNorm algorithm, the best combination of reference gene in each group was further obtained. As shown in Table 3, HcPP2A and HcPTB were the most stable references under abiotic stress condition; HcPP2A and HcACT7 in the different hormonal stimuli and total samples groups for kenaf variety Fuhong992. For Fuhong952, HcRAN, HcPP2A, and HcPTB were listed as the most stable genes in the abiotic stress group; HcPP2A and HcACT7 in the different hormonal stimuli group; HcPP2A and HcPTB, followed by Hc18S rRNA in the total samples group. For GV42 variety, HcPP2A and HcACT7 showed the optimal performance in the abiotic stress group, hormonal stimuli group, and the total samples group (Table 3).

To confirm the reliability of the selected references, the relative expression level of HcWRKY28 was normalized using the most stable and the least stable reference genes by qRT-PCR analysis in kenaf variety Fuhong992. As shown in Figure 3, the relative transcript abundance of HcWRKY28 under different abiotic stress or hormonal stimuli was found to be impartial when HcPP2A, HcPTB, and HcACT7 or their combinations (HcPP2A and HcPTB or HcACT7) were used as normalization factors. Using references HcPP2A, HcPTB, or combination of HcPP2A and HcPTB, the HcWRKY28 gene was significantly up-regulated under salinity and cold stress at 8 h, and was down-regulated under drought stress at 6, 8, and 12 h (Figure 3A). Correspondingly, for different hormonal stimuli, the HcWRKY28 showed relative up-regulation (JA and SA stimuli) or down-regulation (ABA and ET stimuli) at each time point when HcPP2A, HcACT7 or their combination were used for transcript normalization (Figure 3B). However, when the worst reference gene HcUBC used as an internal control, the expression pattern of HcWRKY28 showed strong fluctuations, and the gene expression change level was 5-fold higher than that by the most stable genes under different abiotic stress and hormonal stimuli (Figure 3).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer YQ declared a shared affiliation, though no other collaboration, with several of the authors XN, XH, HC, AT, JX, JQ to the handling Editor, who ensured that the process nevertheless met the standards of a fair and objective review.