Front. Plant Sci. Frontiers in Plant Science Front. Plant Sci. 1664-462X Frontiers Media S.A. 10.3389/fpls.2017.01520 Plant Science Original Research Candidate Reference Genes Selection and Application for RT-qPCR Analysis in Kenaf with Cytoplasmic Male Sterility Background Zhou Bujin 1 Chen Peng 1 Khan Aziz 1 Zhao Yanhong 2 Chen Lihong 1 Liu Dongmei 3 Liao Xiaofang 1 Kong Xiangjun 1 Zhou Ruiyang 1 * 1Key Laboratory of Plant Genetics and Breeding, College of Agriculture, Guangxi University Nanning, China 2Cash Crop Institute of Guangxi Academy of Agricultural Sciences Nanning, China 3College of Biological and Food Science, Shangqiu Normal University Shangqiu, China

Edited by: Roger Deal, Emory University, United States

Reviewed by: Yuhui Chen, Noble Research Institute, LLC, United States; Zhao Su, Pennsylvania State University, United States

*Correspondence: Ruiyang Zhou ruiyangzhou@aliyun.com

This article was submitted to Technical Advances in Plant Science, a section of the journal Frontiers in Plant Science

 01 09 2017 2017 8 1520 25 03 2017 18 08 2017 Copyright © 2017 Zhou, Chen, Khan, Zhao, Chen, Liu, Liao, Kong and Zhou. 2017 Zhou, Chen, Khan, Zhao, Chen, Liu, Liao, Kong and Zhou

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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.

 cytoplasmic male sterility (CMS) reference genes RT-qPCR kenaf gene expression pollen development 31571719 YCBZ2015013 National Natural Science Foundation of China10.13039/501100001809 Innovation Project of Guangxi Graduate Education 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 (Bustin et al., 2005; 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 (Xue et al., 2012), ELF1A (Marum et al., 2012; Yue et al., 2016), GAPDH (Zhao et al., 2013), 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; Wang C. et al., 2016) 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 (Zhao et al., 2013). 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 F1F0-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.

 Materials and methods 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.

Phenotypic characterization of the CMS line 722A and its maintainer line 722B's flower buds. (A,C) Phenotypes of the CMS line's floral buds; (B,D) phenotypes of the maintainer line's floral buds. PMCs, The flower buds of pollen mother cell stage; Ts, The flower buds of tetrad stage; Ms, The flower buds monokaryotic pollen stage; Ds, The flower buds of dual-core pollen stage. Bar = 0.5 cm.

 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™ 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). cDNA was synthesized using 2,000 ng of high-quality total RNA and the PrimeScript™ 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 (R2) 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™ 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.

 Candidate reference gene stability analysis

Four common statistical packages, GeNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), ΔCt (Silver et al., 2006), and BestKeeper (Pfaffl et al., 2004), were implemented to estimate the stability of each candidate reference gene.

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 (Vn/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 Vn/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 F1F0-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).

 Results 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.

Micrographs of anthers at different developmental stages in the maintainer line 722B (A1–A4) and CMS line 722A (B1–B4). (A1–B4) Indicate anther of pollen mother cell, tetrad, monokaryotic, and dual-core stage respectively. The Ep, En, ML, T, PMC, Ts, and Ms, Ds represent epidermis, endothecium, middle layer, tapetum, pollen mother cell, microspore of tetrad, monokaryotic and dual-core stage respectively.

 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 R2 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 R2 between 0.990 and 1.000 (Table 1). These results demonstrated that the primers of these genes were reliable.

The amplicon characteristics of primers, 10 candidate reference genes and 6 target genes.

