Candidate Reference Genes Selection and Application for RT-qPCR Analysis in Kenaf with Cytoplasmic Male Sterility Background

Cytoplasmic male sterility (CMS) is a maternally inherited trait that results in the production of dysfunctional pollen. Based on reliable reference gene-normalized real-time quantitative PCR (RT-qPCR) data, examining gene expression profile can provide valuable information on the molecular mechanism of kenaf CMS. However, studies have not been conducted regarding selection of reference genes for normalizing RT-qPCR data in the CMS and maintainer lines of kenaf crop. Therefore, we studied 10 candidate reference genes (ACT3, ELF1A, G6PD, PEPKR1, TUB, TUA, CYP, GAPDH, H3, and 18S) to assess their expression stability at three stages of pollen development in CMS line 722A and maintainer line 722B of kenaf. Five computational statistical approaches (GeNorm, NormFinder, ΔCt, BestKeeper, and RefFinder) were used to evaluate the expression stability levels of these genes. According to RefFinder and GeNorm, the combination of TUB, CYP, and PEPKR1 was identified as an internal control for the accurate normalization across all sample set, which was further confirmed by validating the expression of HcPDIL5-2a. Furthermore, the combination of TUB, CYP, and PEPKR1 was used to differentiate the expression pattern of five mitochondria F1F0-ATPase subunit genes (atp1, atp4, atp6, atp8, and atp9) by RT-qPCR during pollen development in CMS line 722A and maintainer line 722B. We found that atp1, atp6, and atp9 exhibited significantly different expression patterns during pollen development in line 722A compared with line 722B. This is the first systematic study of reference genes selection for CMS and will provide useful information for future research on the gene expressions and molecular mechanisms underlying CMS in kenaf.


INTRODUCTION
Cytoplasmic male sterility (CMS) is a maternally inherited trait and is widely used in hybrid seed production in many crops (Schnable and Wise, 1998;Li et al., 2014). The phenotype of a CMS plant is similar to that of normal plant but differs in the lack of functional pollen. In flowering plants, pollen development is a complex process that contains a wide variety of genes (Sun et al., 2013).
These genes are mainly involved in energy metabolism, stamen development (B-class and C-class), and biosynthesis. Abnormal expression of these genes causes loss of pollen function. For example, the expression of DcMADS2 and DcMADS3, which are homologs genes with B-class, were significantly down-regulated in "carpeloid" CMS flower of carrot. Research findings show that specific cytoplasmic (mitochondrial) influence the expression of MADS box genes for B-activity during early differentiation of petals and stamens leading to the "carpeloid" CMS phenotype (Linke et al., 2003). Thus, the analysis of the expression patterns of key genes that are involved in pollen development is useful to explore the molecular regulatory mechanisms underlying plant CMS and pollen development.
Numerous methods have been used to investigate gene expression analysis, such as: northern blotting, semi-quantitative reverse transcription PCR, in situ hybridization, real-time quantitative PCR (RT-qPCR), and RNase protection. The RT-qPCR technology is used worldwide for gene expression analysis due to its high sensitivity, specificity, accuracy, and broad quantification range Coito et al., 2012). However, multiple non-specific variations affect the reliability of gene expression analyses, such as RNA quality, amount of input RNA, primer specificities and the efficiency of RNA reverse transcription (Udvardi et al., 2008). RT-qPCR data must be normalized using a reference gene to avoid these non-specific variations (Radonic et al., 2004;Huggett et al., 2005).
The ideal reference gene must have transcription abundance similar to that of the target gene under variable conditions and must not be co-regulated with the target gene (Radonic et al., 2004). Housekeeping genes are considered to be stably expressed regardless of the organism's developmental stage or environmental conditions, and they have been considered as ideal reference genes. Currently, a number of traditional housekeeping genes have been widely used as reference genes for the normalization of quantitative gene expression, including 18S , ELF1A (Marum et al., 2012;Yue et al., 2016), GAPDH , and CYP (Qi et al., 2010). Nevertheless, studies have shown that no single reference gene can always maintain stability in response to changes in the environmental conditions, as noted with 18S (Wan et al., 2010), ELF1A (Xu et al., 2014), GAPDH (Tong et al., 2009; and CYP (Wei et al., 2013;Yeemin et al., 2016). Consequently, a systematic selection and validation of candidate reference genes in each experimental condition should be performed before selecting a reference gene for the normalization of data in a gene expression analysis.
Kenaf (Hibiscus cannabinus L.) is an annual fiber crop having versatile applications e.g., used in cordage, sacking, paper, building materials, animal bedding, carpet backing . Kenaf spontaneous CMS mutant "UG93" and several CMS lines have been reported (Zhou et al., 2008). The CMS lines have been utilized in hybrid seed production for kenaf crops, but little is known regarding the molecular mechanism of the CMS. To understand the molecular mechanisms underlying of kenaf CMS, our previous work had been performed to compare the transcriptional differences between the CMS and its maintainer lines using RNA-seq (Chen et al., 2014), and many differentially expressed genes were cloned including cox3, HcPDIL5-2a, HcMADS-box, and MYB21 (Chen and Zhou, 2011;Jin et al., 2011;Chen et al., 2012;Qian et al., 2016). Analyzing the expression patterns of genes of interest might aid our understanding of the molecular mechanisms underlying CMS in kenaf crop. However, there is no systematic study conducted for reliable reference gene selection at different pollen development phases in kenaf, which limits further the differential gene studies and transcription analysis.
Therefore, this study assessed the expression stabilities of 10 candidate reference genes (ACT3, ELF1A, G6PD, PEPKR1, TUB, TUA, CYP, GAPDH, H3, and 18S) for the normalization of data at three pollen developmental stages in CMS line 722A and maintainer line 722B, using five statistical algorithms. The expression stability of potential reference genes was validated by normalizing the RT-qPCR data with HcPDIL5-2a. Furthermore, the expression patterns of five mitochondria F 1 F 0 -ATPase subunit genes atp1, atp4, atp6, atp8, and atp9, in lines 722A and 722B of kenaf were examined using the newly identified stable reference genes.

