Selection and Validation of Appropriate Reference Genes for Quantitative Real-Time PCR Normalization in Staminate and Perfect Flowers of Andromonoecious Taihangia rupestris

Quantitative real-time reverse transcription-polymerase chain reaction (qRT-PCR) is the most commonly used and powerful method for gene expression analysis due to its high sensitivity, specificity, and high throughput, and the accuracy of this approach depends on the stability of reference genes used for normalization. Taihangia rupestris Yu and Li (Rosaceae), an andromonoecious plant, produces both bisexual flowers and unisexual male flowers within the same individual. Using qRT-PCR technique, investigation of the gene expression profiling in staminate and perfect flowers would improve our understanding of the molecular mechanism in regulation of flower formation and sex differentiation in andromonoecious T. rupestris. To accurate normalize the gene expression level in Taihangia flower, 16 candidate reference genes, including 10 traditional housekeeping genes, and 6 newly stable genes, were selected based on transcriptome sequence data and previous studies. The expressions of these genes were assessed by qRT-PCR analysis in 51 samples, including 30 staminate and perfect flower samples across developmental stages and 21 different floral tissue samples from mature flowers. By using geNorm, NormFinder, BestKeeper, and comprehensive RefFinder algorithms, ADF3 combined with UFD1 were identified as the optimal reference genes for staminate flowers, while the combination of HIS3/ADF3 was the most accurate reference genes for perfect floral samples. For floral tissues, HIS3, UFD1, and TMP50 were the most suitable reference genes. Furthermore, two target genes, TruPI, and TruFBP24, involved in floral organ identity were selected to validate the most and least stable reference genes in staminate flowers, perfect flowers, and different floral tissues, indicating that the use of inappropriate reference genes for normalization will lead to the adverse results. The reference genes identified in this study will improve the accuracy of qRT-PCR quantification of target gene expression in andromonoecious T. rupestris flowers, and will facilitate the functional genomics studies on flower development and sex differentiation in the future.


