Evaluation of reference genes for normalizing RT-qPCR in leaves and suspension cells of Cephalotaxus hainanensis under various stimuli

Background Reverse transcription quantitative real-time PCR (RT-qPCR) is a widely used approach for investigating gene expression levels in plants because of its high reproducibility, sensitivity, accuracy and rapidness. Evaluation of reference genes for normalizing RT-qPCR data is a necessary step, especially in new plant varieties. Cephalotaxus hainanensis is a precious medicinal plant belonging to the family of Cephalotaxaceae and no RT-qPCR studies have been reported on it. Results In this study, 9 candidate reference genes were selected from the transcriptome data of C. hainanensis; 3 statistical algorithms (geNorm, NormFinder, BestKeeper) were applied to evaluate their expression stabilities through 180 samples under 6 stimuli treatments in leaves and leaf-derived suspension cultured cells; a comprehensive stabilities ranking was also performed by RefFinder. The results showed that suitable reference genes in C. hainanensis should be selected for normalization relative to different experimental sets. 18S showed a higher stability than other candidate reference genes which ranked at the top two suitable genes under all experimental setups in this study. Conclusion This study is the first to evaluate the stability of reference genes in C. hainanensis and supply an important foundation to use the RT-qPCR for an accurate and far-reaching gene expression analysis in C. hainanensis. Electronic supplementary material The online version of this article (10.1186/s13007-019-0415-y) contains supplementary material, which is available to authorized users.

the main research direction and studies on functional genes in their biosynthesis pathway has become the principal issue.
Nowadays, substantial researches on medicinal plants aim to increase their active compounds through inducing and regulating the functional genes in their biosynthesis pathways [8]. Many reports indicated that the synthesis and accumulation of plant secondary metabolites have significant correlation with the expression level of functional genes in their biosynthesis pathways [9,10]. In these reports, analysis of functional genes by reverse transcription quantitative real-time PCR (RT-qPCR) has become a widely used approach because of its high reproducibility, sensitivity, accuracy and rapidness [11][12][13]. However, the results of RT-qPCR must to be normalized by reference genes and the stability of reference genes are indeed the critical factor for a reliable result [14,15]. Many reports indicate that ideal stable reference genes do not exist and the stability of reference genes are varied in different plant species, growth environment, growth stage and stress conditions [16,17]. Therefore, systematic experimental design for the selection of stable reference genes is the first step towards applying RT-qPCR on a new species. Thus far, stable reference genes have been identified and verified in many plant species including rice [18], potato [19], rhododendrons [20], Taihangia flower [21], sugarcane [22], grapevine [23], radish [24], lettuce [25], Euscaphis konishii Hayata [26], Baphicacanthus cusia (Nees) Bremek [27], creeping bentgrass [28], Setaria viridis [29], etc. The reference genes used for these species include ACT, NAC, UBQ, F-box, PP2C, TUA, TUB, UBC, 18S, and so on. However, the suitable reference genes vary among these plants and even among treatments in same species. At present, the selection and evaluation of reference genes in Cephalotaxaceae remains unreported. It is necessary to perform a multifactorial analysis to identify the stability of reference genes which will greatly facilitate mining the functional genes in Cephalotaxus ester alkaloids biosynthesis pathway by comparative transcriptomics.
In this study, nine candidate reference genes (ACT, NAC, UBQ, F-box, PP2C, TUA, TUB, UBC, 18S) of C. hainanensis were selected based on the transcriptome sequencing data by the SMRT (Single-Molecule Real-Time) technology on PacBio Sequel (Unpublished). Six treatments (ABA, Ethylene, Mannitol, MeJA, NaCl, SA) were set to identify their expression level by RT-qPCR in leaves and suspension cells of C. hainanensis. The cycle threshold (Ct) values which were detected by RT-qPCR, indicated the expression levels of candidate reference gene directly, lower Ct values represented higher expression levels. Three statistical algorithms (geNorm, NormFinder, BestKeeper) were applied to evaluate their expression stability for normalization. Moreover, comprehensive stability ranking was also performed by Ref-Finder. Norcoclaurine synthase gene (NCS) catalyzes the first step in the biosynthesis of a diverse class of benzylisoquinoline alkaloids (BIAs) and may also be involved in the biosynthesis of Cephalotaxus ester alkaloids in C. hainanensis. Therefore, ChNCS was selected as the target gene, and its expression was used to verify the reliability of the selected reference genes. Different algorithms were used to screen most stable reference genes and a reliable set of reference genes were provided. Moreover, the criteria for a good reference gene was also discussed.
In this study, compound leaves excised from five-yearold C. hainanensis were used for stress treatments. The leaves were immersed into MeJA (100 μM), SA (100 μM), Mannitol (100 μM) NaCl (100 μM), ABA (10 μM) and ethylene (250 μM) respectively for about 1 s and then placed in a moisture chamber. For suspension cultured cells, the corresponding solutions were added into the medium directly to the indicated concentration. After treatment, leaves and suspension cells were harvested at time point 0 h, 2 h, 6 h, 12 h and 24 h and frozen in liquid nitrogen immediately then stored at −80 °C for RNA extraction. Three biological replicates were prepared for each treatment.