Plant Materials
In this study, kenaf CMS line 722A and its maintainer line 722B were used. Both lines 722A and 722B had similar nuclear genetic background, but with different cytoplasm. Lines 722A and 722B were sown in the experimental field of Guangxi University under natural conditions in the year 2015. Standard cultural practices (irrigation, weeding, insecticide, and fertilization) were carried out according to the crop's demand over the whole growth period. The line 722A plants had floral organs with aberrant anthers (Figures 1A, C). However, line 722B plants had floral organs with normal anthers, and mature pollen grains were produced during pollen development phasse (Figures 1B, D). Anther samples were collected at tetrad (floral bud with length of 2.5-3.5 mm), monokaryotic (3.5-6.0 mm) and dual-core (greater than 6.0 mm) stages of pollen development in lines 722A and 722B, which represented the early, middle and late stages of pollen abortion in line 722A, respectively. Plant samples were quickly frozen in the liquid nitrogen and stored at −80 • C for RNA extraction. Three replications were used and each replicate comprised of 30 individuals.

Semi-Thin Sections and Microscopy
Different lengths of floral buds from the CMS and its maintainer lines were placed in a cold Carnoy's fixative solution (ethanol: acetic acid = 3:1) for 24 h at 4 • C. The floral buds were dehydrated using graded ethanol series (Lee et al., 2008). Dehydrated floral buds were embedded in paraffin and sectioned into 12-µm thick slides. Serial sections of the floral bud tissues were mounted on slide glass and stained with pissophane-hematoxylin. Sectioned floral buds were observed and photographed through microscope (Leica DMI3000B).

RNA Extraction and cDNA Synthesis
Total RNA was extracted from each of the samples using the Quick Plant RNA Isolation Kit (Huayueyang, Beijing, China, Cat. No: 0416-50GK). DNase I (RNase free) (TransGen, Beijing, China, Cat. No: GD201-01) was used at 37 • C for 30 min to isolate gDNA. The RNA integrity was examined by 1% agarose gel. The concentration and quality of the total RNA were measured using Nanodrop TM 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). cDNA was synthesized using 2,000 ng of high-quality total RNA and the PrimeScript TM RT reagent Kit (TaKaRa, Dalian, China, Cat. No: RR037A). The cDNA was synthesized according to the kit's protocol and subsequently diluted 5-fold with molecular biology grade sterile water.

