Screening and validation of optimal mirna reference genes in different developing stages and tissues of lilium henryi baker

ABSTRACT Dynamic miRNA detection using the qRT-PCR technique requires appropriate reference genes to ensure data reliability. Previous studies have screened internal reference genes in

plants during embryonic development and various stress treatment, involving relatively few tissues and organs. There is no relevant miRNA study in _Lilium henryi_ Baker and limited research

on the optimal miRNA reference genes in lilies, such as 5S, 18S, U6 and Actin. Twelve genes were selected as candidate reference genes whose expression stability was analyzed in petals at

different developmental stages and other tissues using various algorithms, such as geNorm, NormFinder, BestKeeper, and Delta CT. The results revealed that the optimal combination of

reference genes for _Lilium henryi_ Baker petals at different developmental stages was osa-miR166m and osa-miR166a-3p, while that for different tissues of _Lilium henryi_ Baker was

osa-miR166g-3p and osa-miR166a-3p.Four important genes related to growth and development regulation, namely, osa-miR156a, osa-miR395b, osa-miR396a-3p, and osa-miR396a-5p, were selected for

validation. The findings of the present study could contribute to future investigations onmiRNA expression and the related functions in _Lilium henryi_ Baker while providing important

references for the normalization of the miRNA expression in other varieties of lily. SIMILAR CONTENT BEING VIEWED BY OTHERS SELECTION AND VALIDATION OF REFERENCE GENES IN ALFALFA BASED ON

TRANSCRIPTOME SEQUENCE DATA Article Open access 21 February 2025 VALIDATION OF SUITABLE REFERENCE GENES BY VARIOUS ALGORITHMS FOR GENE EXPRESSION ANALYSIS IN _ISODON RUBESCENS_ UNDER

DIFFERENT ABIOTIC STRESSES Article Open access 15 November 2022 SCREENING AND VALIDATION OF OPTIMAL REAL-TIME PCR REFERENCE GENES FOR _ABELMOSCHUS MANIHOT_ Article Open access 01 April 2025

INTRODUCTION Lilies (_Lilium_ spp.) are considered among the most important floral crops worldwide1,2. Lilies are recognized for their diverse flower colors and forms and are, therefore,

widely used in garden beautification and cut flower markets3,4,5. The growth and development of lilies involve a series of biochemical processes in which miRNAs (microRNAs) play important

roles, for instance, in cell proliferation, cell differentiation, growth, development, and response to external environmental stimuli6,7,8. Plant miRNAs were identified in _Arabidopsis_ and

_Oryza sativa_ in 20029, and since then, hundreds of similar small RNA molecules, which are collectively referred to as miRNAs, have been identified in several model organisms and cells. The

above studies have revealed the complexity of gene regulatory networks in plants, due to which miRNAs have become a research hotspot in recent years. An unstable expression of a reference

gene could lead to biased outcomes10. Among the various detection methods available, qRT-PCR is the preferred method for gene detection owing to its high sensitivity, strong specificity, and

convenient operation11,12,13. Currently, qRT-PCR based on stem-loop and poly(A) methods, which use primers designed specifically for this purpose, is commonly used for miRNA detection. The

stem-loop method14 offers greater specificity and sensitivity compared to the poly(A) method15,16,17. The selection of miRNA-related reference genes varies among different cultivars, and

lilies encompass multiple varieties. For instance, in a study on _Lilium pumilum_ DC. Fisch. For somatic embryos, the reference genes selected were lpu-miR159a and the F-box family protein

(FP)18. In _Lilium_ × _formolongi_, 5S RNA and EF were selected as reference genes19. In the OT (Oriental × Trumpet) hybrid ‘Robina’, U4 was reported as the reference gene20. Among the

bicolor cultivars of Asiatic hybrid lily, namely, Lollypop, Cancún, Vermeer, Sugar Love, and Enjoys, as well asnthree full-color cultivars of the species, namely, Côte d’Azur, Vivaldi, and

Montreux, U6 was identified as the reference gene21. In _Lilium leichtlinii_ var. _maximowiczii_, 18S rRNA was reported as the reference gene5. _Lilium henryi_ Baker is a species endemic to

China and located primarily in regions such as Hubei and Guizhou. This species has characteristic orange flowers, which confer high ornamental value. The bulbs of _Lilium henryi_ Baker are

used for human consumption and also in traditional medicine22. The species exhibits strong adaptability23 and excellent overall resistance. It is the only wild species in the Lilium genus

with a vining growth habit, owing to which it is considered an important breeding material in horticulture. The species is at the state of NT (Near Threatened) according the IUCN standard

