A disease-free and vigorous “Chang’e” O. fragrans tree from Xianning (Hubei, China) was selected as the mother tree. In May 2023, semi-lignified branches measuring 12–16 cm, each with 2–3 lateral buds, were harvested as cuttings. The cuttings were first soaked in distilled water for 12 hours, then disinfected with 1% calcium hypochlorite for 10 minutes, and subsequently rinsed three times with distilled water. Following disinfection, the cuttings were immersed in a 0.1 g L^–1 GGR rooting powder solution for 4 hours. They were then transplanted into a mixed soil matrix composed of peat, perlite, vermiculite, and yellow sand in a 1:1:1:1 ratio. The cuttings were cultivated at the O. fragrans base in Xianning.

Samples of O. fragrans were collected from various tissues, including roots, leaves, seeds, and flowers. A healthy and pest-free O. fragrans plant from Huazhong Agricultural University (114°21’W, 30°29’N) was selected for the study. Flower tissues were harvested at six developmental stages: S1 (stalk stage), S2 (early flowering stage), S3 (mid-flowering stage), S4 (full flowering stage), S5 (late flowering stage), and S6 (petal shedding stage) (Chen et al., 2021). Samples were collected at 10:00 am, with each stage represented by three biological replicates. The samples were rapidly frozen in liquid nitrogen and stored at –80°C for subsequent analysis.

In March 2024, healthy and uniformly growing O. fragrans plants from the campus of Hubei University of Science and Technology (114°19’52’’E, 29°51’19’’N) were subjected to acclimation culture. In April 2024, 81 plants demonstrating uniform growth were selected for nine different stress treatments. For abiotic stress treatments, the plants were exposed to 4°C to simulate cold stress, and treated with 300 mM NaCl and 20% PEG-6000 to induce salt and drought stress, respectively (Zhang et al., 2021c, d). For hormone treatments, the plants were sprayed with 300 μM ABA, 300 μM MeJA, and 5 mM ethephon, respectively (Zhang et al., 2021c, d). Metal ion treatments were induced by spraying the cuttings with 3 mM CuSO₄·5H₂O, 3 mM AlCl₃·6H₂O, and 3 mM FeSO₄, respectively (Khalid et al., 2020; Huang et al., 2021; Soares et al., 2022). Except for the cold stress treatment, each plant received 200 mL of the respective treatment solution, ensuring that all leaves were thoroughly wetted. The plants were then cultured in a light-controlled growth chamber at 25°C with a 12-hour light/12-hour dark cycle and 60% humidity. Three plants were used for each treatment, with three independent biological replicates (3 × 3 plants). Samples were collected at 0, 3, 6, 12, 24, and 72 hours post-treatment (Zhang et al., 2021c, d) and stored at –80°C for further analysis.

Common RGs such as 18S, ACT11, α-tubulin 5 (TUA5), U6, and UBQ4 previously identified in other species were locally blasted against the whole genome data of O. fragrans to identify candidate RGs with high homology. qRT-PCR primers were designed using Primer Premier 5.0 (Premier Biosoft International, Palo Alto, CA, USA), adhering to the following parameters: a PCR product length of 80–250 bp, a melting temperature (T_m) of 58–62°C, and a GC content of 40–60%. Primer specificity was subsequently evaluated using e-PCR (https://yanglab.hzau.edu.cn/OfIR/tools/epcr/).

Based on small RNA-seq data, nine miRNAs with relatively high expression levels and fold changes < 1.4 were selected as candidates for RGs. These include five miRNAs (ofr-miR159b-3p, ofr-miR168b-5p, ofr-miR171a-3p, ofr-miR395e, and ofr-miR403-3p), as well as four novel miRNAs (novel2, novel3, novel8, and novel33) ( Supplementary Table S1 ). In the miRNA sequences, uracil (U) was replaced with thymine (T), and modifications were made to achieve a T_m of 63–65°C, either by adding “G/C” bases at the 5’ end or deleting bases at the 5’/3’ ends. Reverse primers were provided by the miRcute miRNA qPCR Detection Kit (SYBR Green) (TianGen Biotech, Beijing, China), while the remaining primers ( Table 1 ) were synthesized by Tsingke Biotech Co., Ltd. (Nanjing, Jiangsu Province, China).

Total RNA was extracted from each sample using the HiPure Plant RNA Mini Kit (Magen Biotechnology, Guangzhou, Guangdong, China), following the manufacturer’s instructions. RNA integrity was evaluated using 1% agarose gel electrophoresis, and concentration was determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). RNA samples that met the quality criteria were then used for cDNA synthesis. Specifically, 0.8 µg of RNA was reverse transcribed into cDNA using the miRcute miRNA First-Strand cDNA Synthesis Kit (Tiangen Biotech). The resulting cDNA was stored at –20°C.