Gene name (abbreviation) Accession NO. Forward primer sequence (5′-3′)/Reverse primer sequence (5′-3′) PCR Efficiency R2 Product size (bp)
18S ribosomal RNA (18S) KX816334 AGAAACGGCTACCACATC/TACTCATTCCAATTACCAGACTC 0.953 0.999 122
Actin 3 (ACT3) KX689256 GTGAGGATATTCAACCCCTTGTCT/CATCTTTCTGTCCCATACCAACC 0.983 0.999 150
Cyclophilin (CYP) KX689255 TCATCTGCACCGCCAAAA/CTTTCTCCACGGCTCTCACC 0.928 0.999 96
Eukaryotic elongation factor 1-alpha (ELF1A) KX703003 GAACATGATCACGGGGACCT/GAGTGAAGGCAAGCAGAGCA 0.932 0.999 128
Glucose-6-phosphate dehydrogenase (G6PD) KX703004 ACGAATTCTCGAAAAGTAGCCAAG/TCCGAGTCCATCCACCAAG 0.972 0.999 141
Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) KX783037 AACGAAAAGGAATACAAGCCAGAG/AAGACCCTCAACAATGCCAAA 0.924 0.999 117
Histone 3 (H3) KX703005 GTGGAGTCAAGAAGCCTCACAG/ATGGCTCTGGAAACGCAAA 0.902 0.999 164
Phosphoenolpyruvate Carboxylase-Related Kinase 1 (PEPKR1) KX703006 TGCCATGAGAATCGCCAAC/GGACACCAACCAAAAGCACA 1.003 0.999 149
Alpha-tubulin (TUA) KX703007 ATTGGCGGAGGTGATGATG/TGGAAGAGTTGGCGGTATGTT 0.936 0.999 140
Beta-tubulin (TUB) KX703008 TTTTCCGACCCGACAACTTC/AGTTCCGCTCCTTCCGTGT 0.930 0.999 85
Protein disulfide isomerase-like5-2a of kenaf (HcPDIL5-2a) HQ638208 CGTTGCTCCTGATGTGTCTATTCT/CTGTCTCATTCAAGCCAAAACCT 0.913 0.999 110
ATP synthase subunit 1 (atp1) KC686836 TAAAAGCGGTAGATAGCCTGG/AATCGCTACATAGACACAATACAA 0.909 0.999 173
ATP synthase subunit 4 (atp4) MF511050 TCTTCCCATCGCATCCGTCTT/TTATTCTTACCCCACCCGAACCAT 0.926 1.000 159
ATP synthase subunit 6 (atp6) HM535783 GCTAATCTCTTATTGTTTTCGCGCA/ATAGCATAGTCCAAGCGAAGCCAC 0.934 1.000 113
ATP synthase subunit 8 (atp8) KM281806 TTCTGGTTATGCCTTTTCTCTT/TGTTTCTTTCTATTCCGCGTGAG 0.942 1.000 319
ATP synthase subunit 9 (atp9)a JX206830 ATGAATGATAAAGCGCGTGACGAG/CGGTTAGAGCAAAGCCCAAAATG 0.995 1.000 208

The primers sequences of this gene was according to Zhao et al. (2013).

 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).

The threshold cycle (Ct) distribution of each candidate reference genes in all experimental samples. A lines across the box of Ct value represent the median values. Lower and upper boxes show the 25 and 75% percentiles. Whiskers represent the maximum and minimum Ct values. The y-axis represents the Ct values of three biotechnology replicates of all experimental samples, while the x-axis represents the 10 candidate reference genes.

 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).

Stability ranking and pairwise variation (V) analysis of the 10 candidate reference genes in different sample sets using GeNorm. The average expression stability values (M) of the 10 candidate reference genes performed in all sample (A), 722A sample set (B), 722B sample set (C), The most stable genes are on the right, and the least stable gene are on the left. The optimal number of reference genes required for accurate normalization (D). GeNorm was used to analyze pairwise variation (Vn/Vn+1) between the normalization factors NFn and NFn+1, V-value < 0.15 denotes that additional reference genes will not significantly improve the accuracy of the results.

Additionally, GeNorm was used to calculate the optimal number of candidate reference genes using the pairwise variation (Vn/n+1) between two sequential normalization factors (NFn and NFn+1). In this study, when all of the samples were combined, the pairwise variation values were above the proposed 0.15 cut-off, and V2/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 V2/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.

