Screening and verification of real-time fluorescence quantitative PCR reference genes related to flowering regulation of tree peony

Paeonia suffruticosa Andrews. is a perennial deciduous subshrub belonging to the family Paeoniaceae, genus Paeonia L., section Moutan DC., and represents a plant resource endemic to China^[1]. The flowering period in tree peony is characteristically short and highly synchronized: under natural conditions, the period from bud break to petal senescence spans only 50–60 days; individual flowering typically last 3–5 days, while the population-level flowering window extends for approximately 10–15 days. Moreover, the majority of cultivars are classified as mid-season flowering types, with early- and late-flowering cultivars accounting for only a small proportion of the available germplasm^[2]. These phenological traits have a direct influence on both the ornamental appeal and economic value of tree peony^[3]. Meanwhile, flowering time represents a critical developmental transition in tree peony. The regulatory mechanism of flowering in this species is complex and highly sensitive to environmental cues such as temperature, photoperiod, and endogenous hormonal balance^[4]. Precise regulation of flowering time not only determines the synchronization of blooming and market timing but also influences breeding efficiency and genetic improvement strategies. Therefore, elucidating the molecular mechanisms underlying flowering-time regulation in tree peony is of significant importance for both fundamental research and horticultural applications^[5].

Gene expression analysis is a key method for understanding the molecular basis of biological processes and identifying functional genes. The common housekeeping genes used in such studies include Tubulin, Actin, 18S ribosomal RNA (18S rRNA), and glyceraldehyde-3-phosphate dehydrogenase (GAPDH)^[6,7]. However, accumulating evidence indicates that many of these traditional housekeeping genes exhibit considerable expression variability across different tissues, developmental stages, and experimental treatments, thereby undermining their reliability as reference genes for normalization. Consequently, the suitability of several widely used housekeeping genes as internal controls has been increasingly questioned^[8,9]. Therefore, to minimize the bias in gene expression analyses arising from multiple experimental variables, it is essential to select appropriate reference genes as internal controls for normalization^[10]. Real-time quantitative PCR (qRT-PCR) is widely employed for gene expression analysis due to its high sensitivity, specificity, and reproducibility. However, the reliability of qRT-PCR results strongly depends on the use of appropriate reference genes for data normalization^[11]. An ideal reference gene should exhibit stable expression across different tissues, developmental stages, and experimental conditions. In tree peony, gene expression patterns are highly dynamic during floral bud differentiation and flower development, and can vary considerably among cultivars. This complexity, combined with environmental influence, poses significant challenges for the identification of stable internal control genes. Furthermore, although RNA-Seq datasets provide a valuable resource for candidate gene selection, the criteria for screening and validating reference genes are insufficiently standardized, particularly in flowering-related studies of perennial ornamental species^[12]. When assessing the differential expression of a target gene, factors such as initial template concentration, reverse transcription efficiency, primer specificity, and amplification efficiency should ideally be consistent among samples^[13]. However, it is practically difficult to achieve simultaneous consistency in these conditions. Internal reference genes, which exhibit relatively stable expression across tissues within a species, are primarily used to correct unavoidable technical variations—particularly those introduced during sample handling and pipetting—thereby ensuring the reliability and reproducibility of experimental results^[14]. Consequently, the careful selection of suitable reference genes is indispensable for accurate normalization and reliable relative quantification of gene expression levels^[15]. Ideally, reference genes should maintain stable expression across different treatments, tissues, and even across species^[16].

Based on RNA-Seq data from early- and late-flowering tree peony varieties during floral development^[17], 23 candidate reference genes (EF1-α-1, UBi-lp, UBQ, UBC, UPL, GAPDH, HIS, CYP, RPL, CSLD2, KOR, AGD, ARP, IWS, GTP, PP2A, eIF2, eIF3, PPR, HDH, and PPC) were selected for evaluation. These candidates, together with two previously identified optimal reference genes (EF1-α and Tubulin-α) from drought-stress RNA-Seq studies in tree peony, were systematically re-assessed^[18]. Their expression stability was analyzed across flowering regulation and floral developmental stages using geNorm, NormFinder, BestKeeper, and RefFinder algorithms to identify the most suitable reference genes for normalization.

