Selection of endogenous reference miRNAs for quantitative Real-Time PCR analysis under drought and salt stress in Paeonia suffruticosa 'Luoyang Hong'

Tree peony (Paeonia suffruticosa Andr.), a species in section Moutan DC. of the family Paeoniaceae, is a highly valuable multipurpose plant. It is widely appreciated for its ornamental flowers, oil-rich seeds, and medicinal roots, demonstrating significant added value^[1].

Salt and drought stress critically constrain plant growth by disrupting key physiological processes such as photosynthesis, respiration, and stomatal regulation. These disruptions lead to reduced germination, stunted growth, leaf wilting, and metabolic disorders^[2]. In tree peony, severe stress induces leaf wilting, reactive oxygen species (ROS) accumulation, and damage to cell membranes and chloroplasts, thereby impairing photosynthesis and plant development^[3]. Furthermore, such stress conditions inhibit leaf and flower growth, suppress the synthesis of organic compounds and anthocyanins, and ultimately reduce seed yield.

However, the molecular mechanisms underlying stress responses in tree peony remain poorly understood, with many stress-related genes and miRNAs yet to be identified. Accurately quantifying the up- or down-regulation of gene expression under stress conditions is essential, which necessitates the use of reliable reference genes for normalization. MiRNAs play crucial roles in plant responses to drought and salt stress. MiRNAs contribute to plant tolerance against salt and drought stress by regulating key processes such as hormone signalling (e.g., ABA and auxin), antioxidant defense, osmoregulation, and photosynthesis^[4,5].

Quantitative real-time PCR (qRT-PCR) is a highly sensitive and cost-effective method widely used to quantify miRNA expression^[6−8]. Accurate normalization with stably expressed ER genes is a crucial prerequisite for reliable results, as it minimizes technical errors in the quantification process^[9,10]. Since the suitability of ER genes varies depending on plant species and specific stress conditions^[11], it is necessary to systematically identify appropriate reference genes for tree peony under drought and salt stress.

Currently, the identification of ER genes in tree peony has been reported in contexts such as flower development, seed development, hormone treatments, and cultivar differences^[12−15]. In contrast, studies on ER miRNAs remain limited, with existing work focusing on early- and late-flowering varieties, and various tissues during bud and flower development^[16,17]. Evidence suggests that under specific conditions, ER miRNAs can demonstrate greater stability than conventional ER genes^[18]. Compared to traditional housekeeping genes such as U6, 5S, 5.8S^[11,19], the use of miRNA as an ER gene offers several distinct advantages in plant gene expression analysis. The primary advantage lies in its exceptional expression stability: many miRNAs remain consistently expressed across different tissues, developmental stages, and abiotic stress conditions such as drought and salinity, with minimal influence from external treatments. This provides a more reliable baseline for quantitative analysis. Additionally, miRNAs exhibit high technical consistency: when studying other miRNA targets, using a miRNA as an ER ensures uniformity across the entire experimental process, including RNA extraction, cDNA synthesis via the stem-loop primer method, and qPCR amplification. This effectively avoids the systematic errors that can arise when normalizing miRNA data using mRNA reference genes^[20].

In this study, leaves of P. suffruticosa 'Luoyang Hong' subjected to varying degrees of drought and salt stress were used as materials. Based on miRNA sequencing data^[21], 17 miRNAs with relatively stable expression were selected as candidates. Using the commonly applied snRNA U6 as a reference, their expression stability was assessed by qRT-PCR and systematically evaluated with four algorithms: geNorm, NormFinder, BestKeeper, and RefFfinder.

The five-year-old potted P. suffruticosa 'Luoyang Hong' was used as the experimental material (Fig. 1). Before the experiment began, all plants were fully watered to ensure normal growth. For the drought treatment, watering was withheld from 0 to 12 d to allow the plants to dry naturally, and leaf samples were collected at 9:00 a.m. on 0, 6, and 12 d. After the sampling on 12 d, rewatering was carried out by fully watering all plants. The rewatering treatment lasted for 6 d, after which the final sampling was conducted. For the salt stress treatment, on 0 d, each potted plant was irrigated with 1 L of 2.5 mol·L⁻¹ NaCl solution. Leaf samples were collected at 9:00 a.m. on 0, 1, 2, and 3 d. Each experimental treatment was conducted with three biological replicates, wherein three leaves from the middle section of each plant were pooled to form one analytical sample. Immediately after sampling, the leaves were frozen in liquid nitrogen and stored at −80 °C for subsequent total RNA extraction.