RNA isolation and cDNA synthesis
Total RNA was isolated from all prepared samples using RNAprep Pure Kit (Polysaccharides & Polyphenolicsrich) (Tiangen, Beijing, China) and genomic DNA was removed with RNase-free DNase I according to the manufacturer's instructions. The RNA concentration and purity were determined using a Nano Photometer P-Class instrument (Implen, Munich, Germany), and the RNA integrity was also checked on 1% agarose gels. Total RNA (1.0 μg) was used for reverse transcription with a FastQuant RT Kit (Tiangen, Beijing, China) in a 20 μL reaction volume according to the manufacturer's instructions.

Primer design and RT-qPCR conditions
Sequences of candidate reference genes (Additional file 5) were mined from our full-length transcriptome database obtained by the SMRT (Single-Molecule Real-Time) technology on PacBio Sequel platform (Unpublished, Novogene, Beijing, China). Specific primer pairs were designed using Beacon Designer 8 software according to primer sequences of 18-24 nucleotides, amplicon length of 75-150 bp, melting temperature (Tm) of 55-60 °C and GC content of 40-60%. All primer pairs were synthesized by a commercial supplier (Sangon, Shanghai, China) and tested by regular PCR and the products were analyzed by electrophoresis on 1.0% agarose gels before RT-qPCR. In addition, amplification efficiency (E) and correlation coefficients (R 2 ) were calculated by a standard curve with a series of 5 different cDNA dilutions. The primer sequences, amplicon length, Tm, GC content, amplification efficiency and correlation coefficients of nine candidate reference genes are listed in Table 1.
RT-qPCR was carried out in 384-well plates with a QuantStudio 6 Flex real-time PCR system (ThermoFisher, MA, USA) using SYBR Green-based PCR assay. The final reaction volume for each reaction was 10 μL with the following components: 1 μL diluted cDNA template (1 μg), 5 μL SYBR Premix Ex TaqII (TAKATA, Dalian, China), 1.2 μL forward primer (2.5 μM), 1.2 μL reverse primer (2.5 μM), 1.6 μL ddH 2 O. The reaction was conducted under the following conditions: 95 °C for 7 min, followed by 40 cycles of denaturation at 95 °C for 10 s, and annealing/extension at 56 °C for 30 s. The melting curve was obtained by heating the amplicon from 65 °C to 95 °C with increasing 1.0 °C/s. Each RT-qPCR analysis was performed with three technical replicates.

Data analysis of gene expression stability
Four statistical tools, geNorm, NormFinder, BestKeeper and RefFinder were used to analyze the candidate reference gene's stability based on their own algorithms. For geNorm, the expression stability value (M) of each reference gene was calculated based on the average pairwise variation (V) between all genes tested [30]. For Nor-mFinder, an ANOVA-based model of each reference gene was used to calculate the expression stability value by determining inter-and intra-group variation, the gene with the lowest value has the most stable expression [31]. For BestKeeper, the standard deviation (SD) and coefficient of variance (CV) were used to calculate the expression stability of candidate reference genes with raw Ct values, the lowest CV representing the highest stability [32]. For RefFinder (http://150.216.56.64/refer enceg ene. php), a comprehensive ranking was generated with the data from GeNorm (M-values), NormFinder (stability values) and BestKeeper (SD and CV).