INTRODUCTION
Quantitative real-time reverse transcription-polymerase chain reaction (qRT-PCR) has become the most common and powerful method for detecting and measuring transcripts abundance due to its high sensitivity, specificity, accuracy, and high throughput (Gachon et al., 2004). However, to obtain reliable and accurate results, qRT-PCR requires specific strategies to control for possible variability related to the serials of steps of the experimental procedure, such as discrepancy in initial sample amount, quantity of mRNA templates, and qRT-PCR amplification efficiency (Bustin, 2002;Huggett et al., 2005). Relative quantification is a widely accepted procedure to evaluate gene expression by normalization with one or more internal reference genes, which are required for stable expression regardless of experimental conditions Thellin et al., 2009).
Traditionally, the housekeeping genes, such as actin (ACT), elongation factor 1 alpha (EF-1α), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), and ubiquitin (UBQ), involved in primary metabolism or other basic cellular processes are regarded as steadily expressed genes, (Thellin et al., 1999;Kozera and Rapacz, 2013), and thus these genes are commonly used for normalizing qPCR data without any validation. For non-model plant species, the reference genes are mainly identified by the search for orthologous sequences of common housekeeping genes reported in model plant species, due to limited genetic and sequence information (Gonzalezaguero et al., 2013). However, numerous studies have reported that the expression levels of traditional housekeeping genes vary considerably alterations in different biological sample and variable experimental condition (Hu et al., 2009;Zhang et al., 2017). Consequently, the use of improper traditional housekeeping genes for normalization in relative quantification of gene expression could result in bias of expression data and incorrect conclusion (Gutierrez et al., 2008).
With the development of high-throughput sequencing technology, a large number of transcriptome sequence data provide abundant information about gene expression profiling in many plant species. Benefitting from this technique, the growing plant transcriptome datasets have been widely used as new resource for identification of appropriate reference genes, especially in non-model plant (Zhuang et al., 2015;Ma et al., 2016). Furthermore, several statistical algorithms, such as geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), BestKeeper (Pfaffl et al., 2004), and RefFinder (Xie et al., 2012) have been developed to screen the most appropriate reference gene for qRT-PCR normalization in given biological samples. Therefore, next generation sequencing technology such as high-throughput RNA-seq provides new opportunities to explore appropriate reference gene candidates for non-model plant species (Demidenko et al., 2011).
Taihangia rupestris Yu and Li (Rosaceae), an ancient perennial herb belonging to the tribe Dryadeae, is distributed on cliff faces in the southern part of the Taihang Mountains of China (Yu and Li, 1980). T. rupestris produces both bisexual flowers and unisexual male flowers within the same individual, forming the andromonoecious sexual system (Yu and Li, 1983). The Taihangia staminate flowers are bisexual at initiation and become unisexual by arresting the pistil development in consequent developmental stages (Lu, 1996). By using qRT-PCR technique investigation of the gene expression profiling in staminate and perfect flowers would improve our understanding of the molecular mechanism in regulation of flower formation and sex differentiation in T. rupestris. In previous studies, a housekeeping gene Actin was used as the reference gene for normalizing qRT-PCR analysis in T. rupestris flowers (Du et al., 2008;Lu et al., 2010). However, little information is available concerning appropriate reference genes for qRT-PCR normalization, and no systematic selection of reference genes have not been performed in this plant species. Recently, we used de novo RNA sequencing to compare the transcriptome profiles of staminate and perfect flowers at early and late developmental stages (Li et al., 2017), and transcriptome sequences data available could be a potential source for identification of suitable reference genes in andromonoecious T. rupestris.
To accurate normalize the gene expression level in T. rupestris flower, 15 candidate reference genes, including nine traditional housekeeping genes and six novel genes, were selected based on the floral transcriptome datasets. In addition, the housekeeping genes Actin used as reference gene in previous studies was also included. The expression stability of these candidate genes were assessed by qRT-PCR analysis in 30 staminate and perfect flower samples across developmental stages and 21 different floral tissue samples from mature flowers. Moreover, two MADS-box genes TruPI and TruFBP24 involved in control of floral organ identity were used to validate the stable reference genes selected in this study. The aim of this work was to identify reference genes appropriate for transcript normalization in andromonoecious Taihangia flowers, and then validate their stability in staminate flowers, perfect flowers, and different floral tissues. This work will provide the basis for further research in exploring the molecular mechanism of flower development and sex differentiation in andromonoecuous plants.

Plant Material
The plant materials were transplanted from Zhuyufeng (35 • 27 ′ N, 113 • 22 ′ E), Henan, China in April 2014, and were grown in the greenhouse of Henan Polytechnic University. During flowering phrase in early spring of 2016, the staminate and perfect flowers at different developmental stages were harvested from the same plant, respectively ( Figure S1). Flowers were classified into five developmental stages: young flower buds (stage 1, < 0.5 cm), elongated buds (stage 2, 0.5-1 cm), pre-anthesis (stage 3), fully opened flowers (stage 4), and mature flowers (stage 5). Three independent samples of staminate and perfect flowers were taken from six different individuals at each of the developmental stage, respectively. In addition, floral tissues including sepals, petals, stamens, and carpals were collected from fully opened staminate and perfect flowers, respectively. All samples were immediately frozen in liquid nitrogen and stored at −80 • C for RNA isolation.