Petals at full bloom were collected from tree peony cultivars with uniform growth, free of pests and diseases. These cultivars represented distinct flowering times, including 21 early-flowering and 21 late-flowering cultivars (Fig. 1). Additionally, floral buds or petals of Paeonia ostii 'FengDan' and Paeonia suffruticosa 'Lian He' were collected at 16 different stages of floral bud differentiation and flower development (Fig. 2). Each developmental stage was sampled in triplicate. All samples were immediately flash-frozen in liquid nitrogen and stored at −80 °C for subsequent analyses.

Figure 1. 

Floral characteristics of 21 early-flowering and 21 late-flowering peony cultivars at full-bloom stage. A: Paeonia suffruticosa 'Taoyan Hong', B: Paeonia suffruticosa 'Huolian Jindan', C: Paeonia suffruticosa 'Pinghu Qiuyue', D: Paeonia suffruticosa 'JiaoHong', E: Paeonia suffruticosa 'Dapeng Zhanchi', F: Paeonia suffruticosa 'CangJiao', G: Paeonia suffruticosa 'LanJu', H: Paeonia suffruticosa 'HaiBo', I: Paeonia suffruticosa 'Lan Hudie', J: Paeonia suffruticosa 'Lan Yueliang', K: Paeonia suffruticosa 'ZhaoFen', L: Paeonia suffruticosa 'Yucui Caidie', M: Paeonia suffruticosa 'Jianshi Fen', N: Paeonia suffruticosa 'Xishi Fen', O: Paeonia suffruticosa 'Mantian Xing', P: Paeonia suffruticosa 'JingYu', Q: Paeonia suffruticosa 'Suxin Bai', R: Paeonia suffruticosa 'Xueyuan Hongxing', S: Paeonia rockii 'Ziban Bai', T: Paeonia suffruticosa 'ErQiao', U: Paeonia ostii 'FengDan', a: Paeonia suffruticosa 'Feiyan Hongzhuang', b: Paeonia suffruticosa 'Shouan Hong', c: Paeonia suffruticosa 'Haitang Zhengrun', d: Paeonia suffruticosa 'Daduo Lan', e: Paeonia suffruticosa 'Doulv', f: Paeonia suffruticosa 'Shenhei Zi', g: Paeonia suffruticosa 'Zihong Zhengyan', h: Paeonia rockii 'MingMou', i: Paeonia rockii 'Xian Emao', j: Paeonia rockii 'Bai Zhangbing', k: Paeonia rockii 'Bai Yanwei', l: Paeonia rockii 'Yuban Xiuqiu', m: Paeonia rockii 'Bingxin Fenhe', n: Paeonia suffruticosa 'Fen Zhouchou', o: Paeonia suffruticosa 'Xiuqiu Hong', p: Paeonia suffruticosa 'Zilou Xiangcui', q: Paeonia suffruticosa 'ZiYan', r: Paeonia suffruticosa 'JinGe', s: Paeonia suffruticosa 'JinZhi', t: Paeonia suffruticosa 'Baiwang Shizi', u: Paeonia suffruticosa 'Lian He'.

Figure 2. 

The flower development of Paeonia ostii 'FengDan' occurred in 16 stages: (a) overwintering dormancy, (b) bud break, (c) bud swelling, (d) bud elongation, (e) erect bud stage, (f) small bell stage, (g) large bell stage, (h) round peach stage, (i) flat peach stage, (j) color exposure, (k) bud cracking, (l) initial flowering, (m) half flowering, (n) full bloom, (o) onset of senescence, (p) full senescence. (a)–(e): bud stage, (f)–(i) flower bud stage, (j)–(p) flowering stage.

Total RNA extraction and cDNA synthesis

Total RNA was extracted from tree peony petals using a polysaccharide- and polyphenol-specific plant RNA extraction kit (TIANGEN, Beijing) according to the manufacturer's protocol. The integrity, concentration, and purity of the extracted RNA were assessed by performing 1% agarose gel electrophoresis and quantified using a NanoDrop™ One UV-Vis spectrophotometer. Qualified RNA samples were treated with 2 U of gDNA Eraser at 42 °C for 2 min to eliminate the residual genomic DNA contamination. Subsequently, reverse transcription was performed using the PrimeScript™ RT reagent kit with gDNA Eraser (Perfect Real Time; TaKaRa, Japan) to synthesize sufficient cDNA for downstream gene expression analyses.

Selection of candidate reference genes and primer design