Stability ranking of 10 candidate reference gene using NormFinder.

Rank All samples 722A 722B
Gene Stability Gene Stability Gene Stability
1 TUB 0.110 ELF1A 0.067 TUB 0.141
2 CYP 0.113 H3 0.115 CYP 0.172
3 PEPKR1 0.136 TUB 0.120 18S 0.215
4 ELF1A 0.192 PEPKR1 0.174 ELF1A 0.249
5 18S 0.184 CYP 0.205 GAPDH 0.253
6 H3 0.195 ACT3 0.236 PEPKR1 0.252
7 ACT3 0.198 18S 0.254 ACT3 0.271
8 G6PD 0.222 GAPDH 0.316 H3 0.301
9 TUA 0.223 TUA 0.318 TUA 0.338
10 GAPDH 0.247 G6PD 0.363 G6PD 0.468

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).

Stability ranking of 10 candidate reference gene using ΔCt.

Rank All samples 722A 722B
Gene Mean SD Gene Mean SD Gene Mean SD
1 TUB 1.486 H3 1.301 TUB 1.420
2 CYP 1.542 TUB 1.319 CYP 1.440
3 PEPKR1 1.555 ELF1A 1.334 18S 1.497
4 H3 1.796 PEPKR1 1.341 PEPKR1 1.601
5 ELF1A 1.916 CYP 1.423 GAPDH 1.828
6 ACT3 2.158 18S 1.595 H3 1.907
7 GAPDH 2.159 ACT3 1.735 ELF1A 2.112
8 G6PD 3.129 GAPDH 1.988 ACT3 2.182
9 18S 3.504 G6PD 3.201 G6PD 3.091
10 TUA 3.946 TUA 4.460 TUA 3.328

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).

Stability ranking of 10 candidate reference gene using BestKeeper.

Rank All samples 722A 722B
Gene SD Gene SD Gene SD
1 H3 0.4 ELF1A 0.39 H3 0.21
2 CYP 0.66 H3 0.48 GAPDH 0.57
3 TUB 0.7 CYP 0.52 18S 0.64
4 GAPDH 0.96 PEPKR1 0.6 CYP 0.68
5 PEPKR1 0.97 TUB 0.7 TUB 0.69
6 18S 1.02 18S 0.93 PEPKR1 1.26
7 ELF1A 1.38 ACT3 1.22 G6PD 1.73
8 G6PD 1.9 GAPDH 1.33 ELF1A 2.28
9 ACT3 2.01 G6PD 2.08 ACT3 2.33
10 TUA 3.44 TUA 3.73 TUA 3.16
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 ranking for the overall final ranking. For all sample set, TUB, CYP and PEPKR were the three most stable candidate reference genes, and the comprehensive stability rankings were: TUB, CYP, PEPKR1, H3, 18S, ELF1A, GAPDH, ACT3, G6PD, and TUA (Table 5; Table S2). For the 722A sample set, ELF1A and H3 were the two most stable, and the comprehensive stability rankings were: ELF1A, H3, TUB, PEPKR1, CYP, ACT3, 18S, GAPDH, G6PD, and TUA (Table 5; Table S2). For the 722B sample set, TUB and CYP were the two most stable, and the comprehensive stability rankings were: TUB, CYP, 18S, H3, GAPDH, PEPKR1, ELF1A, ACT3, G6PD, and TUA (Table 5; Table S2).

Comprehensive stability ranking of 10 candidate reference genes.