Primer Design
Total, 10 commonly used reference genes, ACT3, ELF1A, G6PD, PEPKR1, TUB, TUA, CYP, GAPDH, H3, and 18S, were identified from the published literature. The corresponding gene sequences were cloned and sequenced based on kenaf transcriptome sequencing data (Chen et al., 2014), and the corresponding gene sequences has been submitted to NCBI (www.ncbi.nlm.nih.gov/). On the basis of the corresponding gene sequences, primers were designed for RT-qPCR using the Primer 3 Plus program (www. primer3plus.com/cgi-bin/dev/primer3plus.cgi). The primers of corresponding gene were tested using Primer Select program within the Lasergene 9 software (Burland, 2000) for dimerization and hairpin formation.

Standard Curve Construction and RT-qPCR
Standard curves were generated using a 5-fold dilution series of the linearized TA-cloning plasmid containing sequences of candidate reference genes and target genes. The correlation coefficient (R 2 ) and slope for each gene were estimated using standard curves. The slope was used to assess the amplification efficiency using the formula E = 10 (−1/slope) −1 (Pfaffl et al., 2004). RT-qPCR reactions were performed in a CFX96 Real Time PCR System (Bio-Rad, California, USA) using a SYBR Premix Ex Taq TM Kit (TaKaRa, Cat. No: RR820A). Each RT-qPCR system contained 1.5 µl of template, 40 nM each primer, and 7.5 µl of SYBR Green PCR master mix, having volume of 15 µl. The RT-qPCR conditions were as follows: 30 s at 95 • C and 40 cycles of 15 s at 95 • C for denaturation and 30 s at 60 • C for annealing. Fluorescence was measured at 65-95 • C. Agarose gel electrophoresis and sequence were performed for the amplified products' specificity analysis. All RT-qPCR reactions were performed in triplicate technical replicates.
For GeNorm and NormFinder, the Ct-values were transformed into relative quantities using the formula 2 − Ct ; Ct is the corresponding Ct-value minus the minimum Ct-value (Chen et al., 2015). GeNorm was used to analyze gene stability (M-value) and pairwise variation (V n / n+1 ) value by pairwise comparisons among all of the reference genes. The M-value correlated with the expression stability, and the lowest M-values were considered the most stable. The threshold of 0.15 was the criterion value of V n / n+1 . A value less than 0.15 are not necessary to use a reference gene for the normalization of the RT-qPCR analysis (Vandesompele et al., 2002). The NormFinder approach is used to rank the candidate reference genes by estimating intra-and inter-group gene expression variations (Andersen et al., 2004). The reference gene with the lowest value was also considered to be the most stable. The Ct approach, which compares standard deviation (SD) values of each gene, was used to rank the candidate reference genes (Silver et al., 2006). BestKeeper ranked the candidate reference genes using the SD and the coefficients of variance (CVs), and the lowest SD value indicated the most stable gene. Based on the ranking from the four statistical approaches, RefFinder (http://www.leonxie.com/ referencegene.php) assigned an appropriate weight to each gene and calculated the geometric mean of their weights for the overall final ranking (Kim et al., 2003).

Validation of Reference Genes
The target gene protein disulfide isomerase-like 5-2a (HcPDIL5-2a), was up-regulated in the CMS line compared with maintainer line of kenaf (Jin, 2011), was used to validate the reliability of the potential reference genes determined by RefFinder and GeNorm approaches across all sample sets. The relative expression levels of HcPDIL5-2a at three stages of pollen development in lines 722A and 722B were normalized with four different reference gene strategies: (1) the most stable reference gene; (2) the geometric mean of the two most stable reference genes; (3) the geometric mean of the three most stable reference genes; and (4) the least stable reference gene. The relative expression of the target gene was calculated according to 2 − Ct method. Student's t-test (P ≤ 0.05) was applied to the variance analysis.

Expression Pattern Analysis of ATP Synthase Subunit Genes
The combination of TUB, CYP, and PEPKR1, which was recommend by the RefFinder and GeNorm methods across all sample sets, was used as an internal control to differentiate the expressions of five mitochondria F 1 F 0 -ATPase subunit genes (atp1, atp4, atp6, atp8, and atp9) in lines 722A and 722B. The expression patterns of these genes were calculated using 2 − Ct method. Least significant differences were determined by Duncan's test (P ≤ 0.05) using SPSS Statistics 21.0 (IBM).

Determining the Stage of Anther Collection Based on a Cytological Analysis
In order to determine the anther collection stage, pollen development at different stages were observed in CMS line 722A and maintainer line 722B of kenaf (Figure 2). The results revealed that the pollen abnormality was initiated after the tetrad stage in line 722A, although no obvious difference was observed between lines 722A and 722B from the pollen mother cell (PMC) to tetrad stage. At the PMC stage, the PMCs were surrounded by the tapetum, and the anther parietal cells were differentiated into the epidermis, endothecium, middle layer and tapetum (Figures 2A1, B1). At the tetrad stage, PMCs formed tetrads after two rounds of meiosis and were wrapped by callose, and tapetum cells expanded and then condensed (Figures 2A2, B2). However, most of the pollen was seriously distorted in line 722A and only a small portion of the pollen was structurally intact compared with line 722B at the monokaryotic stage (Figures 2A3, B3). At the dual-core stage, the microspores were malformed in line 722A (Figure 2A4), although pollen development was normal, the cytoplasm was condensed, and fertile pollen spores were produced in line 722B ( Figure 2B4). Based on the cytological analysis of pollen development in lines 722A and 722B, the anther samples were collected at tetrad, monokaryotic and dual-core stages in lines 722A and 722B,that represented the early, middle and late stages of pollen abortion, respectively, in line 722A.

Evaluation of PCR Primer Specificity and Amplification Efficiency
The specificity of the primer was identified by the presence of a single DNA fragments with the expected size for each primer pair in 3% agarose gel electrophoresis ( Figure S1). These DNA fragments were cloned and sequenced, and the sequences had high identities with expected sequences (Supplementary Material). Additionally, the primer specificity was further validated by the presence of a single peak in the melting curve analysis ( Figure S2). Subsequently the amplification efficiencies (E) and R 2 values were calculated with the standard curve of each primers pair. The E-values of these genes ranged between 0.902 and 1.003, and the R 2 between 0.990 and 1.000 (Table 1). These results demonstrated that the primers of these genes were reliable.

Expression Analysis of the Candidate Reference Gene
To assess the influence of pollen development on the expression stability of the candidate reference genes, the RT-qPCR data of 18 samples were divided into three sets. The first set included all of the samples. The second set contained nine samples from 722A, and the third set was composed of nine samples from 722B.
The transcription levels of 10 candidate reference genes were detected by RT-qPCR assays at three stages of pollen development (tetrad, monokaryotic and dual-core) in lines 722A and 722B of kenaf (Figure 3; Table S1). The Ct-values indicated that the transcripts of these genes were present in a relative widely abundance range, from 6.06 to 33.79, in the 18 tested samples. Among them, TUA showed the lowest transcription abundance, with a mean Ct-value of 28.89, whereas 18S showed the highest abundance, with a mean Ct-value of 7.97. The transcription levels with the least variation were displayed by TUB, CYP, PEPKR1, and H3 (SD < 1.20), while TUA was the most variable (SD = 3.92).

Stability Ranking of Candidate Reference Genes
To evaluate the expression stability of the candidate reference genes at the three stages of pollen development in lines 722A and 722B, four statistical approaches, GeNorm, NormFinder, Ct and BestKeeper, were used.
GeNorm determines stability rankings based on M-values of candidate reference genes. The lowest M-value indicates the most stable expression, while the highest M-value indicates the least stable expression. The most stable candidate reference gene was not the same in the three sample sets, but the least stable candidate reference gene was consistently TUA with M-values of 2.15, 1.97, and 2.04 (Figures 4A,C). In all of the sample sets and the 722B sample set, TUB and CYP were the most stable candidate reference genes, with M-values of 0.44 and 0.13, respectively (Figures 4A,C). In the 722A sample set, ELF1A and H3 were the most stable reference genes, with M-values of 0.26 ( Figure 4B).
NFn+1). In this study, when all of the samples were combined, the pairwise variation values were above the proposed 0.15 cutoff, and V 2/3 was the lowest with V-values of 0.165 ( Figure 4D). This result suggested that at least three candidate reference genes were required for the normalization of a gene expression analysis. However, if the value of V 2/3 was below 0.15 in the 722A and 722B sample sets, then two reference genes were sufficient for normalizing RT-qPCR data ( Figure 4D).
The NormFinder approach calculates gene expression stability based on their minimal combined inter-and intra-group expression variations, which are based on normalization factor calculations. The expression stability of 10 candidate reference genes was calculated with NormFinder ( Table 2). In all of the samples and the 722B sample sets, TUB was the most stable candidate reference genes. In the 722A sample set, ELF1A was the most stable candidate reference genes.
The Ct approach identifies gene expression stability based on the mean SD values of each gene set of pairwise combinations. For all of the samples and the 722B sample set, TUB was the most stable candidate reference gene (Table 3). However, for 722A sample set, H3 was the most stable candidate reference gene ( Table 3). BestKeeper relies on a pairwise correlation analysis of each gene, the lowest SD values indicated the most stable gene. The highest ranking candidate reference gene was H3 in all sample and 722B sample sets (Table 4). ELF1A was the most stably expressed gene with the lowest SD value in 722A sample set ( Table 4).