Based on existing tree peony RNA-Seq data generated in our laboratory^[17], 23 candidate reference genes exhibiting both high expression levels and stable expression patterns were selected. During the screening process, the expression abundance and stability of candidate genes were systematically evaluated. The fragments per kilobase of transcript per million mapped reads (FPKM) value was used as a normalized measure of gene expression level, and genes with a mean FPKM < 1 were excluded. Based on this criterion, a total of 23 candidate reference genes were selected for subsequent validation by qRT-PCR: EF1-α-1, UBi-lp, UBQ, UBC, UPL, GAPDH, HIS, CYP, RPL, CSLD2, KOR, AGD, ARP, IWS, GTP, PP2A, eIF2, eIF3, PPR, HDH, PPC, EF1-α, and Tubulin-α. Specific primers for these genes were designed using Primer Premier 5.0 software (Table 1) and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China).

Real-time fluorescence quantification PCR

qRT-PCR reactions were performed using the TB Green™ Premix Ex Taq™ II (Tli RNaseH Plus) kit (TaKaRa, Japan). Each sample was run in triplicate technical replicates under the following cycling conditions: initial denaturation at 95 °C for 30 s; 45 cycles of denaturation at 95 °C for 5 s and annealing/extension at 60 °C for 30 s. The 20 μL reaction mixture comprised 10 μL of 2×SYBR Green Pro Taq HS Premix, 0.4 μL each of forward and reverse primers (10 μM), 1 μL of cDNA template, and nuclease-free water to volume.

Standard curve and amplification efficiency analysis

cDNA samples were serially diluted tenfold across a range of 10² to 10⁶ to construct standard curves. The amplification efficiency (E) for each reference gene was calculated using the formula E = 10^−1/slope. Each dilution point was assayed in triplicate technical replicates.

Processing and analysis of test data

The expression stability of candidate reference genes was evaluated using three well-known algorithms: geNorm, NormFinder, and BestKeeper. Prior to analysis, raw Ct values were directly imported into each algorithm without any additional data transformation or normalization, as all three algorithms are designed to accept raw Ct values for stability calculations.