Validation of reference gene stability
Norcoclaurine synthase (NCS) is the first committed enzyme and catalyses a central precursor in the biosynthesis pathway of alkaloids derived from tyrosine and phenylalanine [33,34]. In this study, we used the NCS gene in C. hainanensis (ChNCS) as target gene to confirm the reliability of the potential reference genes in RT-qPCR. The relative expression level of ChNCS under ethylene stress treatment was determined and normalized using the most and least stable reference genes according to the statistical software in the same RT-qPCR conditions mentioned above. The relative expression data was calculated by 2 −ΔΔCt method and three technical replicates were performed for each sample.

Primers verification and expression levels of candidate reference genes
Primer specificity amplification of all candidate reference genes were verified by regular PCR and RT-qPCR. All primers pairs amplified their specific amplicon based on agarose gel electrophoresis and a single peak with the melting curve analysis (Additional file 1). All candidate reference gene names and abbreviation, accession number, primer sequences, amplification length, efficiency (E) and correlation coefficient (R 2 ) are listed in Table 1.
Their amplification efficiency varied from 91.1% for TUB to 101.3% for F-box, and correlation coefficients ranged from 0.988 for ACT to 0.999 for PP2C and 18S. The expression levels of candidate reference genes were determined by Ct values directly and showed in Fig. 1.

Expression stability of candidate reference genes geNorm analysis
For geNorm analysis, the stability of all candidate reference genes was evaluated by M-values below the threshold of 1.5, which were calculated by the mean variation of a gene relative to all others. Lower M-value represented higher gene expression stability [30]. Based on the geNorm analysis, M-values were calculated for leaf and suspension cell samples subjected to different treatments respectively. The ranking of the reference genes was found to differ between the different experimental conditions. 18S and TUB had the lowest M-values and thus were most stable in most leaf treatments; 18S and TUA had the lowest M-values and thus were most stable in most suspension cells treatments. Contrarily, ACT and NAC had higher M-values in most leaf samples and PP2C, F-box and NAC had higher M-values in most suspension cell samples (Fig. 2). geNorm algorithms also can determine the optimal number of reference genes for normalization calculated by the pairwise variation Vn/n + 1. Ideal pairwise variation (V) score of below 0.15 was recommend. But in this study, most results of pairwise variation calculated by geNorm were more

NormFinder analysis
For NormFinder analysis, the stability of all candidate reference genes was evaluated by taking into account intragroup and inter-group variations. The stability values were used to rank candidate reference genes with lower values indicating more stability. As shown in Table 2, the most stable reference genes presented in most leaf samples were 18S and UBQ, the most stable reference genes presented in suspension cells samples were 18S, UBC and TUA . Relatively, the least stable reference genes were TUA , NAC and ACT in leaf samples and NAC, F-box and PP2C in suspension cells samples.

BestKeeper analysis
For BestKeeper analysis, the stabilities of candidate reference genes were evaluated using the SD and CV values which were calculated by raw Ct value data directly, lower SD and CV value represented higher gene expression stability, especially when the value SD > 1 indicated the reference gene was unstable and cannot be used for normalization. The results obtained with BestKeeper analysis are shown in Table 3 Table 3, any reference gene with value SD > 1 also presented lower stability and cannot be used for normalization.

RefFinder analysis
The RefFinder approach was used to determine the comprehensive rankings of candidate reference genes based on the results of common analysis programs (geNorm, NormFinder, BestKeeper). The comprehensive ranking calculated by RefFinder is shown in Table 4