Total RNA Isolation and cDNA Synthesis
Total RNA was extracted from 100 mg homogenized flowers and floral tissues using an RNeasy Mini Kit (Qiagen, Germany) according to the manufacturer's manuals. RNA quality was initially assessed by 1.5% agarose gel electrophoresis. After that, RNA was quantified with a UV-visible spectrophotometer (UV-2550, Japan), and were adjusted to the same concentration for each sample after measuring the RNA concentration. The RNA samples with A260/A280 ratios between 1.8 and 2.0 and A260/A230 ratios greater than 2.0 were used for subsequent analyses. To eliminate DNA contamination, 1 µg total RNA was digested using gDNA Eraser (Takara, Japan) at 42 • C for 2 min in a 10 µl mixture. Then, the RNA sample was used to synthesis the first strand cDNA with the Perfect Real Time RT reagent Kit (Takara, Japan) in a 20 µl reaction mixture according to the manufacturer's protocol. The cDNA samples were diluted 1:10 with nuclease-free water prior to the qRT-PCR analysis.

Identification and Selection of Candidate Reference Genes
In our previous study, four cDNA libraries, including male floral bud, hermaphroditic floral bud, male flower, and hermaphroditic flower, were constructed and sequenced by using the Illumina RNA-Seq method. A total of 24,753 unigenes were annotated in the NCBI non-redundant protein database, and further used for mining candidate reference genes based on expression stability (Li et al., 2017). To estimate expression stability of each gene, the values of MV, SD, and CV, were calculated for each gene based on read count, according to the previous described methodology (Czechowski et al., 2005;Gonzalezaguero et al., 2013). The genes that had both a mean of read counts above 500 and a CV below 0.3 were considered to be stably expressed. To facilitate the identification of putative function, the stably expressed genes were further screened with E-value ≤ 10 −50 . Candidate reference genes were selected from the homologous of traditional housekeeping genes previous used for flower development according to gene Nr annotation (Mallona et al., 2010;Yeap et al., 2013;Li et al., 2016;Wang et al., 2016). Meanwhile, some novel genes with a minor variation in expression level were also considered. In addition, the housekeeping gene Actin (ACT2) used as reference gene in previous studies was included (Du et al., 2008).

PCR Primer Design and qRT-PCR Analysis
Primers were designed using the Primer Premier 5 software according to the following parameters: primer lengths 18-23 bp, GC content 40-60%, melting temperatures 55-60 • C, and amplicon lengths 100-300 bp. In order to confirm the primer specificity, the specificities of all primer pairs were initially tested by standard RT-PCR, and amplification product for each gene was verified by electrophoresis on a 1.5% agarose gel followed by ethidium bromide staining.
The qRT-PCR reactions were performed with the MiniOpticon Real-Time Detection System (Bio-Rad) using the SsoFast EvaGreen Supermix RT-PCR kit (Bio-Rad Laboratories). The PCR reaction mixture (20 µl) contained 10 µl of EvaGreen Supermix, 2.0 µl of 1:10 diluted cDNA, 0.4 µl of each primer (10 mM), and 7.2 µl of water. The reactions were incubated under following cycling conditions: 2 min at 95 • C, 40 cycles of 95 • C for 15 s, a specific annealing temperature (Ta) for 15 s, and 72 • C for 30 s, and finally 72 • C for 2 min with a single melt cycle from 60 to 95 • C in 5 s intervals. The specificity of primer pair was verified by the presence of a single peak in the melt curve analysis during qRT-PCR. Three independent biological replicates and three technical repetitions were performed for each of the quantitative PCR experiments.
The threshold cycle (Ct) was measured automatically, and was used to define the expression level of each reference gene. The correlation coefficients (R 2 ) and the slope were calculated from standard curve based on 5-fold series dilution of the cDNA templates, and the corresponding qRT-PCR efficiencies (E) for each gene were determined from the given slope (Ginzinger, 2002).