geNorm^[19] calculates the stability value (M) of each candidate gene based on the relative expression levels derived from Ct values [converted to relative quantities using the formula 2^−ΔCt and determines the optimal number of reference genes required by computing the pairwise variation coefficient (V) between genes. NormFinder^[20] analyzes raw Ct values directly by considering both intra- and inter-group variances to derive a stability value for each gene, enabling a robust ranking of gene stability. No additional data transformation was applied for the values imported into NormFinder, as it operates on the linear scale of raw Ct values. BestKeeper^[21] evaluates each candidate gene individually by calculating the pairwise correlation coefficient (r), coefficient of variation (CV), and standard deviation (SD) directly from raw Ct data. It also performs comparative analyses among candidate genes to assess their expression consistency and suitability as reference controls. Similar to NormFinder, BestKeeper does not require or apply any data transformation prior to analysis. Evaluating the stability of candidate reference genes using multiple analytical algorithms is essential, as their apparent stability can vary substantially depending on the computational approach and underlying assumptions of each program. RefFinder^[22] is a comprehensive web-based platform specifically developed to assess and select suitable reference genes from large experimental datasets. It integrates the major currently available algorithms—geNorm, NormFinder, BestKeeper, and the comparative ΔCt method—to systematically compare and rank candidate reference genes. Based on the ranking outcomes generated by each algorithm, RefFinder assigns an appropriate weight to each gene and calculates the geometric mean of these weights to generate a final comprehensive stability ranking. The input data for RefFinder consist of stability rankings produced by the individual algorithms; no direct transformation of raw Ct values was performed by this tool.

Verification of internal reference genes

The key flowering regulatory gene PsFRL5 was selected as the target gene to validate the reliability of reference gene stability results. The selection of PsFRL5 was based on its putative involvement in flowering regulation. FRL genes are homologs of the FRI gene family, which are known to play important roles in the regulation of flowering time through modulation of FLC expression in model plants. Previous transcriptome studies on tree peony indicated that PsFRL5 exhibited differential expression patterns during floral bud differentiation and development stages, suggesting its potential role in flowering regulation. Therefore, PsFRL5 was considered an appropriate candidate gene for validating the stability and suitability of the selected reference genes in this study. Moreover, expression patterns of PsFRL5 were analyzed across early- and late-flowering cultivars, as well as during floral bud differentiation and flower development stages. For normalization, two of the most stable reference genes and the least stable reference gene were used to assess the impact of reference gene selection on expression profiling of target genes.

Total RNA was extracted from the petals of 21 early-flowering and 21 late-flowering tree peony cultivars at full bloom, as well as from the petals of the early-flowering cultivar 'FengDan' and the late-flowering cultivar 'LianHe' at various stages of floral bud differentiation and flower development. RNA integrity was assessed by performing 1% agarose gel electrophoresis, which showed distinct, intact 28S and 18S rRNA bands. The 28S band was approximately twice as intense as the 18S band (Fig. 3), confirming the high RNA quality. Quantification using Qubit 2.0 Fluorometer showed A260/A280 ratios between 2.0 and 2.1 for all samples, confirming a high RNA purity, which is suitable for downstream experiments.

Figure 3. 

Electrophoresis detection of total RNA. Marker: DNA marker DL2000; 1–21: flower development period of a tree peony.

Primer specificity test of internal reference genes

Primer specificity for the 23 candidate reference genes was first evaluated. Gradient PCR analyses were performed to determine the optimal annealing temperature for each primer pair, and a uniform annealing temperature of 60 °C was selected for all candidates. The resulting PCR products were subsequently examined by electrophoresis on 1% agarose gels. Clear and distinct amplification bands were observed for all 23 genes, with product sizes consistent with the expected lengths, ranging from 140 to 250 bp. No nonspecific amplification, primer–dimer formation, or additional bands were detected (Fig. 4). Furthermore, melting curve analyses revealed a single, sharp peak for each gene, confirming the high specificity of the primer pairs (Fig. 5).

Figure 4. 

Amplification results of each primer. Marker: DNA marker DL2000; 1–23: amplification results of the 23 reference genes to be tested.

Figure 5. 

The melting curve of each primer. a–w: Melting curve plots corresponding to each of the 23 tested candidate reference genes, each showing a single specific amplification peak without non-specific products or primer dimers.

Delta CT method analysis

Ct (quantification cycle): the number of amplification cycles at which the fluorescence signal exceeds the background threshold during qRT-PCR detection. A lower Ct value indicates a higher initial transcript abundance of the target gene. The expression levels of 23 candidate reference genes were detected by qRT-PCR, and their expression stability was analyzed using the ΔCt method. The results showed that the Ct values of all candidate reference genes were distributed within a certain range. A lower Ct value indicated a higher transcriptional expression level, and a smaller variation in Ct values among samples indicated more stable expression. In the full-bloom petals of early-and late-flowering tree peony cultivars, EF1-α-1 had the lowest average Ct value and the highest expression level, whereas RPL had the highest average Ct value and the lowest expression level. Tubulin-α exhibited the largest variation in Ct values among different cultivars and the poorest expression stability, while PPC showed the smallest variation in Ct values and the most stable expression (Fig. 6a). During the 16 stages of floral bud differentiation and flower development in Paeonia ostii 'FengDan' and Paeonia suffruticosa 'LianHe', EF1-α-1 had the lowest average Ct value and the highest expression level, while RPL had the highest average Ct value and the lowest expression level. UBi-lp displayed the largest fluctuation in Ct values across developmental stages and the worst stability, whereas GTP showed the smallest fluctuation in Ct values and the most stable expression (Fig. 6b).

Figure 6. 

Expression levels of the 23 candidate reference genes in a tree peony. (a) Expression levels in full-bloom petals of early- and late-flowering cultivars. (b) Expression levels during 16 stages of floral bud differentiation and flower development.

Standard curve and amplification efficiency analysis

Amplification efficiencies of all primers were evaluated using standard curves, yielding efficiencies between 92.8% and 110.1%, with correlation coefficients (R²) > 0.97 (Table 2). These results demonstrate that all primers exhibit high specificity and efficient amplification, meeting the essential criteria for qRT-PCR assays, and are thus suitable for subsequent experiments.

Screening of optimal internal reference genes for studying the flowering time of tree peony

The geNorm algorithm evaluates the stability of candidate reference genes using an internal calculation method, in which a lower M value indicates greater expression stability. We applied geNorm to assess the stability of 23 candidate reference genes across 21 early-flowering and 21 late-flowering cultivars. The stability ranking derived using geNorm had the following order: PPC/CSLD2 > KOR > AGD > eIF2 > GTP > UBi-lp > UPL > PPR > HDH > GAPDH > UBQ > UBC > eIF3 > EF1-α > HIS > IWS > ARP > PP2A > EF1-α-1 > CYP > Tubulin-α > RPL (Fig. 7a). Pairwise variation analysis showed a V_2/3 value of 0.265, exceeding the default threshold of 0.15 (Fig. 7b). Similarly, the values for V_3/4, V_4/5, and V_5/6 were 0.272, 0.191, and 0.163, respectively, all above the threshold. However, V_6/7 dropped below the threshold at 0.139, indicating that the top six genes—PPC, CSLD2, KOR, AGD, eIF2, and GTP—constitute a stable set of reference genes for this experimental condition.

Figure 7. 

Results of early- and late-flowering variety experiments using geNorm software. a: Average expression stability value (M) ranking of 23 candidate reference genes (genes arranged from least stable on the left to most stable on the right); b: Pairwise variation (V) analysis to determine the optimal number of reference genes for normalization.

The selection criteria of the NormFinder program align with those of geNorm, both identifying the reference gene with the lowest stability value as the optimal internal control. However, NormFinder uniquely selects a single best reference gene. According to NormFinder analysis, the stability ranking of the 23 candidate genes was AGD > PPC > eIF2 > GTP > KOR > CSLD2 > UBi-lp > UBQ > UPL > HDH > eIF3 > PPR > UBC > GAPDH > EF1-α > HIS > IWS > ARP > PP2A > EF1-α-1 > CYP > Tubulin-α > RPL, as per which AGD was identified as the optimal reference gene (Table 3). Notably, the top-ranked genes by NormFinder largely concurred with those identified by geNorm.

According to the BestKeeper program, the stability ranking of the 23 candidate reference genes was eIF3 > GTP > UBQ > eIF2 > KOR > UBC > UBi-lp > AGD > PPC > CSLD2 > HIS > ARP > PP2A > HDH > EF1-α-1 > EF1-α > GAPDH > UPL > PPR > Tubulin-α > CYP > IWS > RPL. The optimal reference gene identified by BestKeeper was eIF3 (Table 4), which differed from the results obtained by geNorm and NormFinder, likely due to differences in the underlying algorithms of the software. However, PPC, CSLD2, KOR, AGD, eIF2, and GTP remained highly ranked in the BestKeeper analysis, further confirming their stability across different algorithms.

The stability of 23 candidate reference genes was evaluated using the RefFinder program, with the following ranking: PPC > AGD > eIF2 > KOR > GTP > CSLD2 > eIF3 > UBQ > UBi-lp > HDH > UPL > UBC > PPR > GAPDH > HIS > EF1-α > ARP > PP2A > IWS > EF1-α-1 > CYP > Tubulin-α > RPL. Based on these calculations, PPC was identified as the most stable reference gene (Table 5), making it the most suitable reference gene for the study of early and late flowering in tree peony.

Screening of optimal reference genes in flower bud differentiation and flower development of tree peony

The stability ranking of 23 candidate reference genes during floral bud differentiation and flower development in Paeonia obtained using geNorm was as follows: EF1-α-1/UBi-lp > HDH > CYP > Tubulin-α > AGD > PPC > IWS > GTP > EF1-α > GAPDH > PPR > UBQ > eIF3 > UBC > KOR > CSLD2 > UPL > HIS > eIF2 > ARP > RPL > PP2A (Fig. 8a). Pairwise variation analysis yielded a V_2/3 value of 0.138, which is below the default threshold of 0.15 (Fig. 8b), indicating that the two most stable reference genes for this experimental condition are EF1-α-1 and UBi-lp.

Figure 8. 

Results of bud differentiation and flower development experiments using geNorm software. a: Average expression stability value (M) ranking of 23 candidate reference genes; genes arranged from left to right represent the least stable to the most stable. b: Pairwise variation (V) analysis to determine the optimal number of internal reference genes for normalization; the threshold value of V is set as 0.15.

The stability ranking of 23 candidate reference genes during floral bud differentiation and flower development in Paeonia according to the NormFinder analysis was GTP > EF1-α > GAPDH > EF1-α-1 > Tubulin-α > HDH > PPR > UBQ > IWS > PPC > UBi-lp > UBC > CYP > AGD > KOR > eIF3 > CSLD2 > ARP > HIS > UPL > eIF2 > RPL > PP2A, as per which GTP was identified as the optimal reference gene (Table 6).

According to the BestKeeper analysis, the stability ranking of 23 candidate reference genes during floral bud differentiation and flower development in Paeonia was RPL > GTP > EF1-α > UBC > PP2A > ARP > IWS > EF1-α-1 > CYP > AGD > HDH > UBQ > Tubulin-α > UBi-lp > KOR > GAPDH > eIF3 > PPC > HIS > PPR > UPL > CSLD2 > eIF2. Although RPL was identified as the optimal reference gene by BestKeeper (Table 7), it showed relatively poor stability in geNorm and NormFinder analyses.

The stability of 23 candidate reference genes was assessed using the RefFinder program, yielding the following ranking: GTP > EF1-α-1 > EF1-α > HDH > UBi-lp > GAPDH > Tubulin-α > CYP > IWS > UBC > RPL > AGD > PPC > UBQ > PPR > ARP > KOR > PP2A > eIF3 > CSLD2 > HIS > UPL > eIF2. Based on these results, GTP was identified as the most stable reference gene (Table 8), making it the most suitable reference gene for studies on flower bud differentiation and development in tree peony.

Stability verification of tree peony internal reference gene

To further validate the stability of the optimal reference genes, we selected the most important floral development-related gene PsFRL5 for verification during flower development in different tree peony cultivars. The results are shown in Fig. 9a and b. When using the most stable reference genes, PPC and GTP, as internal controls, the expression pattern of the target gene PsFRL5 was consistent across cultivars. However, when the least stable reference gene, PP2A (Fig. 9c), was used, the expression pattern of PsFRL5 differed.

Figure 9. 

Results of expression verification analysis of PsFRL5 in different varieties of peony flowers at 7 stages of development. a: PPC; b: GTP; c: PP2A. Notes: MU: mutant of Paeonia ostii 'Fengdan'; FD: Paeonia ostii 'Fengdan'; LH: Paeonia suffruticosa 'Lianhe'.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR) is a widely used method for analyzing gene expression patterns due to its sensitivity, specificity, and reproducibility^[23,24]. However, the accuracy of qRT-PCR results depends heavily on the selection of appropriate reference genes. No single reference gene is universally stable across all plant species, tissues, developmental stages, or experimental conditions^[25,26]. Therefore, it is essential to select the most appropriate reference genes tailored to the specific experimental conditions^[27,28]. The commonly used Actin gene has been employed as a reference gene in species such as Arabidopsis thaliana^[29], Populus^[30], and Oryza sativa^[31] to study various functional genes. In Oryza sativa, researchers analyzed public microarray databases and screened 10 candidate reference genes from four widely expressed gene families: GAPDH, ACT, UBI, and CYP (cytochrome P450). After conducting qRT-PCR experiments and rigorous statistical analysis, they found that a ubiquitin isoform gene performed best as a single reference gene, while a combination of another ubiquitin gene and two cyclophilin isoforms provided the most stable expression, offering a reliable basis for monitoring differential expression of drought-related candidate genes^[32]. Paeonia is a valuable ornamental and medicinal plant endemic to China, characterized by diverse flowering phenology and a complex floral development process, which poses significant challenges for selecting stable reference genes. In this study, we systematically evaluated the expression stability of 23 candidate reference genes across petals of different flowering cultivars, as well as during floral bud differentiation and development stages, using four widely accepted algorithms: geNorm, NormFinder, BestKeeper, and RefFinder. Our aim was to identify reliable internal controls to support future molecular investigations of flowering regulation in Paeonia.