Reference genes validation
The ethylene induced expression level of ChNCS in leaf and suspension cell samples was normalized to validate the selected reference genes. in C. hainanensis. According to the comprehensive analysis of geNorm, Nor-mFinder, BestKeeper and RefFinder, two sets of reference genes were selected. The most stable reference genes were 18S and UBQ for Leaf-Ethylene samples and were assumed as '1' and we used 2 ( −ΔΔCt ) to calculate its relative expression in samples at other time points. As shown in Fig. 3, results for Leaf-Ethylene samples showed that when the most stable reference genes (18S and UBQ) were used for normalization, the relative expression of ChNCS was 0.51 and 0.46 times higher for 2 h samples; 2.24 and 1.60 times higher for 6 h samples; 5.77 and 3.52 times higher for 12 h samples; and 0.11 and 0.09 times higher for 24 h samples. No large difference was observed in the expression trend of both 18S and UBQ normalization results. However, when the least stable reference genes TUB was used for normalization, the relative expression of ChNCS showed a different result especially for the 12 h samples, where it is only 1.39 times higher than the 0 h samples. Similar results were observed for the Cell-Ethylene samples, when the most stable reference genes were used for normalization, the relative expression of ChNCS was 5.68 and 4.95 times higher for 2 h samples. Conversely, a large difference was evident in the change patterns when the least stable reference gene was used for normalization, the relative expression of ChNCS was only 0.66 times higher for 2 h samples.

Discussion
The endangered C. hainanensis harbors unique genes related to medicinally active Cephalotaxus ester alkaloids. Quantifying these genes with stable reference genes will greatly facilitate to decipher the biosynthesis pathway. Although leaves have a comparatively lower concentration of interesting metabolites than bark, the fact that it is more renewable and can give a higher biomass makes it a potentially attractive source tissue from this slow growing tree [35]. Moreover, suspension culture cells have shown the potential to produce interesting secondary metabolic compounds by plant cell fermentation techniques [36] and genetic engineering [37]. Therefore, leaves and suspension cultured cells were selected as materials to evaluate the stabilities of reference genes in this study. The stability of reference genes was analyzed comprehensively in C. hainanensis under NaCl, mannitol stresses and under challenge of stress-related signaling molecules ABA, ethylene, MeJA and SA. NaCl and mannitol, which have been widely observed to play key roles in the reponse to salt and osmotic stress, can be applied exogenously to mimic such stresses [38].
Previously published studies have pointed out that the stable reference genes varied among species and even changed under different experimental treatments in same species [17]. Since there is still no report on the stability of reference genes in C. hainanensis, we selected 9 candidate reference genes for evaluation. By consulting the public studies in other species as well as our transcriptome database, 9 traditional reference genes were extracted as candidates which include protein synthesis (18S), cytoskeleton structure (ACT , TUA , TUB), biological metabolic processes (UBQ, UBC) and novel reference genes of protein phosphatase (PP2C), F-box family protein (Fbox), domain-containing protein (NAC). All these 9 candidate reference genes had been identified and evaluated in other plant species. In Baphicacanthus cusia, 18S was found to be the most stable gene under ultraviolet irradiation and hormonal stimuli (MeJA and ABA) and UBC was the best suitable gene for different plant organs [27]. In Lycoris aurea, the comprehensive ranking results presented that UBC was the most stable reference gene when challenged by NaCl and cold stress subsets [39]. In Panax ginseng, ACT was one of top three-ranked genes in seedlings treated with heat [40]. In Rhododendrons, ACT and 18S were found to be the top choices for different tissues, whereas TUB was not found to favor RT-qPCR normalization in these tissues and NAC also was not the suitable reference gene for different tissues [20]. In Gentiana macrophylla, F-box, UBQ and UBC were tested as reference genes, but only UBC could be considered as reference gene as only the expression of UBC was stable enough in multiple tissues and environments [41].
Since the raw Ct values of candidate reference genes are also the direct readout of expression level [16,28], we then firstly evaluated stability of Ct values in leaf and suspension cell samples. However, the raw Ct values showed large differences under same treatment samples between leaf and suspension cells, which can be seen obviously in Fig. 1, there are approximately 5-7 Ct values in leaf samples more than the same treatment in suspension cell samples. Hence, the leaf and suspension cell samples were separated in the following analysis by geNorm, NormFinder, BestKeeper and RefFinder, which are the most popular statistical algorithms and widely used in recent studies. Numerous reports had confirmed that the ranking results obtained by different statistical algorithms were not completely identical since their different calculation methods [20,21,25,27,41]. In this study, the similar conclusions were obtained after analysis by statistical software. For example, the top three ranking results under Leaf-Ethylene were UBQ, F-box and 18S by geNorm algorithms; 18S, Fbox and UBQ by NormFinder algorithms; 18S, TUA and ACT by BestKeeper algorithms. Relatively, the three least stable reference genes were also different under Leaf-Ethylene, PP2C, ACT and TUB by geNorm algorithms; PP2C, ACT and NAC by NormFinder algorithms; NAC, PP2C and F-box by BestKeeper algorithms. In order to obtain a relatively objective result, the RefFinder software was used for comprehensive ranking and based on this the 18S was identified as the most stable reference gene for Cephalotaxus hainanensis. Furthermore, the same analysis results were presented under other leaf and suspension cell    treatment sets. In addition, as shown in Table 4, Additional files 3 and 4, the comprehensive ranking demonstrates that 18S had a higher ranking number and more stable expression level among all the samples. Many reports indicated that multiple reference genes used for normalization would obtain more accurate results by RT-qPCR and the optimal number could be calculated by geNorm algorithms with ideal pairwise variation (V) score below cut-off value of 0.15 [25,27]. As shown in Additional file 2, most treatment sets under leaf and suspension cells didn't result the ideal value (0.15) and most of their V scores were higher than 0.15. However, the M-values, which can also be calculated by geNorm algorithms to evaluate the stability of single reference gene, were all below the threshold of 1.5. Though the V scores in this study might intuitively present the pairwise variation Vn/Vn + 1 value, the alternative M-values indicated that the cut-off value 0.15 of V scores must not be considered as the only criterion. With M-values, it is even possible to obtain accurate results by single reference gene. The golden rules of qRT-PCR have suggested to take at least 4 reference gene to comprise the deviation by single reference gene [15], here, we provide the priority of selection of reference genes, and suggest that even one high ranking reference genes may be a better choice. We validated the selected reference genes with relative expression level of ChNCS under ethylene induction. The normalization results in Fig. 3 had shown obviously, the relative expression trend of ChNCS had small changes when using the most stable reference genes to normalize both in leaf and suspension cell treatment samples. By contrast, large changes were noted when the least stable reference genes were used for normalization, even a wrong trend were presented in 2 h suspension cell treatment samples. These validation results confirmed the applicability and correctness of the reference genes selected and evaluated in C. hainanensis, also indicated stable reference genes selection and evaluation represent a crucial issue for the proper normalization of the RT-qPCR data.

Conclusion
In conclusion, this study examined the selection and evaluation of 9 candidate reference genes for RT-qPCR normalization under 6 abiotic stresses treatments in C. hainanensis leaves and suspension cells. To the best of our knowledge, this work is the first to validate reference genes in C. hainanensis for the normalization of the RT-qPCR data. Based on the above results, we recommend 18S as a suitable reference gene for normalizing expression levels in C. hainanensis. Other reference genes can be selected as well but which are most suitable differs depending on the tissue and experimental treatment. These selected stable reference genes collectively supply an important foundation to use the RT-qPCR for an accurate and far-reaching gene expression analysis in C. hainanensis.

Additional files
Additional file 1. Gene specificity and amplicon size. Melting curves of 9 reference genes showing single peak.
Additional file 2. Pairwise variation (V) analysis of the 9 candidate reference genes.
Additional file 3. Expression stability ranking of the 9 candidate reference genes in Cephalotaxus hainanensis leaf samples.

Additional file 4.
Expression stability ranking of the 9 candidate reference genes in Cephalotaxus hainanensis suspension cells samples.
Additional file 5. The sequences of 10 candidate genes in Cephalotaxus hainanensis.