Accurately measured 2 μL of cDNA from each sample were combined to assess primer specificity through RT-PCR amplification using the Fast PCR Kit (Vazyme). The reaction mixture comprised 1.2 μL each of forward and reverse primers (10 μM), 1.2 μL of cDNA, 10 μL of 2× Rapid Taq Master Mix, and 6.4 μL of ddH₂O. The amplification protocol involved an initial denaturation at 95°C for 3 minutes, followed by 35 cycles of denaturation at 95°C for 15 seconds, annealing at 59°C for 15 seconds, and extension at 72°C for 15 seconds, with a final extension at 72°C for 5 minutes. Primer specificity was subsequently verified using 2.0% (w/v) agarose gel electrophoresis.

The cDNA template was serially diluted in 5-fold increments (1:4, 1:24, 1:124, 1:624, 1:3124; cDNA:water, v:v). qRT-PCR was performed using a 20 μL reaction mixture, prepared according to the instructions of the miRcute Plus miRNA qPCR Detection Kit (SYBR Green). The reaction mixture comprised 10 μL of 2× miRcute Plus miRNA Premix (SYBR&ROX), 0.4 μL each of forward and reverse primers (10 μM), 2 μL of 10-fold diluted miRNA first-strand cDNA, and 7.2 μL of RNase-free ddH₂O. The qRT-PCR was performed on the Tianlong Gentier 96E system (Tianlong Technology Co., Ltd., Xi’an, China) using the following program: an initial denaturation at 95°C for 15 minutes, followed by 40 cycles of denaturation at 94°C for 20 seconds, and annealing and extension at 60°C for 30 seconds. A melting curve analysis was performed from 60°C to 95°C to detect primer dimers and other amplification artifacts. Each gene assay included a non-template control, with reactions performed in triplicate biological replicates and triplicate technical replicates. Correlation coefficients (R²) and amplification efficiencies (E values) were calculated from the qRT-PCR data (Whistler et al., 2010).

cDNA samples were diluted 10-fold (1:9, cDNA:water, v:v) and analyzed by qRT-PCR to obtain the raw Ct values. The expression stability of the RGs was assessed using the following algorithms: delta-Ct (Silver et al., 2006), geNorm (v3.5) (Vandesompele et al., 2002), NormFinder (v0.953) (Andersen et al., 2004), and BestKeeper (v1.0) (Pfaffl et al., 2004). Specifically, geNorm and NormFinder converted the Ct values to 2^–ΔCt (where delta-Ct = Ct value – minimum Ct value for each group) for stability assessment. BestKeeper evaluated stability based on the E values derived from the Ct values using the LinRegPCR program, considering the coefficient of variation (CV) and standard deviation (SD). Additionally, geNorm determined the optimal number of RGs by calculating the pairwise variation (V_n/V_n+1) between consecutive normalization factors.