The screening results indicated that PPC was the most stable reference gene for distinguishing between early- and late-flowering tree peony varieties. This conclusion was consistent across both geNorm and NormFinder analyses. Notably, AGD ranked first in NormFinder and PPC ranked second, both showing high stability. The agreement between the algorithms confirms the reliability of these genes as reference genes for comparative studies of early and late flowering varieties. PPC, involved in carbon metabolism and stress response, shows stable expression across different flowering varieties, potentially reflecting the conserved metabolic balance required for petal development^[33]. The PPC gene has been identified as the optimal internal reference gene for comparison between flowering periods, which is supported by solid biological evidence. PPC is also a key enzyme involved in carbon fixation and primary metabolism, playing a critical role in maintaining cellular energy homeostasis during flower development. Its stable expression across early-flowering and late-flowering varieties suggests that, regardless of variation in flowering time, the fundamental metabolic processes are tightly regulated, making it a reliable reference gene for comparative studies. CSLD2, a member of the cellulose synthase gene family, plays a crucial role in cell wall biosynthesis and cell expansion, both of which are essential processes during petal formation^[34]. The stable expression of these functional genes across varieties suggests their potential as species-specific reference genes for tree peony.

For the dynamic process of flower bud differentiation and development, geNorm identified EF1-α-1 and UBi-lp as the most stable reference genes, while both NormFinder and BestKeeper confirmed GTP as one of the optimal candidate genes. EF1-α-1 is a widely recognized reference gene across various plant species due to its conserved role in protein synthesis and its stable expression throughout different developmental stages^[35]. UBi-lp is involved in the ubiquitination pathway, which regulates protein turnover and is crucial for developmental transitions, explaining its stable expression during the development of tree peony flowers^[36]. GTP plays a role in signal transduction and energy metabolism, making it a reliable reference gene for tracking gene expression dynamics throughout the flowering process^[37]. Furthermore, GTP is involved in signal transduction and energy metabolism. Its stability during the dynamic process of flower bud differentiation reflects the conserved demand for cellular signal transduction and energy supply throughout the developmental transition. The stable expression of these functional genes across different varieties suggests their potential as species-specific reference genes for peony.