Data Analysis
The Ct values were used to evaluate the stability of 16 candidate reference genes in staminate flowers, perfect flowers, floral tissues, and total floral samples, across developmental stages by using geNorm, NormFinder, and BestKeeper algorithms. Furthermore, the web-based tool RefFinder (http://150.216.56. 64/referencegene.php), which integrates these methods, was used to evaluate and screen the appropriate reference genes for qRT-PCR normalization.

Validation of Reference Genes
To confirm the reliability of the candidate reference genes, the relative expression profiles of two MADS-box transcription factors TruPI and TruFBP24, involved in flower organ identity (Krizek and Meyerowitz, 1996;De Folter et al., 2006), were detected by using RT-qPCR analysis. The expressions of the target genes TruPI and TruFBP24 were normalized using with the most and the least stable reference genes, respectively, as suggested by RefFinder. The qRT-PCR amplification conditions were the same as described above. The relative expression of each gene was calculated using the 2 − CT method (Livak and Schmittgen, 2001).

Selection of Candidate Reference Genes and PCR Amplification
Based on transcriptome data from Taihangia staminate and perfect flowers, a total of 3,628 genes were identified as stable expressed genes with average read counts above 500 and CV< 0.3 (Supplementary Material Data Sheet 1). Based on Nr annotation, we identified nine traditional housekeeping genes, such as GAPDH, EF-1α, UBQ, and ACT, common used as reference genes for flower development (E-value ≤ 10 −50 ). Meanwhile, we selected six novel reference genes including actin-depolymerizing factor 3 (ADF3), ubiquitin fusion degradation protein 1 (UFD1), iron-sulfur cluster assembly protein (ISP), thiosulfate/3-mercaptopyruvate sulfurtransferase (THS), transmembrane protein 50 (TMP50), and vacuolar protein sorting-associated protein 32 (VAP32), with relatively stable expression profiles among staminate and perfect flower across the developmental stages. In addition, the housekeeping gene Actin (ACT2) used as the reference gene in previous Taihangia studies was also included. The results showed that 16 candidate reference genes involved in many aspects of primary metabolism or other basic cellular processes such as cytoskeleton (ACT, ACT2, ADF3, and TUA), transport of ions (THS), transport in vacuoles (VAP32) or membranes (TMP50), glycolysis (GAPDH), energy metabolism (ISP), RNA processing and modification (SPF), protein translation (EF-1α), protein posttranslational modification (UBC, UFD1, and UBQ), cell signaling (PP2A), and chromatin structure, and dynamics (HIS3). Based on Taihangia transcriptome sequence data, we found varied transcriptional profiles for reference gene candidates. Average read count of each gene ranged from 703.25 (ACT) to 7,940.25 (GAPDH), while the CV value varied from 0.007(UFD1) to 0.203 (GAPDH).
The specificity of RT-PCR amplification for each reference gene was verified by 1.5% agarose gel electrophoresis and melting curve analysis ( Figure S2). The results showed that all 16 primer pairs generated single fragments with expected size, ranging from 100 bp for ACT2 to 254 bp for ISP, indicating that no primer dimers or non-specific products amplified in RT-PCR reactions. The amplification specificity of each reference gene was further confirmed by the presence of single-peak melting curve in qRT-PCR analysis. For each primer pair, the amplification efficiencies of ranged from 90.4% for PP2A to 104.5% for TUA, while correlation coefficients (R 2 ) ranged from 0.9901 (ACT2) to 0.9999 (ISP). The candidate reference genes, primer sequences, and characteristics of PCR amplifications were summarized in Table 1.

Expression Stability of the Candidate Reference Genes
To detect overall expression profiles of the 16 reference genes in Taihangia flowers, the 30 floral samples, including staminate and perfect flowers at different developmental stages, and 21 floral tissue samples (sepal, petal, stamen, and carpel) were used for the RT-qPCR analyses. The results showed that the mean Ct values of the 16 reference genes varied from 19.11 ± 0.85 (HIS3) (mean ± SD) to 23.96 ± 0.77 (SPF) for highest and lowest expression levels, respectively, in all samples (Figure 1). In addition, standard deviation (SD) of Ct values, which represented as expression stability, ranged from 0.63 to 1.45. The genes with higher SD of Ct values indicated more variable expression To further select the most appropriate reference gene for qRT-PCR analysis in investigation of flower development and sex differentiation within andromonoecious Taihangia, we divided these floral samples into four groups: staminate flower, perfect flower, floral tissues, and total samples, and four statistical approaches (geNorm, NormFinder, BestKeeper, and RefFinder) were used to analyze the stability of each reference gene from different floral groups.

Genorm Analysis
The geNorm program was employed to assess the expression stability (M) for each reference gene based on the ratio of average pair-wise variation within all gene candidates. According to geNorm algorithm, the reference gene with M value below 1.5 is considered to be stably expressed, and lower M value represents the more stable expression (Vandesompele et al., 2002). In this study, all of the tested candidate genes showed an M value < 1.5, indicating that the 16 selected genes should be considered relatively stable. ADF3 and HIS3 were the most stable reference genes for both staminate and perfect flowers, while UBQ and SPF were identified as the most stable in sepal, petal, stamen, and carpel samples. For all samples tested, HIS3 and UFD1 were recommended as the most stable genes. In contrast, TUA in staminate flowers and GAPDH in perfect flowers, floral tissues, and total samples with the highest M value were identified as the least stable reference genes, respectively (Figure 2).
The optimal number of the reference genes was also determined based on the pairwise variation between sequential ranked genes (Vn/Vn+1) with the cut-off value of 0.15 in geNorm program. If the value of Vn/n+1 was below 0.15, it is not necessary to append additional genes for accurate normalization.
The V2/V3 values were below 0.15 in both staminate and perfect flowers, suggesting that adding an extra gene to obtain accurate results was not necessary for normalization. For the floral tissues and total samples, V2/3 value was 0.219 and 0.185 showed that additional reference genes may be required (Figure 3).

NormFinder Analysis
As a mathematical model-based algorithm, NormFinder is widely applied in determination of the stability value of reference genes, based on inter-and intra-group variance estimation in given serials of specific samples (Andersen et al., 2004). According to the stability value for each reference gene, ranking order was arranged from smallest to largest, with the lower stability value indicating the higher stability. In staminate flowers, UFD1 was ranked first for expression stability, whereas TUA with the largest stability value was the least stably expressed. For perfect flowers and floral tissues samples, HIS3 was identified as the most suitable reference gene, while GAPDH was the most unstably expressed. When total samples were taken into considered, UFD1 was ranked first as the most stable expressed gene, followed by HIS3 and ADF3 ( Table 2).

BestKeeper Analysis
BestKeeper program evaluates the stability of reference genes by means of the standard deviation (SD) and coefficient of variance (CV) of the Ct values, with the lowest SD value indicating the most stable reference gene (Pfaffl et al., 2004). The genes with a SD [± CP] value below 1.0 are considered to be stably expressed and used for the gene expression normalization. According to these criteria, UFD1 was recommended as the most stable gene in staminate flowers, floral tissues, and all samples, while SPF was the most stable in perfect flowers. On the contrary, TUA in staminate flower, PP2A in perfect flower, and GAPDH in floral tissues and total samples with the highest SD value displayed the least stable expression as determined by Bestkeeper software, respectively ( Table 3).

RefFinder Analysis
RefFinder was applied to calculate the geometric mean of weights for the comprehensive ranking order recommended by Ct, geNorm, NormFinder, and BestKeeper (Xie et al., 2012). Based on RefFinder analysis, UFD1 and ADF3 were the optimal reference genes for staminate flowers, while HIS3 and ADF3 were the most stable in perfect floral samples. For floral tissues samples tested, HIS3 and UFD1 were recommended as the suitable reference genes. When total samples were taken into considered, two novel genes (UFD1 and ADF3) and one traditional housekeeping gene (HIS3) were identified as the most stable expressed genes. Conversely, several traditional housekeeping genes showed unstable expression profiles in T. rupestris flowers. GAPDH was the least stable in floral tissues, perfect flowers, and total samples, while TUA was unstable for staminate flowers. The rankings of the four algorithms were integrated by RefFinder and the results are shown in Table 4.

Validation of the Reference Genes
In order to examine reliability of reference genes for normalization, the relative expression patterns of two MADS-box genes TruPI and TruFBP24, which are required to specify flower organ identity, were evaluated in floral tissues, staminate and perfect flowers using qRT-PCR analysis (Krizek and Meyerowitz, 1996;De Folter et al., 2006). According to geNorm pairwise variation (Vn/Vn+1) in staminate flowers, perfect flowers, and floral tissues, we selected the most stable reference genes (UFD1, ADF3, and UFD1/ADF3 for staminate flowers, HIS3, ADF3, and HIS3/ADF3 for perfect flowers, and HIS3/UFD1, HIS3/TMP50, and HIS3/UFD1/TMP50 for floral tissues, respectively) to normalize TruPI gene expression in staminate and perfect flowers across developmental stages. Meanwhile, the expression levels of TruPI were also investigated in different floral tissues of staminate and perfect flowers such as sepal, petal, stamen and carpel. Moreover, we evaluated target genes expression following normalization with the least stable reference gene (TUA for staminate and GAPDH for perfect flowers and floral tissues, respectively) for a comparative analysis. When the stable reference gene(s) were used for normalization, TruPI transcripts increased at early developmental stage with maximum at stage 2, and then steadily decreased with developmental stage in both staminate and perfect flowers (Figure 4). By using different stable reference genes, either single or the combination, the expression patterns of TruPI were similar with minor differences at stage 2. Based on geNorm pairwise variation, the combination of UFD1/ADF3 was recommended to be the optimum pairs of reference gene for staminate flowers, and HIS3/ADF3 was the suitable pairs of reference gene for accurate normalization in perfect flowers. However, the relative expression levels of target gene TruPI were underestimated in Taihangia flowers at developmental stage 2 and 3, when normalized using the most unstable reference genes. Normalization by TUA in staminate and GAPDH in perfect flowers showed decrease trends in the relative expression levels of TruPI at early developmental stages, indicating that the adverse effect of using an inappropriate reference gene. For floral tissues of staminate and perfect flowers, TruPI was expressed in all of the floral organs, including sepals, petals, stamens, and carpels, but at different levels ( Figure 5). When using the combination of HIS3/UFD1, HIS3/TMP50, and HIS3/UFD1/TMP50 for normalization, the qRT-PCR analyses showed that TruPI was strongly expressed petals and stamens, and weakly in sepals and carpels. However, different expression profiles of TruPI were observed in floral tissues, especially in petals and stamens, by using GAPDH for normalization. Despite the higher expression levels of TruPI also detected in petals and stamens, the overestimation of relative expression was observed.
To further explore the stability of reference genes for normalization qRT-PCR between staminate and perfect flowers, the expression levels of TruFBP24 was investigated for all floral samples across five developmental stages. Based on geNorm pairwise variation analysis in total samples (V3/V4 value below 0.15), we selected single and combinations of two or three reference genes to normalize TruFBP24 gene expression suggested by RefFinder. As shown in Figure 6, when using  the combinations of UFD1/ADF3, UFD1/HIS3, ADF3/HIS3, and ADF3/HIS3/UFD1, as reference genes for normalization, similar expression profiles of TruFBP24 were detected in both staminate and perfect flowers at each developmental stage. When the least stable gene GAPDH was used for normalization, TruFBP24 was 7 times higher expression in perfect flower than that in staminate flower at stage 1. Moreover, relative gene expression of TruFBP24 (20.07 ± 4.12) at stage 5 was remarkably lower FIGURE 6 | The expression level of the TruFBP24 gene in staminate and perfect flowers across five developmental stages by using different reference genes for normalization. The single or combinations of UFD1, ADF3, and HIS3, were used for normalization as the most stable reference genes, GAPDH was used as worst stable reference gene. than those obtained by normalizing with the suitable reference genes, indicating the low expression stability of GAPDH as reference gene for normalization qRT-PCR in all floral samples. Overall, our analysis suggested that the target gene's expression profiles are strongly affected by the choice of the reference genes.

DISCUSSION
Andromonoecious plant species, which produces both bisexual flowers and unisexual male flowers within the same individual, provides an excellent model system for studying establishment and development of bi-and unisexual flowers under a uniform genetic background. Using qRT-PCR technique, investigation of gene expression profiles in staminate and perfect flowers not only better understands the underlying molecular mechanism in regulation of flower formation, but also provides insights into the complex regulatory networks in flower sex differentiation. To obtain reliable and accurate quantification results in qRT-PCR analysis, it is fundamental to select and validate of appropriate reference genes for normalizing qRT-PCR data. In this study, ten traditional housekeeping genes and six newly candidate reference genes were selected for evaluation of expression stability based on T. rupestris transcriptome data and previous studies (Du et al., 2008;Li et al., 2017). To the best of our knowledge, this is the first report on the identification and validation of suitable reference genes for normalizing qRT-PCR analysis in andromonoecious plant species.
Ideally, an accurate reference gene should be expressed stably regardless of varied organ, tissue type or developmental stages. By using geNorm and NormFinder, and BestKeeper software, the 16 reference gene candidates showed various performances in expression stability. However, the ranked orders generated by three algorithms were not completely identical. The results of geNorm and NormFinder were similar in some cases, but they showed quite differences from the results obtained from Bestkeeper. For example, SPF was ranked first by BestKeeper in perfect flowers, while it was ranked at a medium position by geNorm and NormFinder. The results were similar to many previous studies, and this apparent inconsistency was probably due to the different principles in the three statistical algorithms (Xiao et al., 2015;Qi et al., 2016). Considered as an integrative statistical program, RefFinder has been widely applied to evaluate the overall stability of reference gene expression and determine appropriate reference genes for diverse plant species (Kim et al., 2015;Qi et al., 2016). In order to obtain an integral assessment of the most suitable reference gene for Taihangia, a comprehensive tool RefFinder was used to determine the final overall ranking by a comparison of different algorithms (Xie et al., 2012). Based on comprehensive RefFinder analysis, ADF3, HIS3, and UFD1 were identified as the most appropriate reference genes for normalization in Taihangia floral tissues and flowers. As a whole, the integrated results obtained from different programs, could lead to better accuracy for each gene in this study, suggesting that more than two algorithms should be used for stability evaluation of reference gene.
Increasing evidences indicated that the expressions of classical housekeeping genes do not always express stably not only across developmental stages but also in variable biological samples Niu et al., 2017). In this study, the expression stability of commonly used housekeeping genes such as EF-1α, GAPDH, ACT, UBQ, TUA, and HIS3, displayed remarkably different expressed patterns in Taihangia floral tissues and flowers. Of these selected housekeeping genes, HIS3 was identified as one of the most stable reference genes for staminate flowers, perfect flowers, and different floral tissues. As a major component of chromatin, HIS3 is thought to be crucial for maintenance of a stable chromatin structure (Oliver et al., 2013). In T. rupestris, they may remain expressed constitutively, and showed minimal changes in mRNA transcription in staminate and perfect flowers across developmental stages. For majority of common used housekeeping genes, including ACT, EF-1α, and UBQ, showed relatively stable expression with middle ranking orders calculated by RefFinder analysis. For example, the ACT gene, used as reference gene in previous study, is not the best choice for normalization in T. rupestris. A possible explanation was that that reference genes are implicated in multiple cell processes. Although TUA and GAPDH were also identified as stable reference genes for flower development in other plant species Karuppaiya et al., 2017), their expressions were the least stable in Taihangia flowers. In addition, the transcriptional profiles of target genes showed a strong bias via qRT-PCR validation, when these genes were used as reference gene for normalization gene expression in staminate and perfect flowers of T. rupestris. Similar to our analysis, a number of reports have demonstrated that the expressions of several traditional housekeeping genes, such as GAPDH and TUA, varied across flower developmental stage and were unsuitable for normalizing qRT-PCR data as internal control gene (Jin et al., 2013;Yuan et al., 2014). These results reported here supported that the traditional reference transcripts may not always show constitutive expression as has often been assumed, and thus it is prerequisite for verification of their expression stability before qRT-PCR normalization in specific biological samples.
Compare to traditional reference genes, the newly discovered reference genes could perform better in qRT-PCR normalization for non-model plants (Narsai et al., 2010;Gonzalezaguero et al., 2013). In this study, UFD1 showed the most stable expression profile in staminate flowers, while ADF3 was identified as the second stable gene in both staminate and perfect flowers. ADF are small actin-binding proteins, which involved in regulation of actin dynamics, play an essential role in plant growth and development . In a previous study, ADF was identified as a reference gene candidate for grapevine flower development based on transcriptome data because of stable expression during anthesis (Gonzalezaguero et al., 2013).In rubber tree, ADF was also identified as a candidate reference gene because its transcription abundance maintained relatively stable duration of latex flow (Chao et al., 2016). For andromonoecious T. rupestris, the stability of ADF3 in staminate and perfect flowers probably attributed to its function for basic cellular processes, such as cell proliferation and regulation of actin dynamics in cells . In plants, UFD1 is an essential ubiquitin recognition component in the ubiquitin-mediated degradation pathway (Wei et al., 2009). Given that possible basic function in cell processes, UFD1 should be suitable as reference gene for qRT-PCR normalization in T. rupestris, although it was not reported as reference gene in plant species. It was note worth that ADF3, HIS3, and UFD1 showed high stability among staminate and perfect flowers when total samples were taken into consideration, implying that those genes could provide accurate RT-qPCR analysis in comparison gene expression between uni-and bisexual flowers within andromonoecious system.
It is widely acceptable that the normalization with combination of multiple reference genes could provide more accurate relative expression quantification than that with a single gene in qRT-PCR analysis (Reid et al., 2006;Expositorodriguez et al., 2008). Based on validation of target gene expression, UFD1 combined with ADF3 for staminate flowers, the combination of ADF3/HIS3 for perfect flowers, HIS3/UFD1/TMP50 for floral tissues, and UFD1/HIS3 for total samples were recommended as appropriate reference genes for qRT-PCR normalization. As the three stable reference genes ADF3, HIS3, and UFD1 belong to the different functional classes, they should be used together as suggested by Vandesompele et al. (2002). Therefore, our results supported that the more accurate quantification of gene expression could be obtain when multiple reference genes were applied in qRT-PCR analysis. Interestingly, the combinations, including one traditional housekeeping gene and one new reference gene, were identified as the most stable reference genes for majority of group tested, suggesting that the selection of internal reference genes with traditional house-keeping gene combined novel reference gene may be a good strategy for normalizing qRT-PCR on T. rupestris flowers.
In summary, 16 reference gene candidates were selected based on our transcriptome sequence data and previous reports, and their expression stability were assessed in staminate and perfect flowers within andromonoecious plant T. rupestris. ADF3 combined with UFD1 were identified as the optimal reference genes for staminate flowers, while combination of HIS3/ADF3 was recommended to be the best one in perfect flowers, respectively. For floral tissues, combination of HIS3/UFD1/TMP50 was the most suitable reference genes for normalization. The stable reference genes reported in this study will be helpful to improve the accuracy of qRT-PCR analysis in andromonoecious T. rupestris, and will facilitate the functional genomics studies on flower development and sex differentiation in the future.

AUTHOR CONTRIBUTIONS
WL and YWZ conceived and designed the work. LZ, YDZ, and GW performed experiments and analyzed data. WL wrote the paper, DS partially revised the manuscript. All authors read and approved the manuscript.