Figure 1. 

Morphological differences of P. suffruticosa 'Luoyang Hong' plants under drought and salt stress. The numbers indicate days after treatment. The horizontal bar represents 30 cm.

miRNA isolation and cDNA synthesis

Total miRNA was extracted using the miRcute Plant miRNA Isolation Kit (Omega, Guangzhou, China). The quality of miRNA extraction was detected by 2% agarose gel electrophoresis using DL2000 DNA Marker (AG, Hunan, China) as an indicator. The concentration of miRNA was measured using the NanoDrop 1000 spectrophotometer (Implen, Germany). The first-strand cDNA of miRNA was synthesized using the miRNA cDNA First-Strand Synthesis Kit (AG, Hunan, China), which reverse-transcribed total miRNA into cDNA.

Potential reference genes selection

The qRT-PCR primers were designed using Primer 5.0 software, and then synthesized by Sangon Biotech Co., Ltd. (Shanghai, China) (Table 1).

RT-PCR and qRT-PCR conditions

Using the cDNA from untreated P. suffruticosa 'Luoyang Hong' leaves as the template, PCR amplification was performed with qRT-PCR primers using Prime STAR LongSeq DNA Polymerase (Takara, Japan). The reaction system (25 μL) consisted of: Prime STAR LongSeq DNA Polymerase 12.5 μL, Primer F (10 μM) 1 μL, Universal Primer R (10 μM) 1 μL, cDNA 2 μL, and ddH₂O 8.5 μL. The reaction program was as follows: 95 °C for 2 min; 95 °C for 30 s, 55 °C for 30 s, 72 °C for 10 s, for 35 cycles; 72 °C for 2 min. The PCR products were detected on a 2% agarose gel using DL2000 DNA Marker (AG, Hunan, China) as an indicator.

Using the cDNA from P. suffruticosa 'Luoyang Hong' leaves under different days and treatments as the template, qRT-PCR was performed with qRT-PCR primers using the Biorad CFX96 instrument and SYBR Green Premix Pro Taq HS qPCR Kit (AG, Hunan, China). The reaction system was: 2 × PCR SYBR Green Premix Pro Taq HS Mix 5 μL, Primer F (10 μM) 0.2 μL, Universal Primer R (10 μM) 0.2 μL, cDNA 1 μL, and ddH₂O 3.6 μL. The reaction program was: 95 °C for 30 s; 95 °C for 5 s, 60 °C for 30 s, for 40 cycles. After the reaction, the Ct values of each candidate ER miRNAs were analyzed for all cycles (65–95 °C). To ensure the repeatability of the results, reverse transcription-quantitative PCR was carried out in three biological replicates × three technical replicates.

Expression stability analysis of reference miRNAs