Our findings align with recent studies on reference gene selection in other ornamental and perennial species, highlighting both commonalities and species-specific patterns that underscore the novelty of this work. In tree peony, Zhou et al.^[38] evaluated nine reference genes in Paeonia suffruticosa 'Doulv' and identified GAPDH and ACT as the most stable in petaloid stamens, while UBC and MBF1A performed best across tissues. Notably, these genes were not top-ranked in our study, reflecting the condition-dependent nature of reference gene stability. Similarly, in Clematis lanuginosa, Li et al.^[39] found PP2A-2 and UBC34 stable in vegetative tissues, whereas PP2A ranked 19th in our study on early- and late-flowering cultivars. Collectively, these comparisons demonstrate the following: (i) no universal reference gene exists across perennial species or experimental conditions; (ii) PPC and GTP are novel reference genes for tree peony; (iii) the condition-specific reference genes validated here provide essential tools for future functional genomics research on flowering regulation in Paeonia species.

Although this study is relatively comprehensive, several limitations should be addressed. First, our evaluation was focused solely on petal tissue, as petals are the primary target organs for studies on flowering time and flower development. However, the stability of reference genes may vary across different tissues (such as leaves, roots, stems, and reproductive organs), and future studies should expand validation to include a broader range of tissue types to support more extensive applications. Second, while this study covered two key experimental contexts—varietal comparison of flowering time and flower bud differentiation/flower development—we did not assess the stability of reference genes under abiotic stress conditions (such as drought, low temperature, and salt stress) or hormone treatments. Therefore, future research should validate the identified reference genes (PPC and GTP) across a wider array of tissues and experimental conditions, including abiotic stress treatments and plant hormone treatments.

In conclusion, this study, through transcriptomic data mining of tree peony, identified suitable reference genes for the regulation of flowering time and flower development. The most appropriate reference gene for early- and late-flowering varieties was found to be PPC, while GTP was deemed the most suitable reference gene for flower bud differentiation and development. These findings provide a solid foundation for accurate gene expression analysis in tree peony and will contribute to advancing future molecular breeding and functional genomics research.

This study systematically screened and validated 23 candidate reference genes for qRT-PCR normalization in tree peony during flowering regulation and floral development. Comprehensive analysis using geNorm, NormFinder, BestKeeper, and RefFinder revealed that PPC was the most stable reference gene for comparing early-and late-flowering cultivars, while GTP performed best across floral bud differentiation and flower development stages. Validation using the flowering-related gene PsFRL5 confirmed the reliability of these reference genes. These results provide robust, condition-specific internal controls for accurate gene expression analysis in tree peony. The identified reference genes will facilitate future molecular studies on flowering time regulation, floral development mechanisms, and functional genomics in Paeonia species. This work also establishes a reproducible pipeline for reference gene selection in other woody perennial ornamentals with similar flowering characteristics.