The geometric mean ranking was calculated based on the average rankings of genes across various treatments, tissues, or all samples, as determined by the four algorithms. Additionally, the stability analysis results of the RGs were comprehensively validated using RefFinder (https://blooge.cn/RefFinder/) (Xie et al., 2023).

To evaluate the reliability of the selected RGs, we used the most stable RG combination and the least stable RG to normalize the expression levels of two miRNAs, ofr-miR166e-5p and ofr-miR396b-3p, under each experimental condition using qRT-PCR. Primers for these miRNAs were designed using the same method as for the candidate reference miRNAs, with the specific primers listed in Table 1 . The relative expression levels of the miRNAs were calculated using the 2^–ΔΔCt method (Livak and Schmittgen, 2000).

In this study, a total of 14 candidate RGs were selected for gene normalization analysis, comprising 5 genes and 9 miRNAs ( Table 1 ). The E values of the candidate RG primers ranged from 91.356% to 111.362%, with R2 values ≥ 0.976 ( Table 1 ). Primer specificity was evaluated using agarose gel electrophoresis and melting curve analyses. The agarose gel electrophoresis results demonstrated that all primers successfully amplified PCR products of the expected size, confirming the correct design and specificity of the primers ( Figure 1A ). Melting curve analysis further verified that each primer produced a single melting peak ( Figure 1B ), suggesting that the amplified products were free from non-specific amplification or primer dimers.

To evaluate the suitability of the candidate RGs, their expression levels were analyzed across all samples using qRT-PCR. These samples included those exposed to various conditions: abiotic stress (low temperature, drought, and salt stress), hormone treatments (ABA, MeJA, and ethephon), metal ion treatments (Fe²⁺, Al³⁺, and Cu²⁺), different tissues (roots, leaves, seeds, and flowers), and during the flower opening and senescence ( Figure 2 ).

Under abiotic stress, the average Ct values of the samples ranged from 15.472 (ofr-miR159b-3p) to 29.802 (TUA5) ( Figure 2A ). Notably, novel8, novel2, and ofr-miR159b-3p exhibited relatively low Ct values, indicating higher initial copy numbers and thus higher expression levels in the samples (Bustin et al., 2009). Among these, novel8 demonstrated the smallest variation in Ct values, with a range of only 2.772 ( Figure 2A ). Under hormone treatments, the average Ct values varied from 15.625 (ofr-miR159b-3p) to 31.180 (TUA5) ( Figure 2B ). Among these, ofr-miR159b-3p displayed the smallest variation in Ct values, with a range of just 2.304 ( Figure 2B ). For metal ion treatments, ofr-miR171a-3p exhibited the least variation in Ct values, with a range of 3.164 ( Figure 2C ). During flower opening and senescence, novel33 had the smallest Ct value variation (1.000) ( Figure 2D ). In different tissues, ofr-miR168b-5p showed the least variation in Ct values (1.718) ( Figure 2E ). Across all samples, U6 had the smallest variation in Ct values, with a range of 6.758 ( Figure 2F ). These results highlight the relative stability of certain genes in their expression levels, suggesting their potential as stable RGs. However, further in-depth analysis is required to confirm their suitability.

To evaluate the stability of the candidate RGs, we employed the delta-Ct method, which compares the relative expression levels of gene across various sample groups. The RGs were ranked according to the reproducibility of the average standard deviation (STDEV) of gene expression differences between samples (Silver et al., 2006). Generally, a lower STDEV value indicates greater stability in gene expression.

Under low temperature stress, abiotic stress, MeJA treatment, hormone treatments, and across all samples, ofr-miR159b-3p exhibited the most stable expression ( Figures 3A, D, G, H, O ). Novel2 was more stable than other RGs under salt stress, Cu²⁺ treatment, and metal ion stress ( Figures 3C, J, L ). Novel3 was the most stable gene under ABA and Al³⁺ treatments ( Figures 3E, I ). Novel33 showed the highest stability during ethephon treatment, Fe²⁺ treatment, and throughout the flower opening and senescence stages ( Figures 3F, K, M ). Under drought stress and across different tissues, ofr-miR168b-5p and ofr-miR395e were identified as the most stable RGs, respectively ( Figures 3B, N ). Notably, TUA5 (10, 66.67%) and ACT11 (10, 66.67%) exhibited the least stability across most conditions ( Figure 3 ; Supplementary Table S2 ).

The geNorm algorithm evaluates the stability of RGs by calculating the mean variation value, commonly known as the M value, where a lower M value indicates greater gene stability (Vandesompele et al., 2002). The results of the geNorm analysis are illustrated in Figure 4 . ofr-miR159b-3p demonstrated the highest stability under cold stress, Fe²⁺ treatment, and metal ion stress ( Figures 4A, K, L ). Novel8 was identified as the most stable RG under drought stress, abiotic stress, MeJA treatment, hormone treatments, Cu²⁺ treatment, and during flower opening and senescence ( Figures 4B, D, G, H, J, M ). Novel3 exhibited the greatest stability under ABA and Al³⁺ treatments ( Figures 4E, I ). Across various tissues and all samples, novel2 displayed superior stability compared to other RGs ( Figures 4N, O ). Under salt stress and ethephon treatment, ofr-miR403-3p and novel33 were the most stable, respectively ( Figures 4C, F ). Notably, consistent with the delta-Ct analysis, TUA5 (10, 66.67%) and ACT11 (10, 66.67%) were the least stable across most conditions ( Figures 4A–O ; Supplementary Table S2 ).

A single RG is generally insufficient for achieving the stability required for accurate normalization. Therefore, it is essential to use two or more RGs to minimize error and ensure more reliable results. The geNorm algorithm determines the optimal number of candidate RGs by calculating the pairwise variation (V_n/V_n+1), with a threshold of 0.15 to identify the ideal number of RGs for normalization. For PEG treatment, hormone treatments, Fe²⁺ treatment, Al³⁺ treatment, and Cu²⁺ treatment, the V₃/V₄ value was below 0.15, indicating that a combination of three RGs is most appropriate ( Figure 4P ). For MeJA treatment and abiotic stress, the V₄/V₅ value was below 0.15, suggesting that the optimal RG combination includes four genes ( Figure 4P ). In metal ion treatments, the use of five RGs is necessary ( Figure 4P ). Across all samples, nine RGs are required for accurate normalization ( Figure 4P ).

NormFinder was employed to evaluate the expression stability of candidate RGs by ranking them according to variance analysis results (Andersen et al., 2004). A low stability value indicates that the gene exhibits significant variability under different experimental conditions, making it unsuitable as a RG. In contrast, a high stability value indicates that the gene shows minimal variation, thus demonstrating higher stability and being more reliable as a RG. The NormFinder analysis identified the most stable RGs across various experimental conditions, as depicted in Figure 5 . Under cold, drought, and salt stress conditions, the most optimal RGs were U6, ofr-miR168b-5p, and novel8, respectively ( Figures 5A–C ). For Al³⁺ and Cu²⁺ treatments, as well as across all samples, 18S, novel2, and novel3 exhibited the highest stability compared to other genes ( Figures 5I, J, O ). Similarly, under abiotic stress, MeJA treatment, and hormone treatments, ofr-miR159b-3p demonstrated the greatest stability ( Figures 5D, G, H ). Under ABA treatment, ethephon treatment, and the flower opening and senescence process, novel33 was identified as the most stable RG ( Figures 5E, F, M ). For Fe²⁺, metal ion treatments, and different tissues, ofr-miR395e emerged as the optimal RG ( Figures 5K, L, N ). It is noteworthy that, consistent with the delta-Ct and geNorm analyses, TUA5 (10, 66.67%) and ACT11 (9, 60.00%) exhibited the lowest stability across most conditions ( Figure 5 ; Supplementary Table S2 ).

BestKeeper evaluates the stability of RGs by calculating the SD and CV of the Ct values (Pfaffl et al., 2004). RGs with an SD of less than 1.0 are considered to exhibit stable expression, as a lower CV generally indicates greater stability. Under cold stress, ABA treatment, and Cu²⁺ treatment, ofr-miR159b-3p was identified as the most stable RG ( Figures 6A, E, J ). Novel3 proved to be the optimal RG under drought stress, salt stress, abiotic stresses, ethephon treatment, MeJA treatment, hormone treatments, Fe²⁺ treatment, and across all samples ( Figures 6B–D, F–H, K, O ). Novel33 exhibited the highest stability during Al³⁺ treatment, metal ion treatments, and throughout the processes of flower opening and senescence ( Figures 6I, L, M ). Additionally, ofr-miR168b-5p displayed superior stability compared to other RGs across various tissues ( Figure 6N ). Consistent with the delta-Ct, geNorm, and NormFinder analyses, TUA5 (12 instances, 75.00%) and ACT11 (10 instances, 66.67%) showed the lowest stability under most conditions ( Figure 6 ).

Different algorithms utilize distinct principles to evaluate gene stability, often resulting in variations in their rankings. To achieve a more comprehensive assessment, we calculated the geometric mean of stability scores across delta-Ct, geNorm, NormFinder, and BestKeeper for each experimental condition. The optimal RG combinations were subsequently selected based on geNorm analysis. The results were as follows: for ethephon treatment: novel33 and ofr-miR159b-3p; for ABA treatment: novel3 and ofr-miR159b-3p; for MeJA treatment: ofr-miR159b-3p, novel8, novel2, and novel3; for hormone treatments: ofr-miR159b-3p, novel8, and novel3; for drought stress: ofr-miR168b-5p, ofr-miR159b-3p, and novel2; for salt stress: novel8 and novel3; for cold stress: ofr-miR159b-3p and ofr-miR403-3p; under abiotic stress: ofr-miR159b-3p, novel8, ofr-miR403-3p, and novel2; for Fe²⁺ treatment: novel3, novel33, and ofr-miR159b-3p; for Al³⁺ treatment: novel3, ofr-miR395e, and novel33; for Cu²⁺ treatment: novel2, novel8, and ofr-miR159b-3p; for metal ion stress: novel3, ofr-miR159b-3p, novel2, novel33, and ofr-miR395e; across various tissues: novel2 and ofr-miR395e; during flowering stages: novel33 and ofr-miR395e; and across all samples: ofr-miR159b-3p, novel3, novel2, novel33, novel8, U6, ofr-miR403-3p, ofr-miR395e, and ofr-miR168b-5p ( Table 2 ). Additionally, the least stable genes identified were as follows: TUA5 for cold stress, salt stress, abiotic stress, ABA treatment, Al³⁺ treatment, Fe²⁺ treatment, metal ion stress, and across all samples; UBQ4 for drought stress, and hormone treatments; ACT11 for MeJA, and ethephon treatments; ofr-miR171a-3p for different tissues, and during flowering stages; and ofr-miR168b-5p for Cu²⁺ treatment ( Table 2 ).

The optimal RG combinations and the least stable genes identified by RefFinder are largely consistent with those obtained from the geometric mean analysis. The only exception is under drought stress, where the most stable RG combination was identified as ofr-miR168b-5p, ofr-miR159b-3p, and novel2 (rather than novel8) ( Table 2 ; Supplementary Table S3 ).

To ensure the accuracy of stable RG expression, we evaluated the relative expression levels of two target genes, ofr-miR166e-5p ( Figure 7 ) and ofr-miR396b-3p ( Figure 8 ), using both stable and unstable RG combinations under various experimental conditions.

When the most stable RGs were employed, the expression patterns and levels of the target genes remained consistent across various experimental conditions ( Figures 7 , 8 ). In contrast, the use of unstable RGs resulted in significant discrepancies in the expression patterns of these target genes ( Figures 7 , 8 ). For example, under cold stress, the use of stable RGs indicated that ofr-miR166e-5p and ofr-miR396b-3p exhibited higher expression levels at 0 and 72 hours ( Figures 7A , 8A ), whereas the use of the unstable RG TUA5 led to peak expression levels of these target genes at 6 hours ( Figures 7A , 8A ). Similarly, under drought treatment, ofr-miR166e-5p and ofr-miR396b-3p reached their highest expression levels at 3 hours with stable RGs, while peak expression occurred at 12 hours when unstable RGs were used ( Figures 7C , 8C ). Across various treatments, although the expression patterns of the target genes were generally similar when using different stable RGs, their expression levels varied between treatments ( Figures 7 , 8 ). For instance, during flowering and senescence, stable RGs such as novel2 and ofr-miR395e, either individually or in combination, showed that the target genes had the highest expression in seeds, with lower levels in flowers and roots. In contrast, using the unstable RG ofr-miR171a-3p resulted in similar expression trends, but expression levels in flowers were only 14.43% and 3.32% of those observed with stable RGs ( Figures 7K , 8K ). Overall, the selection of appropriate RGs is crucial for accurate normalization of target gene expression, as incorrect RG selection can lead to inaccurate estimates of relative expression levels.