(http://www.iplant.cn/rep/prot/Lilium%20henryi). For the above reasons, a detailed investigation of _Lilium henryi_ Baker is of paramount importance. In the same context, the selection of

miRNA reference genes and the establishment of a qRT-PCR detection system for _Lilium henryi_ Baker are imperative. Our research group pioneered miRNA omics studies related to flower color

traits in _Lilium henryi_ Baker (yet to be published). A small RNA sequencing of _Lilium henryi_ Baker was conducted and revealed numerous miRNA sequences. Among those sequences, ten miRNAs

with relatively stable expression were screened out and used as candidate reference genes. In addition, two commonly used reference genes were also selected. The stability of the candidate

reference genes was analyzed using different software tools, such as geNorm24, NormFinder25, and BestKeeper26. The objective of this study was to identify suitable miRNA reference genes for

different growth and development stages and various tissues to provide appropriate internal reference choices for accurate quantification of miRNA gene expression levels during the growth

and development of _Lilium henryi_ Baker. MATERIALS AND METHODS EXPERIMENTAL MATERIALS _Lilium henryi_ Baker was used as the experimental material in the present study (Fig. 1), Fig. 1

generate using Illustrator 2021 software. The plants of this species were cultivated in a greenhouse at the Beijing Academy of Agriculture and Forestry Sciences. Samples were collected at

different stages of flower development (S1: green bud stage; S2: bud stage; S3: coloring stage; S4: early blooming stage; and S5: full blooming stage) and also from the various tissues

(petals at the full blooming stage, stem-root, stem, and leaves). Three biological replicates were performed for each group of samples. After collection, the samples were rapidly frozen in

liquid nitrogen and stored at − 80 °C until use. EXPERIMENTAL METHODS EXTRACTION AND REVERSE TRANSCRIPTION OF THE MIRNAS The required miRNAs (from the petals of _Lilium henryi_ Baker at the

five stages of development and also from the stem-root, stem, and leaf tissues of the species) were extracted using the miRcute miRNA Isolation Kit (Tiangen Biotech Co.). The integrity of

the extracted miRNA was assessed using agarose gel electrophoresis, and the purity and concentration of the miRNA were determined using a NanoDrop 2000. After quality control analysis, the

miRNAs were reverse transcribed using the miRNA 1st Strand cDNA Synthesis Kit (using the stem-loop method) from Nanjing Vazyme Biotech Co. Ltd., following the manufacturer’s instructions.

First-strand cDNA synthesis was performed in a reaction volume of 20 µL, which contained 0.5 µg of the miRNA to be reverse transcribed in an RNase-free centrifuge tube. In this process, the

genomic DNA was first removed from the miRNA sample by heating at 42 °C for 2 min. Afterward, the first-strand cDNA synthesis reaction was conducted using the following temperature program:

25 °C for 5 min, 50 °C for 15 min, and finally, 85 °C for 5 min. SCREENING OF THE CANDIDATE REFERENCE GENES AND PRIMER DESIGN Using our research group’s published small RNA sequencing data

for _Lilium henryi_ Baker (PRJNA1014586), ten miRNAs with relatively stable expression were selected as candidate reference genes. In addition, commonly used reference genes for plant miRNA

qPCR were also selected, along with four miRNAs for validation. The sequences of primers used for all of the above analyses are provided in Table 1. These sequences were synthesized at

Beijing Tsingke Biotech Co. Ltd. QRT-PCR ANALYSIS OF THE CANDIDATE REFERENCE GENES A 50 µL reaction mixture, which contained 5 µL of the cDNA template, 2 µL each (10 µM) of forward and

reverse primers, 25 µL of 2 × miRNA Universal SYBR qPCR Master Mix, and 16 µL of ddH2O, was prepared according to the instructions of the Vazyme miRNA Universal SYBR qPCR Master Mix Kit.

Each reaction was performed in triplicate. The PCR was conducted using the following program: initial denaturation at 95 °C for 5 min, followed by 40 cycles of denaturation at 95 °C for 10 

s, followed by annealing at 60 °C for 30 s, and a final extension step at 72 °C for 15 s. A melting curve analysis was conducted using the following temperature program: 95 °C for 15 s, 60 

°C for 60 s, and 95 °C for 15 s. Each reaction system was run in triplicate. The amplification efficiency of the primers was determined by generating a standard curve using a fivefold

dilution series of the cDNA template in water (v:v) –at 1:4, 1:24, 1:124, 1:624, and 1:3124. After each dilution, qRT-PCR was performed for each primer pair to obtain the Cq values and

generate the standard curve. The amplification efficiency (E) and the R2 value were calculated, and the correction equation was determined to be: E = (5 −1/slope − 1) × 100%27,28. DATA

PROCESSING AND ANALYSIS The data obtained was processed and analyzed using the geNorm, NormFinder, and BestKeeper software programs, the Delta CT method29, and the RefFinder online tool30.

These algorithms were subsequently used to analyze the expression stability of the candidate reference genes during the five different stages of petal growth and development and in the

various tissues of the species, including petals at the full blooming stage, stem-root, stem, and leaves. RESULTS ANALYSIS OF THE SPECIFICITY AND AMPLIFICATION EFFICIENCY OF THE PRIMERS

qRT-PCR analysis was performed using first-strand cDNA synthesized from petal, stem-root, stem, and leaf samples collected at five different stages of growth and development. The resulting

melting curves of all 12 candidate reference genes and those of the 4 validation genes presented a single distinct peak (Fig. 2), indicating the specificity of the primers. The replicates of

each sample exhibited good consistency, further confirming the specificity of the primers. The amplification efficiency of each primer was determined based on the standard curves (Table 2).

All primers exhibited good specificity and amplification efficiency and met the experimental standards for qRT- PCR. Therefore, these primers were used in subsequent experiments. EXPRESSION

ABUNDANCE OF THE CANDIDATE REFERENCE GENES The Cq values of the 12 candidate reference genes at different flowering stages and from various tissues of _Lilium henryi_ Baker were summarized

to evaluate the stability of their expression at the transcriptional level (Fig. 2). Specifically, 18S in the petals of _Lilium henryi_ Baker at different developmental stages presented the

lowest Cq value of 3.27, while novel_20 presented the highest Cq value of 31.37 (Fig. 3a). Among the different tissues of _Lilium henryi_ Baker, the lowest Cq value of 2.76 was again

obtained for 18S, while novel_69 presented the highest Cq value of 36.34 (Fig. 3b). STABILITY ANALYSIS OF THE CANDIDATE REFERENCE GENES THE GENORM ANALYSIS geNorm software was used to

caculate calculates the pairwise variation and average expression stability value (M value) among all the candidate genes. In each iteration, the gene with the highest M value is excluded

until the most stable gene combination is obtained.The results of the geNorm analysis revealed osa-miR166a-3p and osa-miR166m as the most stable miRNA reference genes in the petals at

different developmental stages of _Lilium henryi_ Baker as well as in the different tissues of _Lilium henryi_ Baker (Fig. 4). Among the petals at different developmental stages in _Lilium

henryi_ Baker, all had M- values significantly lower than 1.5, indicating relatively stable expression levels. The three most stable miRNA internal controls were osa-miR166a-3p, osa-miR166m,

and novel_69 (Fig. 4a). The stability of the candidate miRNA internal reference genes in the different tissues of _Lilium henryi_ Baker indicated relatively stable expression levels. The

three most stable miRNA internal controls were osa-miR166a-3p, osa-miR166m, and osa-miR166g-3p. However, the expression stabilities of U6 and 18S were the lowest among the 12 genes (Fig. 

4b). Paolacci et al.31 suggested using multiple internal reference genes for data correction and thereby obtaining accurate results. Therefore, it is considered necessary to combine pairwise

variation to predict the optimal number of internal reference genes by calculating the value of Vn/n + 1. Accordingly, when Vn/Vn + 1 < 0.15, the optimal number of reference genes is n.

When Vn/Vn + 1 > 0.15, the optimal number of reference genes is n + 132. In the present study, all the obtained Vn/Vn + 1 values were less than 0.15 (Fig. 5), of which 2 was the optimal

number of internal reference genes for this experiment. Accordingly, by considering both M and V values, the suitable internal reference genes for the petals at the different developmental

stages of _Lilium henryi_ Baker and also for the various tissues of the species were osa-miR166a-3p and osa-miR166m. NORMFINDER ANALYSIS NormFinder ranks the stability of genes based on

their expression stability values, with smaller stability values indicating a further stable gene expression25.The results of NormFinder analysis revealed that in petals at different

developmental stages of _Lilium henryi_ Baker (Fig. 6a). The most stable miRNA reference gene with a stability value of 0.1 was osa-miR166m, and osa-miR166a-3p was also relatively stable

with a stability value of 0.14. According to this analysis, the least stable reference genes were U6 and osa-miR159f. In the different tissues of _Lilium henryi_ Baker were used (Fig. 6b).

The most stable miRNA reference gene with a stability value of 0.13 was osa-miR166g-3p, and osa-miR166m was also relatively stable, with a stability value of 0.28. The highest expression

stability values were obtained for U6 and novel_69, indicating that their expression stability was the poorest among all the candidate genes, therefore, these genes were not suitable for use

as reference genes. BESTKEEPER ANALYSIS BestKeeper software primarily evaluates the stability of candidate reference genes by comparing the standard deviation (SD) and the coefficient of

variation (CV) of their CT values. Smaller CV and SD values indicate greater stability of the candidate reference gene and its suitability as an internal reference. In general, an SD value

of less than 1 is considered indicative of stable gene expression. BestKeeper analysis of the petals at different developmental stages revealed 8 genes with an SD < 1, while 4 genes had

an SD above 1.The 3 most stable miRNAs were osa-miR166a-3p, osa-miR166m, and novel_1, with osa-miR166a-3p exhibiting the smallest standard deviation and coefficient of variation values,

which indicated that it was the most stable internal reference gene. U6 exhibited the lowest stability, with the largest standard deviation and coefficient of variation (Fig. 7a). BestKeeper

analysis results of the various tissues revealed 7 genes with SDs < 1, and the smallest SD value of 0.2 was obtained for 18S. Only 5 genes had an SD above 1. The order of stability of

the miRNA reference g genes was as follows: 18S > novel_20 > osa-miR167d-5p > osa-miR166g-3p > osa-miR160a-5p > osa-miR166m > osa-miR166a-3p > osa-miR159f. > novel_1 

> novel_77 > novel_69 > U6. The smallest sum of thestandard deviation and coefficient of variation values was obtained for novel_20, while U6 exhibited the lowest stability (Fig. 

7b). DELTA CT ANALYSIS The Delta CT method was used to selects the optimal candidate internal reference gene based on the mean standard deviation (mean SD) value26. A lower mean SDobtained

for an internal reference gene indicates a higher expression stability. Delta CT analysis of petals at different developmental stages. The mean SD values of osa-miR166a-3p, osa-miR166m, and

novel_1 were relatively low, indicating their good expression stability across all the samples. The most stable internal reference gene was osa-miR166a-3p, while U6 presented the largest

mean SD value, indicating that U6 has the poorest stability (Fig. 8a). Among the different tissues of _Lilium henryi_ Baker, 18S, novel_20, and osa-miR167d-5p presented relatively small mean

SD values, indicating good expression stability across all the samples. U6 presented the largest mean SD, indicating that it had the poorest stability (Fig. 8b). REFFINDER ANALYSIS

According to the ranks assigned to the genes based on each of the four programs, RefFinder assigns an appropriate weight to each gene and calculates the geometric mean of these weights. A

resultant lower value indicates higher stability of the candidate internal reference gene under the given experimental conditions. Figure 9 shows the Venn diagram of the 5 most stable

internal reference genes identified based on the geNorm, NormFinder, BestKeeper, and Delta CT analyses. The results of the comprehensive analysis conducted using RefFinder revealed that in

_Lilium henryi_ Baker petals at different developmental stages (Fig. 10a). According to the results of the analysis based on the five evaluation methods, osa-miR166m, and osa-miR166a-3p

exhibited relatively stable expression, with corresponding stability values of 1.19 and 1.41, respectively. U6 exhibited the least stable expression. According to these results of the geNorm

analysis, the optimal number of internal reference genes to be introduced in the experiment was 2: osa-miR166m and osa-miR166a-3p. The different tissues of _Lilium henryi_ Baker were used

(Fig. 10b). According to the results of the analysis based on the five evaluation methods, osa-miR166g-3p and osa-miR166a-3p exhibited relatively stable expression, with corresponding

stability values of 1.41 and 2.74, respectively. U6 exhibited the least stable expression. According to the results of geNorm analysis, the optimal number of internal reference genes to be

introduced in the experiment was again 2: osa-miR166g-3p and osa-miR166a-3p. According to all the algorithms used in the present study, in _Lilium henryi_ Baker, the miR166 family members,

including osa-miR166a-3p, osa-miR166g-3p, and osa-miR166m, consistently ranked high in stability. Therefore, it was inferred that the miR166 family genes could serve as suitable candidates

for use as internal reference genes in both petals at the different developmental stages of _Lilium henryi_ Baker and in the different tissues of the species. VALIDATION OF GENES According

to previous studies, miR156 plays an important role in regulating flowering time, lateral root development, and anthocyanidin synthesis33,32,35. miR396 is also widely involved in the

regulation of growth and development in plant petals, leaves, and other parts36,37.Therefore, the expression levels of miR156 and miR396 were verified. Among the petals at the different

developmental stages of the species, when the most stable reference genes were used, osa-miR166m and osa-miR166a-3p, the expression patterns of the target gene osa-miR156a were consistent

(Fig. 11a). This gene exhibited the lowest expression level during the green bud stage, reached its maximum expression level during the coloring stage as the petals developed, then again

exhibited a decreased expression level during the initial opening stage, and finally, a slight increase. However, for the least stable reference gene, U6, the expression level of osa-miR156a

was quite low during the coloring, initial opening, and full bloom stages. Similarly, the expression pattern of the target gene osa-miR396a-3p was consistent (Fig. 11b) in the different

developmental stages of petals, with slight decreases during the bud stage, increases during the coloring stage, sharp decreases during the initial opening stage, and a maximum expression

level during the full bloom stage. However, when U6, the least stable gene, was used as the reference gene, the expression pattern of osa-miR396a-3p varied and was opposite to that observed

when osa-miR166m and osa-miR166a-3p were used as reference genes. The expression pattern of the target gene osa-miR396a-5p was consistent (Fig. 11c) among the different developmental stages

of petals, with the maximum expression occuring during the full bloom stage. Therefore, osa-miR156a and osa-miR396a may have regulatory effects on the growth and development of petals,

especially during the blooming stage. Due to the gradual formation of spots in the lower part of petals during development and because the main component of the spots is anthocyanins, some

studies suggest that osa-miR395b is a key gene related to anthocyanin synthesis38. Therefore, osa-miR166m and osa-miR166a-3p should be used as reference genes during the five stages of petal

development to determine the expression of osa-miR395b in the lower and upper parts of petals o further verify the reliability of osa-miR166m and osa-miR166a-3p as internal reference

genes.The expression trend of osa-miR395b was consistent when osa-miR166m and osa-miR166a-3p were used as internal reference genes. The expression level of osa-miR395b was extremely low in

the lower part of the bud stage (Fig. 11d) and then significantly increased during the coloring stage. However, in the upper part of the petals (Fig. 11e), the expression level of

osa-miR395b gradually increased during the green bud stage, bud stage, coloring stage, and early bloom stage, while the expression level was extremely low in the full blooming stage. This

further proves the accuracy of osa-miR166m and osa-miR166a-3p as internal reference genes, whether for the entire petal or for different parts of the petal. Among the different tissues, when

the most stable reference genes, osa-miR166g-3p and osa-miR166a-3p, were used, the expression pattern of the target gene osa-miR156a remained consistent (Fig. 11f), with the lowest

expression level in the petals observed during the full bloom stage and the highest expression level observed in the stem. Therefore, osa-miR156a may be involved in the regulation of stem

growth.However, when the least stable reference gene U6 was used, the expression pattern was different, with the highest expression level observed in the stem and the lowest expression level

observed in the petals at the full blooming stage. Similarly, the expression patterns of the target gene osa-miR396a-3p (Fig. 11g) in the petals and the stem were consistent, with

relatively higher expression levels observed, while lower expression levels of this gene were observed in the stem and leaves. However, when the least stable reference gene U6 was used, the

expression level of osa-miR396a-3p was extremely high in the stem. The expression pattern of the target gene osa-miR396a-5p (Fig. 11h) also remained consistent, with the highest expression

level observed in the stem, while relatively lower expression levels were observed in the stem, leaves, and petals at the full bloom stage. Therefore, osa-miR396a-3p and osa-miR396a-5p may

be involved in the growth and development of stem rooting. However, when the least stable reference gene U6 was used, significant differences in the expression of osa-miR396a-5p were

observed among the different tissues, with the expression level in the stem approaching 0. DISCUSSION The application of qRT-PCR for the quantitative analysis of the expression profiles of

key genes is a fundamental and widely employed strategy for deciphering the mechanisms of plant growth and development. The selection of reference genes, particularly the internal control

genes, has a crucial influence on the accuracy of gene expression profiles39,40. The screening and validation of miRNA reference genes have been reported for several plants, including

_Brassica napu_s41, grapevine (_Vitis vinifera_ L. ‘Muscat Hamburg’)42, and rice (_Oryza sativa_)43. In poplar, U6 was identified as the most suitable reference gene for miRNA qRT-PCR

experiments44. In sugarcane buds, miR171 and 18S rRNA were revealed as the most suitable reference genes for individual use16. In Chinese cedar (_Cryptomeria fortunei_), the most stable

reference genes for miRNA expression were identified as novel16, cln-miR6725, novel1, and U645. The expression levels of reference genes vary significantly among different plant varieties,

tissues, organs, or physiological states. For instance, in a study by Yu et al.46, the reference gene used in the salt stress response experiment in the soybean variety ‘Williams 82’ was

gma-miR1520d. However, Liu et al.47 suggested that in Soybean seeds (‘Williams 82’) under different abiotic stress conditions, the most suitable combination of reference miRNAs was miR166a

and miR167d in the leaf samples and miR171a and miR156a or miR167a and miR171a in the root samples. Related research has also been conducted in the field of lilies. Zhang et al.48 used 5S as

an internal reference gene in _Lilium_ × _formolongi_. Gao et al.49 used 18S as an internal reference gene in _Lilium regale_. Moreover, U4 and U6 have also been used as internal reference

genes in lilies50,51. Therefore, it was considered necessary to conduct a study on the screening of miRNA reference genes in _Lilium henryi_ Baker. In this context, the present study was

conducted to determine the suitable reference genes for _Lilium henryi_ Baker. A total of 12 candidate genes were selected, and their performance was evaluated in different samples (flower

petals at different development stages, stems, stem-roots, and leaves) of _Lilium henryi_ Baker. Four commonly used algorithms (geNorm, NormFinder, BestKeeper, and Delta CT) were applied to

analyze the expression stability of these 12 candidate reference genes in the petals at different developmental stages and in various tissues of the species. The objective was to evaluate

and determine stable reference genes. While the top five genes selected based on the results of the different algorithms were generally similar, slight differences were observed in their

stability levels, which was attributable to the different calculation methods employed by each software package. However, despite these differences, the results obtained using these

different algorithms were significantly consistent in terms of the selected best reference genes. According the NormFinder analysis, the stability values of osa miR166m and osa miR66a-3p

were both 0.14, indicating that osa-miR166m and osa-miR166a-3p exhibited the highest stability (Fig. 6a). To further analyze the results obtained using the four algorithms, an online

analysis tool referred to as RefFinder was used. The results obtained using the geNorm, NormFinder, BestKeeper, and Delta CT methods were consistent with the results obtained using

RefFinder, indicating the accuracy of the software analysis. In the case of different tissues, the results obtained using geNorm and NormFinder were generally consistent, while those

obtained using BestKeeper and Delta CT presented certain differences. Similar findings were reported in previous studies and therefore, these differences are acceptable from a comprehensive

perspective52,51,54. Further analysis using RefFinder revealed results consistent with those of geNorm and NormFinder analyses. The best combination of reference genes was determined based

on the comprehensive analysis conducted using RefFinder. The results indicated that suitable reference genes should be selected for different experimental samples, as not all commonly used

reference genes are applicable to different tissue samples or experimental conditions. This finding was similar to that reported for _Euscaphis konishii_ Hayata10. To further validate the

accuracy of the selected reference genes, osa-miR156a, osa-miR396a-3p, and osa-miR396a-5p were selected as validation genes, owing to their important regulatory roles in plant growth and

development. Moreover, during the development of lily petals, spots gradually form in the lower tissue. To investigate whether the selected internal reference genes involved in petal spot

formation are suitable for miRNA transcription level analysis, the expression level of the osa-miR395b gene, which is related to anthocyanin synthesis was detected. The results confirmed the

accuracy of the screened reference genes, indicating that the optimal internal reference genes osa-miR166m and osa-miR166a-3p can be used for analyzing miRNA transcription levels in

different tissues of petals at different developmental stages. Previous studies have demonstrated that U6 is highly conserved across different plant species, implying that its sequence is

relatively stable among different plants. Consequently, the small nuclear RNA of U6 is widely used as an internal reference gene for miRNA studies on various plants, including tea plants55,

rice (_Oryza sativa_)56, and tomato57. However, while U6 is commonly used as a reference gene in several plant species, differences in the plant genomic structures and regulatory mechanisms

may lead to variations in the expression of U6 under certain conditions. For instance, in _Lilium_ species (_Lilium pumilum_ DC. Fisch. and _Lilium davidii_ var. unicolor)51, U6 reportedly

exhibited relatively lower stability ranks, while in sweet potato (_Ipomoea batatas._ L.)58, U6 expression was relatively unstable. The present study revealed that the commonly used

reference gene U6 exhibited unstable expression during the growth and development of _Lilium henryi_ Baker. On the other hand, certain miRNAs from the miR166 family exhibited relatively

stable expression and were more suitable as reference genes. When miR166 is used as a reference gene, it typically includes members of its subfamily, such as miR166a, miR166b, miR166c, and

miR166m. The miR166 family is a group of similar miRNAs that share the same miRNA precursor sequence or have highly conserved mature miRNA sequences. These subfamily members frequently

exhibit similar expression patterns and functions. Therefore, when selecting reference genes, the entire miR166 family or its multiple subfamilies may be considered for detection. This

approach allows for a further comprehensive and accurate assessment of the expression levels of the miR166 family of genes used as reference genes. Similarly, in Tamarillo (_Solanum

betaceum_) callus tissue samples, miR166a was reportedly the most stable miRNA59. Previous studies have screened plant internal reference genes such as EF48, 18S49, U650 and Actin 51 during

embryonic development and various stress treatment, involving relatively few tissues and organs in lily. Therefore, this study screened and validated internal reference genes involved in

various flower development stages and tissue parts of _Lilium henryi_ Baker. A more precise validation was conducted on the selected internal reference genes in the petals, avoiding

interference from specific parts of the petals on the results of selecting the best internal reference genes. Therefore,the present study provides important information for studying the role

of miRNAs in the growth and development of lilies and offers a strategy and method for screening miRNA reference genes in other plants. However, further research is warranted to validate

and optimize the results obtained in the present study. Future studies could consider expanding the sample size and including samples from different developmental stages to further validate

the stability and applicability of the miR166 family. Further investigations based on reference gene screening could reveal the differential expression and functions of miRNAs during the

growth and development of lilies, providing a solid foundation and theoretical support for unraveling the mechanisms underlying the formation of important traits and the process of molecular

breeding in lilies. DATA AVAILABILITY The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. REFERENCES *

Alappat, B. & Alappat, J. Anthocyanin pigments: Beyond aesthetics. _Molecules._ 25, 5500 (2020). Article  CAS  PubMed  PubMed Central  Google Scholar  * Yamagishi, M. How genes paint

lily flowers: Regulation of colouration and pigmentation patterning. _Sci. Hortic-Amsterdam._ 163, 27–36 (2013). Article  CAS  Google Scholar  * Al-Ajlouni, M. G., Othman, Y. A., A’Saf, T.

S. & Ayad, J. Y. Lilium morphology, physiology, anatomy and postharvest flower quality in response to plant growth regulators. _S. Afr. J. Bot._ 156, 43–53 (2023). Article  CAS  Google

Scholar  * Bakhshaie, M., Khosravi, S., Azadi, P., Bagheri, H. & van Tuyl, J. M. Biotechnological advances in lilium. _Plant. Cell. Rep._ 35, 1799–1826 (2016). Article  CAS  PubMed 

Google Scholar  * Wang, Z., Li, X., Chen, M. M., Yang, L. Y. & Zhang, Y. C. Molecular and metabolic insights into anthocyanin biosynthesis for spot formation on lilium leichtlinii var.

maximowiczii flower petals. _Int. J. Mol. Sci._ 24, 1844 (2023). Article  CAS  PubMed  PubMed Central  Google Scholar  * Bajczyk, M. _et al._ Recent insights into plant mirna biogenesis:

multiple layers of mirna level regulation. _Plants-Basel._ 12, 342 (2023). Article  CAS  PubMed  PubMed Central  Google Scholar  * Wang, J. W., Czech, B. & Weigel, D. Mir156-regulated

spl transcription factors define an endogenous flowering pathway in arabidopsis thaliana. _Cell._ 138, 738–749 (2009). Article  CAS  PubMed  Google Scholar  * Willmann, M. R. & Poethig,

R. S. Conservation and evolution of mirna regulatory programs in plant development. _Curr. Opin. Plant. Biol._ 10, 503–511 (2007). Article  CAS  PubMed  PubMed Central  Google Scholar  *

Reinhart, B. J., Weinstein, E. G., Rhoades, M. W., Bartel, B. & Bartel, D. P. Micrornas in plants. _Gene. Dev._ 16, 1616–1626 (2002). Article  CAS  PubMed  PubMed Central  Google Scholar

  * Liang, W. X. _et al._ Selection and evaluation of reference genes for qrt-pcr analysis in euscaphis konishii hayata based on transcriptome data. _Plant Methods._ 14, 1–9 (2018). Article

  ADS  Google Scholar  * Gachon, C., Mingam, A. & Charrier, B. Real-time pcr: what relevance to plant studies?. _J. Exp. Bot._ 55, 1445–1454 (2004). Article  CAS  PubMed  Google Scholar

  * Kou, S. _et al._ Selection and validation of suitable reference genes for mirna expression normalization by quantitative Rt-Pcr in citrus somatic embryogenic and adult tissues. _Plant

Cell. Rep._ 31, 2151–2163 (2012). Article  CAS  PubMed  Google Scholar  * Kozera, B. & Rapacz, M. Reference genes in real-time pcr. _J. Appl. Genet._ 54, 391–406 (2013). Article  CAS 

PubMed  PubMed Central  Google Scholar  * Chen, C. F. _et al._ Real-time quantification of micrornas by stem-loop Rt-Pcr. _Nucl. Acids Res._ 33, e179 (2005). Article  PubMed  PubMed Central

  Google Scholar  * Shi, R. & Chiang, V. L. Facile means for quantifying microrna expression by real-time Pcr. _Biotechniques._ 39, 519–525 (2005). Article  CAS  PubMed  Google Scholar 

* Yang, Y. _et al._ Selection of reference genes for normalization of microrna expression by rt-qpcr in sugarcane buds under cold stress. _Front. Plant. Sci._ 7, 86 (2016). PubMed  PubMed

Central  Google Scholar  * Tashiro, R. M., Philips, J. G. & Winefield, C. S. Identification of suitable grapevine reference genes for Qrt-Pcr derived from heterologous species. _Mol

Genet Genomics._ 291, 483–492 (2016). Article  CAS  PubMed  Google Scholar  * Zhang, J. _et al._ Small Rna and transcriptome sequencing reveal a potential mirna-mediated interaction network

that functions during somatic embryogenesis in lilium pumilum Dc. Fisch.. _Front Plant. Sci._ 8, 566 (2017). Article  PubMed  PubMed Central  Google Scholar  * Zhang, Q., Zhao, Y., Gao, X.

& Jia, G. Analysis of mirna-mediated regulation of flowering induction in lilium × formolongi. _BMC Plant. Biol._ 21, 1–16 (2021). Google Scholar  * Chen, Y. _et al._ Key MicroRNAs and

target genes involved in regulating maturation in lilium. _Ornament. Plant. Res._ 2, 1–14 (2022). Article  Google Scholar  * Yamagishi, M. & Sakai, M. The Microrna828/Myb12 module

mediates bicolor pattern development in asiatic hybrid lily (lilium Spp) flowers. _Front Plant. Sci._ 11, 590791 (2020). Article  PubMed  PubMed Central  Google Scholar  * Yang, C., Shi, L.

& Lei, X. _Fanjingshan research editorial committee, edited by fanjing mountain medicinal plants editorial committee_ (Guizhou Science and Technology Press, 2016). Google Scholar  * Xue,

J. _King of bulb fresh cut flowers - Lily_ (Zhongyuan Farmers’ Publishing House, 2006). Google Scholar  * Hellemans, J., Mortier, G., De Paepe, A., Speleman, F. & Vandesompele, J. Qbase

relative quantification framework and software for management and automated analysis of real-time quantitative pcr data. _Genome Biol._ 8, 1–14 (2007). Article  Google Scholar  * Andersen,

C. L., Jensen, J. L. & Orntoft, T. F. Normalization of real-time quantitative reverse transcription-pcr data: A model-based variance estimation approach to identify genes suited for

normalization, applied to bladder and colon cancer data sets. _Cancer Res._ 64, 5245–5250 (2004). Article  CAS  PubMed  Google Scholar  * Pfaffl, M. W., Tichopad, A., Prgomet, C. &

Neuvians, T. P. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: bestkeeper–excel-based tool using pair-wise correlations. _Biotechnol.

Lett._ 26, 509–515 (2004). Article  CAS  PubMed  Google Scholar  * Dong, X. M., Zhang, W. & Zhang, S. B. Selection and validation of reference genes for quantitative real-time PCR

analysis of development and tissue-dependent flower color formation in _Cymbidium lowianum_. _Int. J. Mol. Sci._ 23, 738 (2022). * Tong, Z. G., Gao, Z.H., Wang, F., Zhou, J. & Zhang, Z.

Selection of reliable reference genes for gene expression studies in peach using real-time PCR. _BMC Mol. Biol._. 10, 71 (2009). * Xie, F., Xiao, P., Chen, D., Xu, L. & Zhang, B.

Mirdeepfinder: A mirna analysis tool for deep sequencing of plant small Rnas. _Plant Mol. Biol._ 80, 75–84 (2012). Article  CAS  Google Scholar  * Yao, L. _et al._ Selection of Q Rt-Pcr

reference genes for amomum tsaoko seeds during dormancy release. _China J. Chin. Mater. Med._ 46, 3832–3837 (2021). Google Scholar  * Paolacci, A. R., Tanzarella, O. A., Porceddu, E. &

Ciaffi, M. Identification and validation of reference genes for quantitative Rt-Pcr normalization in wheat. _BMC Mol. Biol._ 10, 1–27 (2009). Article  Google Scholar  * Vandesompele, J. _et

al._ Accurate normalization of real-time quantitative rt-pcr data by geometric averaging of multiple internal control genes. _Genome Biol._ 3, 34 (2002). Article  Google Scholar  * Asgari,

M., Mirzaie-asl, A., Abdollahi, M. R. & Khodaei, L. Flowering time regulation by the Mirna156 in the beet (beta vulgaris Ssp. maritima). _Plant Cell Tiss Org._ 150, 361–370 (2022).

Article  CAS  Google Scholar  * Yu, N., Niu, Q. W., Ng, K. H. & Chua, N. H. The role of Mir156/Spls modules in arabidopsis lateral root development. _Plant J._ 83, 673–685 (2015).

Article  CAS  PubMed  Google Scholar  * Gou, J., Felippes, F. F., Liu, C., Weigel, D. & Wang, J. Negative regulation of anthocyanin biosynthesis in arabidopsis by a Mir156-targeted spl

transcription factor. _Plant. Cell._ 23, 1512–1522 (2011). Article  CAS  PubMed  PubMed Central  Google Scholar  * Cao, D. Y. _et al._ regulations on growth and development in tomato

cotyledon, flower and fruit via destruction of mir396 with short tandem target mimic. _Plant. Sci._ 247, 1–12 (2016). Article  CAS  PubMed  Google Scholar  * Yang, F. X. _et al._

Organ-specific gene expression reveals the role of the cymbidium ensifolium-Mir396/growth-regulating factors module in flower development of the orchid plant cymbidium ensifolium. _Front.

Plant. Sci._ 12, 7997780 (2022). Article  Google Scholar  * Li, H. S. _et al._ Effects of blue light on flavonoid accumulation linked to the expression of Mir393, Mir394 and Mir395 in longan

embryogenic calli. _PLOS One._ 13, 191444 (2018). Google Scholar  * Derveaux, S., Vandesompele, J. & Hellemans, J. How to do successful gene expression analysis using real-time pcr.

_Methods._ 50, 227–230 (2010). Article  CAS  PubMed  Google Scholar  * Heid, C. A., Stevens, J., Livak, K. J. & Williams, P. M. Real time quantitative Pcr. _Genome. Res._ 6, 986–994

(1996). Article  CAS  PubMed  Google Scholar  * Machado, R. D. _et al._ Comprehensive selection of reference genes for quantitative gene expression analysis during seed development in

_Brassica napus_. _Plant Cell Rep._ 34, 1139–1149 (2015). Article  CAS  PubMed  Google Scholar  * Luo, M. _et al._ Selection of reference genes for Mirna Qrt-Pcr under abiotic stress in

grapevine. _Sci Rep-UK_ 8, 4444 (2018). Article  ADS  Google Scholar  * Verstraeten, B., De Smet, L., Kyndt, T. & De Meyer, T. Selection of mirna reference genes for plant defence

studies in rice (_Oryza Sativa_). _Planta._ 250, 2101–2110 (2019). Article  CAS  PubMed  Google Scholar  * Tang, F. _et al._ Selection and validation of reference genes for quantitative

expression analysis of mirnas and mrnas in poplar. _Plant Methods._ 15, 1–15 (2019). Article  ADS  Google Scholar  * Zhang, Y. T. _et al._ Selection and optimization of reference genes for

microrna expression normalization by Qrt-Pcr in Chinese Cedar (_Cryptomeria fortunei_) under multiple stresses. _Int. J. Mol. Sci._ 22, 7246 (2021). Article  CAS  PubMed  PubMed Central 

Google Scholar  * Yu, Y. H., Tang, J., Liu, C. & Ni, Z. Y. Microrna4359B positively regulates the soybean response to salt stress by targeting the F-Box protein Gmfbx193. _Environ. Exp.

Bot._ 206, 105177 (2023). Article  CAS  Google Scholar  * Liu, W. C. _et al._ Normalization for relative quantification of mrna and microrna in soybean exposed to various abiotic stresses.

_PLOS One._ 11, 155606 (2016). Google Scholar  * Zhang, Q., Zhao, Y., Gao, X. & Jia, G. Analysis of mirna-mediated regulation of flowering induction in lilium × formolongi. _BMC Plant.

Biol._ 21, 190 (2021). Article  CAS  PubMed  PubMed Central  Google Scholar  * Gao, X. _et al._ The Lre-Mir159a-Lrgamyb pathway mediates resistance to grey mould infection in lilium regale.

_Mol. Plant. Pathol._ 21, 749–760 (2020). Article  CAS  PubMed  PubMed Central  Google Scholar  * Yamagishi, M. & Sakai, M. The Microrna828/Myb12 module mediates bicolor pattern

development in asiatic hybrid lily (Lilium Spp.) flowers. _Front. Plant. Sci._ 11, 590791 (2020). Article  PubMed  PubMed Central  Google Scholar  * Zhang, J. _et al._ The use of mirnas as

reference genes for mirna expression normalization during lilium somatic embryogenesis by real-time reverse transcription pcr analysis. _Plant. Cell. Tiss. Org._ 129, 105–118 (2017). Article

  CAS  Google Scholar  * Barbosa Amorim, L. L. _et al._ Cowpea and abiotic stresses identification of reference genes for transcriptional profiling by Qpcr. _Plant. Methods._ 14, 1–17

(2018). Google Scholar  * Hong, S., Seo, P. J., Yang, M., Xiang, F. & Park, C. Exploring valid reference genes for gene expression studies in brachypodium distachyon by real-time Pcr.

_BMC Plant. Biol._ 8, 1–1 (2008). Article  Google Scholar  * Wan, H. _et al._ Selection of appropriate reference genes for gene expression studies by quantitative real-time polymerase chain

reaction in cucumber. _Anal Biochem._ 399, 257–261 (2010). Article  CAS  PubMed  Google Scholar  * Song, H. _et al._ Selection and verification of candidate reference genes for mature

microrna expression by quantitative Rt-Pcr in the tea plant (_Camellia Sinensis_). _Genes-Basel._ 7, 25 (2016). Article  PubMed  PubMed Central  Google Scholar  * Ding, Y., Chen, Z. &

Zhu, C. Microarray-based analysis of cadmium-responsive micrornas in rice (_Oryza Sativa_). _J. Exp. Bot._ 62, 3563–3573 (2011). Article  CAS  PubMed  PubMed Central  Google Scholar  * Koul,

A. _et al._ Carotenoid profiling, in silico analysis and transcript profiling of mirnas targeting carotenoid biosynthetic pathway genes in different developmental tissues of tomato. _Plant.

Physiol. Bioch._ 108, 412–421 (2016). Article  CAS  Google Scholar  * Liu, X. Y. _et al._ Optimization of reference genes for qrt-pcr analysis of microrna expression under abiotic stress

conditions in sweetpotato. _Plant. Physiol. Bioch._ 154, 379–386 (2020). Article  CAS  Google Scholar  * Cordeiro, D., Rito, M., Borges, F., Canhoto, J. & Correia, S. Selection and

validation of reference genes for Qpcr analysis of mirnas and their targets during somatic embryogenesis in Tamarillo (_Solanum Betaceumcav_.). _Plant. Cell. Tiss. Org._ 143, 109–120 (2020).

Article  Google Scholar  Download references FUNDING This research was funded by the China National Natural Science Foundation Projects for Youth (Grant number 32201606), Key research and

development projects of the National Forestry and Grassland Administration (Grant number GLM [2021] No. 72), Special projects for capacity building in scientific and technological innovation

of the Beijing Academy of Agriculture and Forestry (Grant number KJCX20230106),Beijing Innovation Consortium of Agriculture Research Sysem (Grant number BAIC09-2023). AUTHOR INFORMATION

Author notes * These authors contributed equally: Ge Jin and Xiuhai Zhang. AUTHORS AND AFFILIATIONS * Beijing Academy of Agriculture and Forestry Sciences, Key Laboratory of Urban

Agriculture (North China), Ministry of Agriculture and Rural Affairs, Beijing, 100097, China Ge Jin, Xiuhai Zhang, Shiyin Yu, Yunpeng Du & Mingfang Zhang * College of Forestry and

Grassland Science, Jilin Agricultural University, Changchun, 130118, China Ge Jin & Chunli Zhao * College of Landscape Achitecture, Beijing Forestry University, Beijing, 100097, China

Shiyin Yu & Meixian Wang Authors * Ge Jin View author publications You can also search for this author inPubMed Google Scholar * Xiuhai Zhang View author publications You can also search

for this author inPubMed Google Scholar * Shiyin Yu View author publications You can also search for this author inPubMed Google Scholar * Yunpeng Du View author publications You can also

search for this author inPubMed Google Scholar * Meixian Wang View author publications You can also search for this author inPubMed Google Scholar * Chunli Zhao View author publications You

can also search for this author inPubMed Google Scholar * Mingfang Zhang View author publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS All authors

contributed to the study conception and design. Material preparation was performed by G.J., S.Y. and Y.D. Data collection and analysis were performed by G.J. and M.Z. The first draft of the

manuscript was written by G.J. and M.Z. The reviewing and editing were performed by X.Z, M.W. and C.Z. The funding acquisition was performed by M.Z. All authors commented on previous

versions of the manuscript. All authors read and approved the final manuscript. CORRESPONDING AUTHORS Correspondence to Chunli Zhao or Mingfang Zhang. ETHICS DECLARATIONS COMPETING INTERESTS

The authors declare no competing interests. ADDITIONAL INFORMATION PUBLISHER'S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published maps and

institutional affiliations. RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing,

adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons

licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise

in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the

permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Reprints and

permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Jin, G., Zhang, X., Yu, S. _et al._ Screening and validation of optimal miRNA reference genes in different developing stages and tissues of

_Lilium henryi_ Baker. _Sci Rep_ 14, 1545 (2024). https://doi.org/10.1038/s41598-024-51562-1 Download citation * Received: 24 August 2023 * Accepted: 06 January 2024 * Published: 17 January

2024 * DOI: https://doi.org/10.1038/s41598-024-51562-1 SHARE THIS ARTICLE Anyone you share the following link with will be able to read this content: Get shareable link Sorry, a shareable

link is not currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative