Selection of Reference Genes for Quantitative Real-Time PCR during Flower Development in Tree Peony (Paeonia suffruticosa Andr.)

Tree peony (Paeonia suffruticosa) is a perennial plant indigenous to China known for its elegant and vibrantly colorful flowers. A few genes involved in petal pigmentation have been cloned in tree peony. However, to date, there have been few studies on the comparison and selection of stable reference genes for gene expression analysis by quantitative reverse-transcription PCR (qRT-PCR) in this species. In this study, 10 candidate reference genes were evaluated for the normalization of qRT-PCR in three tree peony cultivars. GAPDH and UBC were identified as the top two most stable reference genes in ‘Feng Dan’ and ‘Xi Shi,’ and EF-1α/UBC was recommended to be the best combination for ‘Que Hao.’ The expression stability of various reference genes differed across cultivars, suggesting that selection and validation of reliable reference genes for quantitative gene expression analysis was necessary not only for different species but also for different cultivars. The results provided a list of reference genes for further study on gene expression in P. suffruticosa. However, in any case, a preliminary check on the accuracy of the best performing reference genes is requested for each qRT-PCR experiment.


INTRODUCTION
Gene expression analysis provides improved understanding and insights into the molecular basis underpinning various biological processes . Quantitative real-time polymerase chain reaction (qRT-PCR) is one frequently used platform to quantify transcript abundance, for its sensitivity, accuracy and reproducibility (Gachon et al., 2004). However, the accuracy of qRT-PCR is heavily dependent on the stability of the internal reference genes used to normalize transcript abundances (Huggett et al., 2005). An ideal reference gene must be expressed steadily both in different tissues and at different developmental stages and remain unaffected by any experimental treatment (Czechowski et al., 2005). Many studies have confirmed the fact that the utilization of unstable reference genes may result in significant biases and misinterpretations of data (Ferguson et al., 2010;Mafra et al., 2012). It is thus critical to select suitable reference genes before quantifying the expression level via qRT-PCR.
The tree peony (Paeonia suffruticosa Andr.) is an ornamental plant native to China with striking ornamental and medicinal value. It has a cultivation history going back more than 1600 years in China (Cheng and Li, 1998). Qualitative and quantitative analyses of anthocyanins have been conducted in P. suffruticosa (Zhang et al., 2007;Yin et al., 2012). A few genes involved in petal pigmentation have also been cloned in the tree peony (Zhou et al., 2013;Du et al., 2015;Shi et al., 2015). However, to date, there have been few studies on the comparison and selection of reference genes for qRT-PCR during flower development in tree peony. Only few studies focusing on the bud dormancy release (Zhang et al., 2011), various tissues (Liu et al., 2015), and different treatments in cut flowers (Wang et al., 2012) have been reported, showing that no gene has a constant expression profile under all developmental or experimental conditions. To normalize the gene expression level during flower development in different cultivars of tree peony, 10 candidate reference genes were selected based on the floral transcriptome datasets of tree peony and the stability of their expression level was evaluated. The results provided a list of reference genes suitable for the accurate quantification of gene expression during flower development.

Plant Materials
Tree peonies were grown in the germplasms nursery of Shanghai Chenshan Botanical Garden. Three cultivars ('Feng Dan, ' 'Xi Shi, ' and 'Que Hao')

RNA Extraction and cDNA Synthesis
Total RNA was extracted using an E.Z.N.A. TM Plant RNA Kit-R6827 (Omega) following its manufacturer instructions. To avoid DNA contamination, the RNA samples were treated with an RNase-free DNase Set (Omega). The total RNA concentration, purity and integrity were determined using a NanoDrop2000c Spectrophotometer (Thermo Scientific, U.S.) and visually assessed via 1.5% agarose gel electrophoresis. The RNA samples with absorption ratios of A260/A280 between 1.9 and 2.1 and A260/A230 ratios above 2.0 were used for subsequent analyses. First-strand cDNA was synthesized using 1 µg total RNA with the PrimeScript RT Master Mix (Perfect Real Time; Takara, Japan) according to the manufacturer's protocol.

Candidate Reference Genes Selection and PCR Primer Design
An Illumina/Solexa library of tree peony was constructed in our preliminary study. A total of 15019 unigenes (NR database) were annotated with BLASTX (data was not shown). The indices of mean expression value (MV), standard deviation (SD), and coefficient of variation (CV) value of RPKM, were calculated using the method described by de Jonge et al. (2007) for mining candidate reference genes from transcriptome data. A total of 10 candidate genes were selected and assessed (Table 1).
Primers were designed with Primer Premier 5.0 (Lalitha, 2004), using the following criteria: targeting the 3 ′ -untranslated region, no hairpin/dimer/false priming/cross dimer structures, Tm = 55-60 • C, 18-25 bases in length, GC content between 40 and 60%, and amplicon length from 100 to 300 bp. The designed primer sets were BLASTed against the local floral transcriptional data of tree peony to verify primer specificity. Amplification efficiency (E) was evaluated using a standard curve generated by qRT-PCR using a 5-fold dilution series. The primer specificity was judged by melting-curve analysis and agarose gel electrophoresis analysis. at 95 • C and 30 s at 60 • C, and a melting curve protocol (65-95 • C with fluorescence measured every 0.5 • C). Each reaction was performed in three biological replicates and three technical replicates. Three no-template controls (NTC) were used to monitor possible DNA contamination.

