Plant materials and stress treatments
Seeds of L. × formolongi (‘Razan No.2’) and L. regale were stratified at 4°C for approximately one month and sown in a greenhouse. According to a previous study in our laboratory, seedlings of lily must progress through a vegetative development stage (rosette stage) to reach the adult stage, in which bolting is followed by reproductive development . The fresh leaves of two lilies were obtained from three individual plants at the rosette leaf stage, bolting stage and visible bud stage. All samples were immediately frozen in liquid nitrogen and stored at -80°C.
Different tissues, including fresh leaves, middle stems, bulb scales from the middle of the bulb, roots and young seeds, were also harvested from single plants in the adult stage after flowering with three biological replicates in the two lilies, followed by freezing and storage.
To investigate the stability of the selected candidate reference genes under abiotic conditions, seedlings of lilies were transplanted to pots in chambers at 5 weeks after sowing and grown under a 16 h light (32°C)/ 8 h dark (16°C) cycle with 70% humidity. Fresh leaves of the 2 kinds of lilies were collected at 0 weeks, 12 weeks and 22 weeks after heat treatment.
For biotic stress treatment, leaves were detached from the middle of the stems of L. × formolongi and L. regale at the flower bud stage. After cleaning with distilled water, these leaves were treated with 5 × 104 conidia•mL−1 of Botrytis elliptica via surface spraying  and placed in trays covered with preservative film under 90-100% humidity. The leaves were collected at 0, 6, 12, and 24 hours post inoculation (hpi) and prepared for miRNA extraction.
Small RNA extraction and cDNA synthesis
miRNA was extracted using miRcute miRNA Isolation kits (Tiangen Bio, Beijing, China) according to the manufacturer’s instructions. The enzymatic addition of a poly (A) tail was performed using an E. coli Poly (A) Polymerase kit (Invitrogen, USA) according to the manufacturer’s instructions after measuring the quality and concentration of non-poly (A) RNAs. Then, cDNA was obtained using a Quant Script RT Kit (Tiangen Bio, Beijing, China) and prepared for further qRT-PCR analysis.
Selection of candidate reference genes
Three classic housekeeping genes, 5S, 5.8S and U6, were selected, and their sequences were obtained from previous studies. Many miRNAs were identified in small RNA libraries generated in our laboratory during the vegetative development and initial flowering for L. × formolongi, and extensive sequencing data provided additional potential reference genes. Expression levels and differences in deep-sequencing counts in the form of copy numbers and p-values were also analysed in the libraries . Therefore, we selected mature miRNAs that displayed consistently high expression with low variance as candidate reference genes. miRNAs with copy numbers greater than the average (total copies/ (sample numbers × total miRNAs)) in random samples were defined as showing high expression levels. In addition to conserved miRNAs, several novel miRNAs were also identified, and the secondary structures of their precursors were predicted using RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi). The secondary structures had to satisfy the criteria for the annotation of miRNAs reported previously [5, 35].
Primer design and qRT-PCR analysis
Specific miRNA primers were designed using the Primer-BLAST tool of NCBI (http://www.ncbi.nlm.nih.gov/tools/prime r-blast/) based on the sequences of the mature miRNAs and their precursors in the small RNA libraries according to parameters reported in a previous study . The details of the primers are shown in Table 1.
qRT-PCR was carried out by using THUNDERBIRD SYBR qPCR Mix Without Rox (Toyobo, Shanghai, China) on a Bio-Rad CFX96 system (CFX96 Touch, BIO-RAD, USA). Each reaction volume of 20 μL contained 10 μL of THUNDERBIRD SYBR qPCR Mix, 4.0 μL of each diluted (1:9) primer (10 μM), and 2.0 μL of the diluted (1:9) cDNA template. A melting curve was generated for each primer pair over a temperature range of 50-59.4°C to inspect the single specific amplification products.
Standard curve analysis was conducted to determine the corresponding correlation coefficients (R2) and the efficiency of amplification (E). The cycling program for product amplification was as follows: 95°C for 30 s (hot-start activation), followed by 40 cycles of 95°C for 10 s (denaturation) and 50-59.4°C/ 72°C for 20 s (annealing/ extension). Each reaction was performed in triplicate.
Stability assessment of candidate reference genes
The stability of candidate reference genes was examined utilizing three Excel-based algorithms tools, geNorm , BestKeeper  and NormFinder . Based on Cq values representing the expression levels of each candidate, comparisons between the expression stability values (M) generated by these tools were conducted according to the underlying principles and calculations described in the manual. In addition, the simple delta Ct approach was utilized to compare ‘pairs of genes’ in samples to evaluate the expression stability of candidate reference genes [13, 39]. To determine the best reference genes, the final ranking of the set of candidates in these programs was integrated for weighted rank aggregation analysis employing the R package RankAggreg .