Comprehensive Stability Ranking Based on RefFinder
RefFinder integrates the above four statistical approaches to produce a comprehensive stability value for each candidate reference gene, and calculates the geometric mean of their  . GeNorm was used to analyze pairwise variation (V n /V n+1 ) between the normalization factors NF n and NF n+1 , V-value < 0.15 denotes that additional reference genes will not significantly improve the accuracy of the results.

Reference Gene Validation
The relative expression level of the target gene HcPDIL5-2a was used to validate the reliability of the reference genes that were recommended by the RefFinder and GeNorm approaches. In the present study, the transcript levels of HcPDIL5-2a at three pollen developmental stages were normalized using four strategies: the one [NF1-1 (TUB)], two [NF1-2 (TUB and CYP)], and three [NF1-3 (TUB, CYP and PEPKR1)] most stable reference genes and the least stable reference gene [NF0 (TUA)] across all sample. At tetrad and monokaryotic stages, the expression levels of HcPDIL5-2a in line 722A were up-regulated relative to line 722B when normalized using NF1-1, NF1-2, NF1-3, and NF0. However, the expression levels of HcPDIL5-2a normalized using NF1-3 and NF0 were significantly different (P < 0.01) (Figure 5). When normalized using NF1-1, NF1-2, and NF1-3 at dual-core stage, the expression level of HcPDIL5-2a in line 722A was upregulated compared with line 722B. In contrast, HcPDIL5-2a's expression was down-regulated when normalized using NF0, and the expression levels of HcPDIL5-2a normalized using NF1-3 and NF0 were significantly different (P < 0.01) (Figure 5).