Analysis of Gene Expression Stability and Validation of Reference Genes
Three software programs (GeNorm 3.5, BestKeeper and NormFinder) were used to rank the stability of the 10 selected reference genes in the study. Moreover, RefFinder (http://fulxie. 0fees.us/?type=reference), a comprehensive web-based tool was used to integrate and rank comprehensively the tested candidate genes. The raw Ct values were converted into relative quantities using the formula: E − Ct ( Ct= Ct value of each sample-the minimum Ct value) before data entry (Ramakers et al., 2003). All three of the software programs were run based on the software FIGURE 2 | Pairwise variation (V) of 10 candidate reference genes calculated by GeNorm. The cutoff value to determine the optional number of reference genes for normalization is 0.15. FD, 'Feng Dan'; XS, 'Xi Shi'; QH, 'Que Hao.' manuals to select suitable reference genes. RefFinder generated the final overall ranking of tested reference genes based on the geometric mean of the weights of every gene calculating by each program.
To validate the reliability of the reference genes, the relative expression levels of PsF3H in three cultivars, a key gene involved in the anthocyanin biosynthesis, was analyzed using the two topranking reference genes and the least stable reference gene, as determined by RefFinder, alone or in combination (calculated by geometric mean) for data normalization. The qRT-PCR amplification conditions were the same as described above. The primer specificity of PsF3H was verified as described for reference genes ( Table 1). The relative expression data was calculated according to the 2 − Cq method and presented as fold change (Livak and Schmittgen, 2001).

Amplification Specificity and Efficiency
The primer sequences and amplicon characteristics of 10 candidate reference genes were summarized in Table 1. The primer specificities were assessed by melting-curves and agarose gel electrophoresis analysis, which showed only a single product of the expected size ( Figure S1A) and a single peak melting curves ( Figure S1B) for each primer set. The amplification efficiency (E) of each primer pair varied between 90.51% for β-TUB and 104.74% for SAMS, and the regression coefficient (R 2 ) ranged from 0.991(SAMS) to 0.999 (α-TUB; Table 1). The amplicon size ranged from 100 to 264 bp.
To assess the expression stability of the reference genes in different samples, the transcript abundances of 10 candidate reference genes were presented as their mean Ct values (quantification cycles). The mean Ct values of these candidates varied from 20.9 ± 1.2 (SAMS) to 26.8 ± 1.0 (ACTIN) in different samples. Genes with greater SD of Ct values have more variable expression than those with lower SD. Of the 10 candidate genes, UBC showed the least variation in gene expression (23.7 ± 0.6), and the expression of SAMS was the most variable (20.9 ± 1.2). The candidate genes can be ranked as: UBC > EF-1α/TIP41 > PP2A/CYP > GAPDH > α-TUB/ ACTIN > β-TUB > SAMS (Table S1) based on the standard deviation of the Ct values.

Expression Stability of the Reference Genes
Different software programs (GeNorm, NormFinder, BestKeeper, and RefFinder) were used to assess and rank the expression stability of the candidate genes. The average expression stability value (M-value) and the average pairwise variation of template normalization factor (V n/n+1 value) were used to identify the most stable reference genes and obtain the required reference gene number in GeNorm. The M-value of a gene was inversely correlated with its expression stability, with the M-values less than 1.5 indicating stable expression (Hellemans et al., 2007). The cut-off value of 0.15 was used as a threshold of the V n/n+1 value. If V n/n+1 ≤ 0. 15, it is not necessary to introduce n+1 reference genes as the internal control. Figure 1 showed the ranking order of different candidate genes according to the M-value. For total samples, PP2A and Tip41 were identified as the most stable genes, while β-TUB and SAMS were the least stable reference genes. As depicted in Figure 2, the optional number of reference genes for qRT-PCR data normalization was two in each cultivar. Taking the results of the most stable gene selection into account, the combinations GAPDH/EF-1α, UBC/GAPDH, and UBC/EF-1α were predicted to deliver the most reliable level of normalization for 'Feng Dan, ' 'Xi Shi, ' and 'Que Hao, ' respectively. NormFinder was used to determine the stability of reference genes based on inter-and intra-group variance in expression level. The stability value (SV) was calculated, with the lower stability value indicating the higher stability (Andersen et al., 2004). The results of NormFinder analysis were shown in Table 2. For total samples, UBC and TIP41 were the most stable genes, while SAMS and β-TUB being the least stable. The top two most stably expressed genes were GAPDH and UBC in 'Feng Dan, ' GAPDH and TIP41 in 'Xi Shi,' respectively. The ranking order was slightly different from that of GeNorm.
Coefficient of variance (CV) and the standard deviation (SD) were calculated by Bestkeeper analysis. The gene with the lowest CV and SD was considered to be the most stable reference gene (Chang et al., 2012). As shown in Table 3, PP2A and α-TUB were suggested to be the most stable genes in 'Feng Dan.' UBC and GAPDH had the lowest SD values in 'Xi Shi, ' and CYP and UBC performed the best in 'Que Hao.'

Comprehensive Ranking of the Reference Genes by RefFinder
RefFinder was used to calculate the geometric mean of weights for the comprehensive ranking order recommended by geNorm, NormFinder, and BestKeeper. The comprehensive ranking order recommended by RefFinder was shown in Table 4. In total samples, the ranking order recommended by RefFinder was similar to those of NormFinder and GeNorm: TIP41> UBC> GAPDH> PP2A> EF-1α > CYP> ACTIN> α-TUB> β-TUB >SAMS (from the most stable to the least stable). GAPDH and UBC were ranked as the top two stable genes in 'Feng Dan' and 'Xi Shi.' EF-1α and UBC were the most stable genes in 'Que Hao.' β-TUB and SAMS were supposed to be the least stable genes in three cultivars.

Validation of the Reference Genes
To assess the reliability of the selected candidate reference genes, the transcription levels of PsF3H were evaluated using different reference genes in three cultivars. The two most stable reference genes identified in each cultivar (UBC and GAPDH for 'Feng Dan' and 'Xi Shi, ' EF-1α and UBC for 'Que Hao') and the least stable reference gene (SAMS or β-TUB) were used either singly or in combination in qRT-PCR analyses.
Although the overall expression patterns of PsF3H were similar, differences were found when normalized to different reference genes. In 'Feng Dan' and 'Xi Shi, ' the normalized expression levels of PsF3H using the least stable gene SAMS were dramatically decreased compared to the stable genes GAPDH and UBC (Figures 3A,B). Conversely, in 'Feng Dan' , 'Xi Shi' and 'Que Hao, ' the expression level normalized by β-TUB   'Feng Dan'; XS, 'Xi Shi'; QH, 'Que Hao.' significantly increased compared to EF-1α and UBC ( Figure 3C). When the combinations of GAPDH/SMAS and UBC/SAMS were used as internal control, the expression patterns of PsF3H were similar to that of SAMS (Figures 3A,B). When the combinations of EF-1α/β-TUB and UBC/β-TUB were used, the expression level of PsF3H at stage 6 was significantly overestimated. The combination of GAPDH/UBC was recommended to be the optimum pairs of reference gene for 'Feng Dan' and 'Xi Shi, ' and EF-1α/UBC was the suitable pairs of reference gene for accurate normalization in 'Que Hao.'

DISCUSSION
To ensure the accuracy of analysis, it is important to validate suitable reference genes for gene expression analysis using qRT-PCR. In this study, 10 reference genes were evaluated for the normalization of qRT-PCR in tree peony cultivars. For the cultivars 'Feng Dan' and 'Xi Shi, ' the top two most stable reference genes, GAPDH and UBC, were similar, while EF-1α was identified as the most stable gene in 'Que Hao.' The expression stability of various reference genes differed across cultivars. Similar results were reported in strawberries and Panax ginseng (Galli et al., 2015;Wang and Lu, 2016). The diverse genetic background and biological process between cultivars probably affected the expression stability of reference genes. The results suggested that selection and validation of reliable reference genes for quantitative gene expression analysis was necessary not only for different species but also for different cultivars. Among the candidate genes evaluated, UBC and GAPDH were the two most stable reference genes in both 'Feng Dan' and 'Xi Shi' cultivars. This result was consistent with previous studies in Salicornia europaea (Xiao et al., 2015), Plukenetia volubilis (Niu et al., 2015), and navel oranges . They seemed unstable for normalization at different development stages in petunia hybrid (Mallona et al., 2010). Both ACTIN and β-TUB were involved in basic cellular processes, and have long been considered suitable reference genes in numerous species (Jin et al., 2013;Zhang et al., 2015). However, in this study, β-TUB was identified as the least stable genes in all three cultivars.
With the development of RNA sequencing, increasing amounts of transcriptome data have been successfully utilized to screen reference genes for non-model plants (Demidenko et al., 2011;Xiao et al., 2015;Zhuang et al., 2015). It has been demonstrated that the newly discovered reference genes could perform better than traditional reference genes in Brassica napus  and tea (Hao et al., 2014). The results of this study also showed that, compared to classic reference genes, the newly discovered ones were more stably expressed in tree peony cultivars.
In conclusion, this study evaluated the expression stability of 10 candidate reference genes in three tree peony cultivars. GAPDH/UBC was identified as the optimum pair of reference genes in 'Feng Dan' and 'Xi Shi.' EF-1α/UBC was recommended to be the best combination for 'Que Hao.' The results provided useful information for selection and validation of reference genes for transcript normalization in gene expression studies in P. suffruticosa. In any case, a preliminary check on the accuracy of the best performing reference genes is requested for each qRT-PCR experiment.