Methods 1 2 3 4 5 6 7 8 9 10
ALL SAMPLES
RefFinder TUB CYP PEPKR1 H3 18S ELF1A GAPDH ACT3 G6PD TUA
GeNorm TUB/CYP PEPKR1 18S H3 ELF1A ACT3 GAPDH G6PD TUA
NormFinder TUB CYP PEPKR1 18S ELF1A H3 ACT3 G6PD TUA GAPDH
ΔCt TUB CYP PEPKR1 H3 ELF1A ACT3 GAPDH G6PD 18S TUA
BestKeeper H3 CYP TUB GAPDH PEPKR1 18S ELF1A G6PD ACT3 TUA
722A
RefFinder ELF1A H3 TUB PEPKR1 CYP ACT3 18S GAPDH G6PD TUA
GeNorm ELF1A/H3 TUB PEPKR1 CYP ACT3 18S GAPDH G6PD TUA
NormFinder ELF1A H3 TUB PEPKR1 CYP ACT3 18S GAPDH TUA G6PD
ΔCt H3 TUB ELF1A PEPKR1 CYP 18S ACT3 GAPDH G6PD TUA
BestKeeper ELF1A H3 CYP PEPKR1 TUB 18S ACT3 GAPDH G6PD TUA
722B
RefFinder TUB CYP 18S H3 GAPDH PEPKR1 ELF1A ACT3 G6PD TUA
GeNorm TUB/CYP 18S PEPKR1 H3 GAPDH ELF1A ACT3 G6PD TUA
NormFinder TUB CYP 18S ELF1A PEPKR GAPDH ACT3 H3 TUA G6PD
ΔCt TUB CYP 18S PEPKR1 GAPDH H3 ELF1A ACT3 G6PD TUA
BestKeeper H3 GAPDH 18S CYP TUB PEPKR1 G6PD ELF1A ACT3 TUA
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 up-regulated 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).

Relative expression levels of a target gene, HcPDIL5-2a, at tetrad, monokaryotic and dual-core stage in lines 722A and 722B using different normalization strategies. NF0: the least stable reference gene, TUA; NF1-1: the most stable reference gene, TUB; NF1-2: the most stable two reference gene, TUB and CYP; NF1-3: the most stable three reference gene, TUB, CYP, and PEPKR1; The error bars represent standard deviations of three biological replicates. Significantly up-regulated with respect to control was shown by **P < 0.01, and significantly down-regulated was shown ##P < 0.01.

 Differential expression of mitochondria F<sub>1</sub>F<sub>0</sub>-ATPase subunit genes by RT-qPCR

To explore the effects of the mitochondria F1F0-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).

Expression pattern analysis of mitochondria F1F0-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.

 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 R2 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 (Wang C. et al., 2016). 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 (Zhang et al., 2015; 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 (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 (Zhao et al., 2013; 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 F1F0-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 F1F0-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 F1F0-ATPase (Senior, 1979; Yesodi et al., 1997). The abnormal expression of F1F0-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 F1F0-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 F1F0-ATPase exhibited gradually down-regulated 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.

Conflict of interest statement

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.

This study was supported by National Natural Science Foundation of China (Grant Number: 31571719) and Innovation Project of Guangxi Graduate Education (Grant Number: YCBZ2015013).

 Supplementary material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017.01520/full#supplementary-material

 References Andersen C. L. Jensen J. L. Orntoft T. F. (2004). Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 64, 5245–5250. 10.1158/0008-5472.CAN-04-049615289330 Artico S. Nardeli S. M. Brilhante O. Grossi-de-Sa M. F. Alves-Ferreira M. (2010). Identification and evaluation of new reference genes in Gossypium hirsutum for accurate normalization of real-time quantitative RT-PCR data. BMC Plant Biol. 10:49. 10.1186/1471-2229-10-4920302670 Burland T. G. (2000). DNASTAR's Lasergene sequence analysis software. Methods Mol. Biol. 132, 71–91. 10.1385/1-59259-192-2:7110547832 Bustin S. A. Benes V. Garson J. A. Hellemans J. Huggett J. Kubista M. . (2009). The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 55, 611–622. 10.1373/clinchem.2008.11279719246619 Bustin S. A. Benes V. Nolan T. Pfaffl M. W. (2005). Quantitative real-time RT-PCR a perspective. J. Mol. Endocrinol. 34, 597–601. 10.1677/jme.1.0175515956331 Chang E. M. Shi S. Q. Liu J. F. Cheng T. L. Xue L. Yang X. Y. . (2012). Selection of reference genes for quantitative gene expression studies in Platycladus orientalis (Cupressaceae) using real-time PCR. PLoS ONE 7:e33278. 10.1371/journal.pone.003327822479379 Chen P. Zhou R. Y. (2011). Cloning and sequnence analysis of cox3 gene in cytoplasmic male sterility and its maintainer liner of Hibiscus cannabinus. Biotechnol. Bull. 22, 37–40. 10.3969/j.issn.1009-0002.2011.01.009 Chen P. Li R. Liao J. Zhou R. Y. (2012). Cloning, expression and characterization of a novel MADS-box gene from kenaf (Hibiscus cannabinus L.). J. Animal Plant Sci. 13, 1714–1724. Available online at: http://www.m.elewa.org/JAPS Chen P. Ran S. M. Li R. Huang Z. P. Qian J. H. Yu M. L. . (2014). Transcriptome de novo assembly and differentially expressed genes related to cytoplasmic male sterility in kenaf (Hibiscus cannabinus L.). Mol. Breed. 34, 1879–1891. 10.1007/s11032-014-0146-8 Chen Y. Hu B. Y. Tan Z. Q. Liu J. Yang Z. M. Li Z. H. . (2015). Selection of reference genes for quantitative real-time PCR normalization in creeping bentgrass involved in four abiotic stresses. Plant Cell Rep. 34, 1825–1834. 10.1007/s00299-015-1830-926179072 Cheng Y. Pang X. Wan H. Ahammed G. J. Yu J. Yao Z. . (2017). Identification of optimal reference genes for normalization of qpcr analysis during pepper fruit development. Front. Plant Sci. 8:1128. 10.3389/fpls.2017.0112828706523 Chi X. Y. Hu R. B. Yang Q. L. Zhang X. W. Pan L. J. Chen M. N. . (2012). Validation of reference genes for gene expression studies in peanut by quantitative real-time RT-PCR. Mol. Genet. Genomics 287, 167–176. 10.1007/s00438-011-0665-522203160 Coito J. L. Rocheta M. Carvalho L. Amancio S. (2012). Microarray-based uncovering reference genes for quantitative real time PCR in grapevine under abiotic stress. BMC Res Notes 5:220. 10.1186/1756-0500-5-22022564373 Ding Y. Sun H. Y. Zhang R. Yang Q. Liu Y. T. Zang X. N. . (2015). Selection of reference gene from Gracilaria lemaneiformis under temperature stress. J. Appl. Phycol. 27, 1–8. 10.1007/s10811-014-0423-2 Fan C. Qiu Z. Zeng B. Liu Y. Li X. Guo G. (2016). Selection of reference genes for quantitative real-time PCR in Casuarina equisetifolia under salt stress. Biol. Plant. 2016, 1–10. 10.1007/s10535-016-0670-y Gutierrez L. Mauriat M. Guenin S. Pelloux J. Lefebvre J. F. Louvet R. . (2008). The lack of a systematic validation of reference genes: a serious pitfall undervalued in REVerse transcription-polymerase chain reaction (RT-PCR) analysis in plants. Plant Biotechnol. J. 6, 609–618. 10.1111/j.1467-7652.2008.0034618433420 Hatefi Y. (1985). The mitochondria electron transport and oxidative phosphorylation system. Annu. Rev. Biochem. 54, 1015–1069. 10.1146/annurev.bi.54.070185.005055 Hu R. B. Fan C. M. Li H. Y. Zhang Q. Z. Fu Y. F. (2009). Evaluation of putative reference genes for gene expression normalization in soybean by quantitative real-time RT-PCR. BMC Mol. Biol. 10:93. 10.1186/1471-2199-10-9319785741 Huggett J. Dheda K. Bustin S. Zumla A. (2005). Real-time RT-PCR normalisation; strategies and considerations. Genes Immun. 6, 279–284. 10.1038/sj.gene.636419015815687 Jin G. (2011). Conversion of SNP Markers Associated with CMS Cytoplasm and Differential Expression Analysis in Anther of CMS and Its Maintainer Line of Kenaf (Hibiscus Cannabinus L.). Nanning: Guangxi University. Jin G. Chen P. Tang X. M. Zhou R. Y. (2011). synchronous cloning of 5' cDNA end of two PDIL homologous gene from kenaf (Hibiscus cannabinus L.) using degenerate primers. Biotechnology Bulletin. 10, 109–113. 10.13560/j.cnki.biotech.bull.1985.2011.10.008 Kim B. R. Nam H. Y. Kim S. U. Kim S. I. Chang Y. J. (2003). Normalization of reverse transcription quantitative-PCR with housekeeping genes in rice. Biotechnol. Lett. 25, 1869–1872. 10.1023/A:102629803200914677714 Kuijk E. W. du Puy L. van Tol H. T. Haagsman H. P. Colenbrander B. Roelen B. A. (2007). Validation of reference genes for quantitative RT-PCR studies in porcine oocytes and preimplantation embryos. BMC Dev. Biol. 7:58. 10.1186/1471-213X-7-5817540017 Lee S. L. J. Warmke H. E. (1979). Organelle size and number in fertile and T-cytoplasmic male-sterile corn. Amer. J. Bot. 66, 141–148. jstor: stable/2442516 Lee Y. P. Park S. Lim C. Kim H. Lim H. Ahn Y. . (2008). Discovery of a novel cytoplasmic male-sterility and its restorer lines in radish (Raphanus sativus L.). Theor. Appl. Genet. 117, 905–913. 10.1007/s00122-008-0830-318597066 Li W. Q. Zhang X. Q. Xia C. Deng Y. Ye D. (2010). MALE GAMETOPHYTE DEFECTIVE 1, encoding the FAd subunit of mitochondrial F1F0-ATP synthase, is essential for pollen formation in Arabidopsis thaliana. Plant Cell Physiol. 51, 923–935. 10.1093/pcp/pcq06620435646 Li Y. Liu T. K. Duan W. K. Song X. M. Shi G. J. Zhang J. Y. . (2014). Instability in mitochondrial membranes in Polima cytoplasmic male sterility of Brassica rapa ssp. chinensis. Funct. Integr. Genomics 14, 441–451. 10.1007/s10142-014-0368-124652098 Liao X. F. Zhao Y. H. Chen P. Zhou B. J. Diao Y. Yu M. L. . (2016). A comparative analysis of the atp8 gene between a cytoplasmic male sterile line and its maintainer and further development of a molecular marker specific to male sterile cytoplasm in kenaf (Hibiscus cannabinus L.). Plant Mol. Biol. Rep. 34, 29–36. 10.1007/s11105-015-0894-6 Linke B. Nothnagel T. Borner T. (2003). Flower development in carrot CMS plants: mitochondria affect the expression of MADS box genes homologous to GLOBOSA and DEFICIENS. Plant J. 34, 27–37. 10.1046/j.1365-313X.2003.01703.x12662306 Lovdal T. Saha A. (2014). Reference gene selection in Carnobacterium maltaromaticum, Lactobacillus curvatus, and Listeria innocua subjected to temperature and salt stress. Mol. Biotechnol. 56, 210–222. 10.1007/s12033-013-9697-x24037409 Maroufi A. Van Bockstaele E. De Loose M. (2010). Validation of reference genes for gene expression analysis in chicory (Cichorium intybus) using quantitative real-time PCR. BMC Mol. Biol. 11:15. 10.1186/1471-2199-11-1520156357 Marum L. Miguel A. Ricardo C. P. Miguel C. (2012). Reference gene selection for quantitative real-time PCR normalization in Quercus suber. PLoS ONE 7:e35113. 10.1371/journal.pone.003511322529976 Pfaffl M. W. Tichopad A. Prgomet C. Neuvians T. P. (2004). Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper–Excel-based tool using pair-wise correlations. Biotechnol. Lett. 26, 509–515. 10.1023/B:BILE.0000019559.84305.4715127793 Qi J. N. Yu S. C. Zhang F. L. Shen X. Q. Zhao X. Y. Yu Y. J. . (2010). Reference gene selection for real-time quantitative polymerase chain reaction of mRNA transcript levels in Chinese Cabbage (Brassica rapa L. ssp. pekinensis). Plant Mol. Biol. Rep. 28, 597–604. 10.1007/s11105-010-0185-1 Qian J. H. Zhou B. J. Kong X. J. Li Z. Q. Shi Q. Q. Liao X. F. . (2016). Cloning, expression analysis and vector construction of MYB21 gene in Kenaf. Biotechnol. Bullet. 32, 170–179. 10.13560/j.cnki.biotech.bull.1985.2016 Radonic A. Thulke S. Mackay I. M. Landt O. Siegert W. Nitsche A. (2004). Guideline to reference gene selection for quantitative real-time PCR. Biochem. Biophys. Res. Commun. 313, 856–862. 10.1016/j.bbrc.2003.11.17714706621 Reid K. E. Olsson N. Schlosser J. Peng F. Lund S. T. (2006). An optimized grapevine RNA isolation procedure and statistical determination of reference genes for real-time RT-PCR during berry development. BMC Plant Biol. 6:27. 10.1186/1471-2229-6-2717105665 Schnable P. S. Wise R. P. (1998). The molecular basis of cytoplasmic male sterility and fertility restoration. Trends Plant Sci. 3, 175–180. 10.1016/S1360-1385(98)01235-7 Senior A. (1979). The mitochondria ATPase, in Membrane Proteins in Energy Transduction. ed Capaldi R. (New York, NY: Dekker), 233–278. Silveira E. D. Alves-Ferreira M. Guimaraes L. A. Da S. F. Carneiro V. T. (2009). Selection of reference genes for quantitative real-time PCR expression studies in the apomictic and sexual grass Brachiaria brizantha. BMC Plant Biol. 9:84. 10.1186/1471-2229-9-8419573233 Silver N. Best S. Jiang J. Thein S. L. (2006). Selection of housekeeping genes for gene expression studies in human reticulocytes using real-time PCR. BMC Mol. Biol. 7:33. 10.1186/1471-2199-7-3317026756 Sun M. Lu M. X. Tang X. T. Du Y. Z. (2015). Exploring valid reference genes for quantitative real-time PCR analysis in Sesamia inferens (Lepidoptera: Noctuidae). PLoS ONE 10:e0115979. 10.1371/journal.pone.011597925585250 Sun X. J. Wang Y. Chen Y. L. Xu L. Jiang L. N. Gong Y. Q. . (2013). Differential gene expression profiling of Ogura CMS line and its maintainer in radish (Raphanus sativus L.). Acta Physiol. Plant. 35, 3413–3425. 10.1007/s11738-013-1376-9 Tarze A. Deniaud A. Le Bras M. Maillier E. Molle D. Larochette N. . (2007). GAPDH, a novel regulator of the pro-apoptotic mitochondrial membrane permeabilization. Oncogene 26, 2606–2620. 10.1038/sj.onc.121007417072346 Tong Z. G. Gao Z. H. Wang F. Zhou J. Zhang Z. (2009). Selection of reliable reference genes for gene expression studies in peach using real-time PCR. BMC Mol. Biol. 10:71. 10.1186/1471-2199-10-7119619301 Udvardi M. K. Czechowski T. Scheible W. R. (2008). Eleven golden rules of quantitative RT-PCR. Plant Cell 20, 1736–1737. 10.1105/tpc.108.06114318664613 Vandesompele J. De Preter K. Pattyn F. Poppe B. Van Roy N. De Paepe A. . (2002). Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 3:RESEARCH0034. 10.1186/gb-2002-3-7-research003412184808 Wan H. J. Zhao Z. G. Qian C. T. Sui Y. H. Malik A. A. Chen J. F. (2010). Selection of appropriate reference genes for gene expression studies by quantitative real-time polymerase chain reaction in cucumber. Anal. Biochem. 399, 257–261. 10.1016/j.ab.2009.12.00820005862 Wang C. Cui H. M. Huang T. H. Liu T. K. Hou X. L. Li Y. (2016). Identification and validation of reference genes for RT-qPCR analysis in non-heading Chinese Cabbage Flowers. Front. Plant Sci. 7:811. 10.3389/fpls.2016.0081127375663 Wang S. Wang C. Zhang X. X. Chen X. Liu J. J. Jia X. F. . (2016). Transcriptome de novo assembly and analysis of differentially expressed genes related to cytoplasmic male sterility in cabbage. Plant Physiol. Bioch. 105, 224–232. 10.1016/j.plaphy.2016.04.02727116370 Wei L. B. Miao H. M. Zhao R. H. Han X. H. Zhang T. D. Zhang H. Y. (2013). Identification and testing of reference genes for Sesame gene expression analysis by quantitative real-time PCR. Planta 237, 873–889. 10.1007/s00425-012-1805-923229061 Xu X. Y. Yang Z. P. Sun X. L. Zhang L. G. Fang Z. Y. (2014). Selection of reference genes for quantitative real-time PCR during flower bud development in CMS7311 of heading Chinese cabbage (Brassica rapa L. ssp. pekinensis). Acta Physiol. Plant. 36, 809–814. 10.1007/s11738-013-1437-0 Xue J. Qiao N. Zhang W. Cheng R. L. Zhang X. Q. Bao Y. Y. . (2012). Dynamic interactions between Bombyx mori nucleopolyhedrovirus and its host cells revealed by transcriptome analysis. J. Virol. 86, 7345–7359. 10.1128/JVI.07217-1222532689 Yeemin K. Sariah M. Chailing H. Muiyun W. (2016). Selection of reference genes for quantitative real-time PCR normalization in Ganoderma-infected oil palm (Elaeis guineensis) seedlings. Australas. Plant Path. 45, 261–268. 10.1007/s13313-016-0417-4 Yesodi V. Hauschner H. Tabib Y. Firon N. (1997). An intact F1ATPase α-subunit gene and a pseudogene with differing genomic organization are detected in both male-fertile and CMS petunia mitochondria. Curr. Genet. 32, 348–357. 10.1007/s0029400502879371886 Yue H. Deng P. C. Liu S. Y. Wang M. Song W. N. Nie X. J. (2016). Selection and evaluation of reference genes for quantitative gene expression analysis in broomcorn millet (Panicum miliaceum L.). J. Plant Biol. 59, 435–443. 10.1007/s12374-016-0024-5 Zala D. Hinckelmann M. V. Yu H. Lyra D. C. M. Liot G. Cordelieres F. P. . (2013). Vesicular glycolysis provides on-board energy for fast axonal transport. Cell 152, 479–491. 10.1016/j.cell.2012.12.02923374344 Zhang S. D. An S. H. Li Z. Wu F. M. Yang Q. P. Liu Y. C. . (2015). Identification and validation of reference genes for normalization of gene expression analysis using qRT-PCR in Helicoverpa armigera (Lepidoptera: Noctuidae). Gene 555, 393–402. 10.1016/j.gene.2014.11.03825447918 Zhao Y. H. Chen P. Liao X. F. Zhou B. J. Liao J. Huang Z. P. . (2013). A comparative study of the atp9 gene between a cytoplasmic male sterile line and its maintainer line and further development of a molecular marker specific for male sterile cytoplasm in kenaf (Hibiscus cannabinus L.). Mol. Breed. 32, 969–976. 10.1007/s11032-013-9926-9 Zhou R. Y. Zhang X. Zhang J. Q. Gan Z. H. Wei H. X. (2008). A breakthrough in Kenaf cytoplasmic male sterile lines breeding and heterosis utilization. Sci. Agric. 41, 313–314. Zhu X. Y. Li X. P. Chen W. X. Chen J. Y. Lu W. J. Chen L. . (2012). Evaluation of new reference genes in papaya for accurate transcript normalization under different experimental conditions. PLoS ONE 7:e44405. 10.1371/journal.pone.004440522952972