The stability of the 17 candidate ER miRNAs and U6 was evaluated and ranked using GeNorm, NormFinder, and Bestkeeper software, respectively. Subsequently, the online tool RefFinder (http://blooge.cn/RefFinder/) was utilized to comprehensively evaluate and analyze the results obtained from the aforementioned three software programs, ultimately screening out the most suitable ER miRNA under drought and salt stress in tree peony.

The GeNorm tool identifies the most stable ER genes and their optimal number by analyzing the average expression stability value (M value) and the pairwise variation (V_n/V_n+1 value). Its fundamental principle is that the expression ratio of two ideal ER genes should remain constant across different template groups. Specifically, when the M value is less than 1.5, it indicates relatively stable gene expression. Moreover, there is an inverse relationship between the M value and the expression stability of the gene; that is, the lower the M value, the more stable the expression^[22,23]. The pairwise variation (V_n/V_n+1 value) is used to determine the optimal number of ER genes. Typically, 0.15 is used as the critical value for judgment; when the V_n/V_n+1 value is less than 0.15, it indicates that the optimal number of ER genes is 'n'^[20,24].

NormFinder assesses the stability of ER genes by analyzing both intergroup and intragroup differences in expression levels. The stability value (SV) it calculates is inversely related to the stability of the ER gene^[25]. Bestkeeper evaluates the stability of ER genes by analyzing the coefficient of variance (CV) and standard deviation (SD). Lower CV and SD values indicate stronger stability of the ER gene. RefFinder integrates the results from the above three software programs, calculates the weights using the geometric mean method, and ultimately screens out the most stable ER miRNA under osmotic stress conditions in tree peony.

Using cDNA as the template and employing primers for PCR amplification, the electrophoresis detection results (Fig. 2) show that the amplification results of the 17 candidate ER miRNAs and U6 primers all appear as single bands, with lengths consistent with expectations. This indicates that the primers are specific and can be used for subsequent experiments.

Figure 2. 

Amplified gel electrophoresis of 18 candidate RGs in Table1. M: 2,000 bp Marker; 1−18: miRNAs in Table 1 and the order are kept consistent.

Delta CT method analysis

qRT-PCR was used to analyze the expression levels of the 17 candidate ER miRNAs and U6 under drought stress. The results showed that the Ct values of the 17 candidate ER miRNAs and U6 ranged from 15.45 to 26.25. Among them, the one with the lowest average Ct value was peu-miR2916_R-1_2ss6AG21TG, and the one with the highest average Ct value was csi-miR1515a. The gene with the largest difference in Ct values among different samples under drought stress was mdm-MIR482a-p3_2ss11CG20TA, while the one with the smallest difference in Ct values was ptc-MIR156b-p3 (Fig. 3a). The analysis of the expression levels of the 17 candidate ER miRNAs and U6 under salt stress revealed that their Ct values ranged from 13.94 to 25.11. Among them, the one with the lowest average Ct value was peu-miR2916_R-1_2ss6AG21TG, and the one with the highest average Ct value was csi-miR1515a. The gene with the largest difference in Ct values among different samples under salt stress was ptc-MIR156b-p3, while the one with the smallest difference in Ct values was PC-3p-160_22001 (Fig. 3b).

Figure 3. 

Expression levels of miRNAs in 18 candidate parameters. (a) Expression levels under drought stress. (b) Expression levels under salt stress. Box-whisker plot showing the CT variation among 18 test samples. A line across the box depicts the median. In each box, the upper and lower edges indicate the 25^th and 75^th percentiles. Whiskers go from the minimal to maximal value or, if the distance from the first quartile to the minimum CT value is more than 1.5 times the interquartile range (IQR), from the smallest value included within the IQR to the first quartile.

GeNorm analysis

The analysis results from the GeNorm software indicate that under drought stress, the M values of the 17 candidate miRNAs and U6 are all less than 1.5, demonstrating good stability. Among them, ath-miR858a_L-1R+1 and mtr-miR166a_L+1R-1 exhibit the best stability (Fig. 4a). The pairwise variations of the 18 candidate ER miRNAs under drought stress show that V2/3 = 0.128 < 0.15, indicating that the optimal number of ER miRNAs for the quantitative expression of miRNAs is two, namely ath-miR858a_L-1R+1 and mtr-miR166a_L+1R-1 (Fig. 4b).

Figure 4. 

(a) Expression stability of 18 candidate reference miRNAs analyzed in drought stress. The most stable miRNA is shown on the right and the least stable miRNA is shown on the left. (b) Pairwise variations of 18 candidate reference miRNAs in drought stress. (c) Expression stability of 18 candidate reference miRNAs analyzed in salt stress. The most stable miRNA is shown on the right and the least stable miRNA is shown on the left. (d) Pairwise variations of 18 candidate reference miRNAs in salt stress.

Under salt stress, the M values of the 18 candidate miRNAs are also all less than 1.5, indicating good stability. Among them, mtr-miR166e-5p_1ss10GT and PC-3p-160_22001 exhibit the best stability, while ptc-MIR156b-p3 is the most unstable in expression (Fig. 4c). The pairwise variations of the 18 candidate ER miRNAs under salt stress showed that V3/4 = 0.138 < 0.15, indicating that the optimal number of ER miRNAs for quantitative expression of miRNAs was three, namely mtr-miR166e-5p_1ss10GT, PC-3p-160_22001, and gma-miR403a_2ss19TC20TC (Fig. 4d).

NormFinder analysis

The analysis results from NormFinder show that under drought stress, 16 out of the 17 candidate ER miRNAs have a Stability Value (SV) of less than 1, indicating good stability. Among them, ptc-MIR156b-p3 exhibits the strongest stability. However, mdm-MIR482a-p3_2ss11CG20TA, and mtr-miR396b-3p_L-1 have SV values greater than 1, which indicates poor stability (Fig. 5a).

Figure 5. 

Normfinder stability values of 18 candidate reference miRNAs. (a) Normfinder stability values under drought stress. (b) Normfinder stability values under salt stress. The abscissa represents different Normfinder stability values, and the ordinate represents different candidate reference miRNAs name.

Under salt stress, 17 out of the 18 candidate ER miRNAs have an SV value of less than 1, indicating good stability. Among them, csi-miR1515a exhibits the strongest stability, while ptc-MIR156b-p3 has the poorest stability (Fig. 5b).

Bestkeeper analysis

The Bestkeeper software primarily evaluates the stability of ER genes by calculating the standard deviation (SD) and the coefficient of variation (CV). The analysis results show that under drought stress, 14 out of the 17 candidate ER genes, including U6, have an SD value greater than 1, indicating good stability. Among them, ptc-MIR156b-p3 exhibits the best stability with a coefficient of variation of 1.45% and a standard deviation of 0.37. The next most stable genes are mtr-miR156b-3p and PC-3p-160_22001. In contrast, mdm-MIR482a-p3_2ss11CG20TA has the poorest stability (Table 2).

Under salt stress, 11 out of the 18 candidate miRNAs have an SD value greater than 1. Among them, PC-3p-160_22001 exhibits the best stability with a coefficient of variation of 1.47% and a standard deviation of 0.36. In contrast, ptc-MIR156b-p3 has the poorest stability with a coefficient of variation of 7.26% and a standard deviation of 1.76 (Table 3).

RefFinder analysis

The RefFinder analysis provides a comprehensive evaluation and ranking of the stability of ER genes obtained from the GeNorm, NormFinder, and Bestkeeper software. The analysis results indicate that under drought stress (Table 4), ptc-MIR156b-p3 is the most stable, which is consistent with the results from NormFinder and Bestkeeper. The next most stable miRNAs are mtr-miR156b-3p and mtr-miR166a_L+1R-1, which align with the most stable miRNAs predicted by GeNorm. The least stable miRNA is mtr-miR396b-3p_L-1, which is consistent with the results from GeNorm and NormFinder, and it also has a low score in Bestkeeper. Overall, the high similarity among the results from the four software programs suggests that the findings are quite reliable.

Under salt stress (Table 5), csi-miR1515a was predicted to be the most stable miRNA, which is consistent with the results of the NormFinder analysis. In the Bestkeeper analysis and GeNorm prediction, the most stable gene PC-3p-160_22001 ranked second in the comprehensive scoring by RefFinder. Moreover, ptc-MIR156b-p3 was predicted to be the least stable gene by all four software programs, which also indicates a high degree of consistency in the results.

Previous studies have shown that drought and salt stress significantly alter the expression patterns of numerous genes and their products. Among these, miRNAs as key gene regulators often undergo substantial expression changes under abiotic stress^[26,27]. Quantitative real-time PCR (qRT-PCR) is widely used to analyze such expression changes due to its high sensitivity, specificity, and reproducibility^[28]. However, this high sensitivity also makes the technique susceptible to experimental variations, potentially compromising accuracy^[29]. Consequently, as miRNA research advances, the potential of using specific miRNAs as ERs has garnered increasing attention. In view of this, this experiment used qRT-PCR technology and five methods to conduct gene quantitative analysis on tree peony under drought and salt stress, aiming to screen out the most suitable ER miRNA.

The comparative ΔCt method identifies stable reference genes by analyzing the relative expression of 'gene pairs' across different samples or treatment conditions. The geNorm algorithm ranks candidate reference genes by calculating the average expression stability value (M), and the pairwise variation (V) between each gene and all others; as a result, it may recommend multiple optimal reference genes. In contrast, NormFinder not only assesses the overall expression variation of a gene but also identifies variation between sample subgroups, and estimates the systematic error that might be introduced by its use. BestKeeper determines the most suitable reference genes by analyzing the expression correlation among the candidate genes. RefFinder integrates the ranking results from the four methods mentioned above, assigns a specific weight to each, and ultimately provides a comprehensive stability ranking for the candidate reference genes^[30]. These five methods are all commonly used for reference gene selection and have been extensively reported in studies across various plant and animal species. Due to differences in their algorithms and analytical focus, the results they yield can vary. Therefore, employing a comprehensive analysis that combines all five methods leads to a more accurate and reliable identification of the most stable reference genes^[8,31].

MiR166 is a highly conserved miRNA family that plays crucial roles in plant growth, development, and stress responses^[32−34]. Its exceptional stability as a reference gene has been demonstrated across diverse species and conditions, including different tissues and developmental stages in Lilium henryi^[35] and nutrient-stressed grafted watermelon^[36]. These findings establish miR166 family genes as ideal reference candidates due to their high conservation and stability across plant species and environmental conditions.

Under salt stress, csi-miR1515a, PC-3p-160_22001, and mtr-miR166e-5p_1ss10GT were identified as the most stably expressed miRNAs. Among them, csi-miR1515a demonstrated high stability in this study, consistent with previous findings where miR1515 served as a reliable reference in leguminous plants^[37], and Gma-miR1515a normalized miRNA expression in salt-stressed soybean root tips^[38]. Current research also links miR-PC-3p to plant disease resistance, such as its regulatory role in maize infected by A. lepigone and in the pathogenicity of apple arc bacteria, highlighting its functional relevance beyond stress response.

U6 snRNA is widely employed as a reference for miRNA quantification due to its stable expression in species such as poplar, grapevine, and tea plant^[11,18,39]. Although included in this study as a comparative control, both previous screenings in tree peony and the present experimental results consistently demonstrate that U6 exhibits relatively low stability and is unsuitable as an ER under drought and salt stress conditions.

Analysis of candidate reference miRNAs revealed distinct stability patterns under drought vs salt stress. Under drought stress, ptc-MIR156b-p3, mtr-miR156b-3p, and mtr-miR166a_L+1R-1 exhibited the highest stability. The observed stability of ptc-MIR156b-p3 and mtr-miR156b-3p may be attributed to the high evolutionary conservation of the miR156 family, which targets SPL genes to regulate plant development^[40,41]. Under salt stress, however, csi-miR1515a, PC-3p-160_22001, and mtr-miR166e-5p_1ss10GT demonstrated superior stability. Notably, PC-5p-50_56138 maintained consistent performance across both stress conditions, ranking fourth in RefFinder analysis, suggesting its potential as a universal reference under osmotic stress, though further validation is warranted. These results emphasize the necessity of condition-specific reference gene selection.

The optimal reference miRNAs identified in this study for tree peony under drought and salt stress conditions (ptc-MIR156b-p3 and csi-miR1515a) belong to the same families as previously reported stress-related reference miRNAs in other plants. For example, gma-miR1515a, utilized as a reference in soybean salt-stress research^[38], and miR156, validated in poplar under pathogen stress^[42], both exhibit high expression stability under specific stress conditions, aligning with the fundamental principle that reference genes should be minimally affected by external perturbations. This study is the first to systematically validate the suitability of ptc-MIR156b-p3 and csi-miR1515a as reference miRNAs in tree peony under combined drought–salt stress, further supporting the notion that certain miRNA families retain stable expression under stress and are suitable as internal controls. Moreover, although earlier studies suggested that miR156b is involved in stress response^[43], and may therefore be unsuitable as a reference under stress conditions, integrating findings from other species with the present experimental results indicates that some members of the MIR156 family can maintain stable expression in stressful environments, demonstrating their potential as reliable reference miRNAs.

Based on the miRNA sequencing results, this study identified 17 miRNAs with relatively stable expression. Further stability analysis was conducted on samples subjected to different levels of drought and salt stress using qRT-PCR technology. The stability of 18 candidate ER genes, including U6, was assessed and ranked using three software programs: GeNorm, NormFinder, and Bestkeeper. Ultimately, the online tool RefFinder was employed to comprehensively analyze and evaluate the results from the aforementioned three software programs. The study found that under drought stress conditions, ptc-MIR156b-p3 is the most suitable miRNA to serve as an ER; under salt stress conditions, csi-miR1515a exhibits the highest stability. Additionally, PC-5p-50_56138 has a good ranking under both drought and salt stress. These results provide important references for the further identification of tree peony miRNA expression and offer more options for the use of ER genes (Fig. 6).

Figure 6. 

Conceptual figure. The yellow and blue arrows represent drought stress and salt stress, respectively. The right panels show the optimal reference miRNAs identified by Delta CT, GeNorm, NormFinder, and BestKeeper analyses, respectively. The bottom panel displays the top-ranked reference miRNA determined by a comprehensive analysis with RefFinder.