Differential Expression of Mitochondria F 1 F 0 -ATPase Subunit Genes by RT-qPCR
To explore the effects of the mitochondria F 1 F 0 -ATPase subunit genes (atp1, atp4, atp6, atp8, and atp9) on kenaf CMS, the expression patterns of these genes were differentiated during pollen development between lines 722A and 722B by RT-qPCR. All five genes (atp1, atp4, atp6, atp8, and atp9) exhibited a similar or significantly different expression pattern (Figure 6) in line 722A compared with line 722B after being normalized by the combination of TUB, CYP, and PEPKR1, which was recommend by GeNorm and RefFinder. The expression patterns of atp1 and atp6 were up-regulated gradually (P < 0.05) during pollen development in line 722B, but were down-regulated gradually (P < 0.05) in line 722A (Figure 6). During pollen development, atp9 was down-regulated gradually (P < 0.05) in line 722A, whereas it maintained a steady level (P < 0.05) in line 722B (Figure 6). The expression patterns of atp4 and atp8 were increased (P < 0.05) in lines 722A and 722B (Figure 6).

DISCUSSION
The kenaf CMS line plays a critical role in F1 hybrid seed production and heterosis utilization. Analyses of the expression level of genes of interest that are involved in the regulation of pollen development would aid in our understanding of the molecular mechanisms of CMS in kenaf. RT-qPCR is a useful tool for gene expression studies that relies on a stable reference gene to normalize data. However, many studies have indicated that even the most stable reference gene cannot always be stably expressed under different environmental conditions (Artico et al., 2010). Therefore, it is necessary to assess and validate the stable reference genes under particular experimental conditions instead of using traditional reference genes that have been previously published. The present study is the first systematic survey of the expression stability of 10 candidate reference genes in a CMS and its maintainer lines of kenaf. The expression stability of these genes were evaluated using five statistical approaches, GeNorm, NormFinder, Ct, BestKeeper and RefFinder, to identify the most reliable reference genes. The primer specificities of the target and reference genes should be checked first using several experimental approaches, such as agarose gel electrophoresis, sequencing and melting curves (Bustin et al., 2009;Wang S. et al., 2016). In this study, primer pairs for all 6 target genes and 10 candidate reference genes were validated by agarose gel electrophoresis ( Figure S1), melting curve analysis ( Figure S2) and sequencing. In addition, the primer E and R 2 values were calculated using corresponding standard curves and ranged from 0.902 and 1.042, and 0.990 and 1.000, respectively. These results indicated that all of the material used were suitable for RT-qPCR analysis. Previous research reported that the ranking order of stable reference genes through different statistic approaches . In this study, four commonly used statistic approaches (GeNorm, NormFinder, Ct, and BestKeeper) were used to evaluate the expression stabilities of candidate reference genes. The stability rankings of the 10 candidate reference genes were not the same when using different approaches ( Table 5). For example, TUB was identified as the most stable reference gene by GeNorm, NormFinder and Ct in the sample sets, but H3 was considered to be the most stable reference gene by BestKeeper. Similar results were observed in several studies Fan et al., 2016;Yue et al., 2016). Different statistical approaches, which are based on different principles, may yield contradictory results from the same data FIGURE 6 | Expression pattern analysis of mitochondria F 1 F 0 -ATPase subunit genes atp1, atp4, atp6, atp8, atp9, using RT-qPCR at three stage of pollen development in lines 722A and 722B. The results were shown in relative expression level of the target genes normalized with the most stable three reference genes, TUB, CYP, and PEPKR1. Different lowercase letters indicate significant differences (P < 0.05) in relative expression of each gene among different pollen development stages. The error bars represent standard deviations of three biological replicates. (Hu et al., 2009;Wan et al., 2010). This further confirms that a comprehensive approach (RefFinder) is necessary to evaluate the stability rankings of the 10 candidate reference genes.
It is has been well documented that the accuracy of RT-qPCR results can be significantly improved when multiple reference genes were used (Reid et al., 2006;Gutierrez et al., 2008;Cheng et al., 2017). Thus, in order to explore the optimal number of reference genes at different pollen development stages of kenaf, the Vn/n+1 values which were calculated by GeNorm approach. In this research, all of the Vn/n+1 values were above the proposed 0.15 cut-off across all of the sample (Figure 4D), which indicated that there was no optimal number of reference genes. However, many previous researchers have suggested that three of the most stable genes should be employed in situations when too many, or even no, optimal number of genes were determined (Kuijk et al., 2007;Silveira et al., 2009;Maroufi et al., 2010;Sun et al., 2015). Therefore, TUB, CYP, and PEPKR1, which were recommended by RefFinder, were used as internal controls.
Considering the comprehensive ranking, TUB was identified overall as the most stable candidate reference gene across all of the sample sets. Tubulin plays a critical role in structural support, intracellular transport and DNA segregation. In previously studies, TUB has been confirmed as the most reliable reference genes in CMS7311 and female cabbage during flower development and in Platycladus orientalis under stress conditions (Chang et al., 2012;Xu et al., 2014). By contrast, TUB has been excluded as a good internal control in the sexual tissues of Brachiaria brizantha (Silveira et al., 2009). TUA was the least stable candidate reference gene based on the RefFinder approach in this study, indicating that this gene should be avoided as an internal control when analyzing gene expression patterns in CMS and maintainer lines during pollen development. Our validation experiment indicated that the use of TUA as an internal control led to the misinterpretation of the HcPDIL5-2a expression levels at different stages of pollen development in lines 722A and 722B (Figure 5). By contrast, TUA exhibited a high stability during seed development in peanut (Chi et al., 2012).
GAPDH encodes an important catalyzing enzyme that is involved in the sixth step of glycolysis and has several functions in non-metabolic processes (Tarze et al., 2007;Zala et al., 2013). It was considered as the most stable reference gene under various experimental conditions (Qi et al., 2010) and has been extensively used as a reference gene in kenaf Liao et al., 2016). In contrast, our results indicated that GAPDH was only moderately stable and was not an ideal reference gene for gene expression analyses during pollen development in lines 722A and 722B (Table 3). It is well-known that pollen development/abortion is related to energy metabolism and energy deficiency. The above results may be caused by GAPDH's involvement in energy metabolism, eventually leading to a decrease in expression stability.
Previous studies showed that the application of invalidated reference genes may significantly affect the quantifications of gene expression analyses (Lovdal and Saha, 2014). Even using the single most stable reference gene after assessment, may result in contradictory outcomes (Ding et al., 2015). In the present study, the levels of HcPDIL5-2a expression varied substantially after normalizing with one, two and three of the most and least stable reference genes (Figure 5). The results indicated that using unstable reference genes or a single reference gene is insufficient to obtain reliable and accurate results. Therefore, the expression stability of reference genes should be assessed and multiple reference genes should be used under particular experimental conditions to obtain more accurate and reliable results (Zhu et al., 2012).
Additionally, the combination of TUB, CYP, and PEPKR1 as an internal control was used to normalize the expression pattern of five mitochondria F 1 F 0 -ATPase subunit genes (atp1, atp4, atp6, atp8, and atp9) in lines 722A and 722B. These genes are involved in energy metabolism in plant cells. Their expression patterns in line 722A provided a preliminary impression of the complicated CMS molecular mechanisms in kenaf.
Pollen development is highly energy-consuming process (Lee and Warmke, 1979). Mitochondria are the site of both the tricarboxylic acid (TCA) and oxidative phosphorylation pathway, plays a crucial role in energy and carbon metabolism in eukaryotic cells (Hatefi, 1985). A vital component of all mitochondria is F 1 F 0 -ATPase (Complex V) which is reversibly involved in the synthesis and the hydrolysis of adenosine triphosphate (ATP), depending upon the direction of an electrochemical gradient that is formed by the passage of protons through F 1 F 0 -ATPase (Senior, 1979;Yesodi et al., 1997). The abnormal expression of F 1 F 0 -ATPase subunit gene could cause destruction of membrane electronic potential, which led to mitochondria dysfunction and eventually pollen defect (Li et al., 2010). In the present study, five subunit genes (atp1, atp4, atp6, atp8, and atp9) of mitochondria F 1 F 0 -ATPase had significantly increased or stable expression pattern in line 722B (Figure 6), indicating that these genes had basic function under pollen development for maintaining the normal membrane electronic potential. In contrast, three subunit genes (atp1, atp6, and atp9) of mitochondria F 1 F 0 -ATPase exhibited gradually downregulated expression patterns in line 722A, indicating that the normal membrane electronic potential may be was destroyed, which led to mitochondria dysfunction and eventually pollen defect in kenaf. However, it is necessary to study them further and to verify their correlation with kenaf ' CMS.

CONCLUSION
This study evaluated the expression stability of 10 candidate reference genes across 19 samples using five common statistical approaches. A gene combination (TUB, CYP, and PEPKR1) was predicted to be reliable and was validated using the reference genes to normalize the expression levels of HcPDIL5-2a using RT-qPCR, during the pollen development stage in lines 722A and 722B. Furthermore, the combination of TUB, CYP, and PEPKR1 was used as an internal control for differentiating the expression patterns of five mitochondria F1F0-ATPase subunit genes (atp1, atp4, atp6, atp8, and atp9) during the pollen development in lines 722A and 722B. Three genes (atp1, atp6, and atp9) exhibited significant expression trends in the line 722A compared with line 722B. The data suggest that these genes may be associated with the CMS of kenaf crop. These data will provide useful information for future research on gene expression analysis.

AUTHOR CONTRIBUTIONS
RZ conceived and designed the research. PC and AK revised the manuscript. BZ, YZ, and DL analyzed data and wrote the manuscript; LC, XL, and XK assist during the experiments. All authors read and approved the manuscript.