The Asiatic lily cultivar ‘Tiny Padhye’ (a white-purple bicolor cultivar) was used in this study. Plants were grown in a greenhouse at the Chinese Academy of Agricultural Sciences (Beijing, China). No artificial light was provided. A total of 13 samples were collected. The expression stability of candidate reference genes was analyzed in the following four experimental sets. Samples from the first experimental set A represented inner tepals at four different developmental stages. Upper parts and bases of inner tepals were collected separately. Flower bud development stages (S) were as follows: stage 1 (S1; bud length of about 1.5 cm and no anthocyanin pigment is visible on tepals); stage 2 (S2; anthocyanin pigment becomes visible on tepal bases); stage 3 (S3; the day before anthesis, lower halves of tepals are fully pigmented); and stage 4 (S4; 0 day post-anthesis) (Figure 1). In experimental set B (shading treatment), 1.5-cm flower buds of Asiatic ‘Tiny Padhye’ were covered with silver paper, which was removed when the flowers developed to S3. Upper parts and bases of inner tepals were then immediately collected from these S3 inner tepals. Plants under natural light conditions were used as a control group. In experimental set C (plant tissues), upper parts and bases of S3 inner tepals, fresh bulb scales, stems, and leaves, were collected from untreated lily plants during flowering. The fourth experimental set D (all) was composed of all samples. The collected samples were frozen in liquid nitrogen immediately and stored at -80°C until use. Samples were obtained from 10 individual plants and three independent biological replicates were collected for each sample.

Total RNA was extracted using an RNAprep Pure Plant Kit (Tiangen, Beijing, China), according to the manufacturer’s instructions. RNA integrity was confirmed by denaturing 1.0% agarose gel electrophoresis. RNA quality and quantity were determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). Only the RNA samples with absorbance ratios at OD₂₆₀/OD₂₈₀ between 1.9 and 2.2 and OD₂₆₀/OD₂₃₀ greater than 2.0 were used for cDNA synthesis. First-strand cDNA synthesis was performed using a SuperScript III reverse transcription kit (Invitrogen, Carlsbad, CA, USA), according to the manufacturer’s instructions.

Ten reference genes (commonly used in qRT-PCR) and two target genes were selected from the transcriptome of the Asiatic lily cultivar ‘Tiny Padhye’ (Xu et al., 2016). Primers for ACT were previously described (Xu et al., 2016). Primer pairs for genes were designed using Primer Premier 5.0 software¹ with the following parameters: melting temperatures (T_m) of 55–65°C, primer lengths of 17–25 bp, and amplicon lengths of 100–300 bp (Table 1). Primer specificity was assessed using gel electrophoresis and melting-curve analyses.

qRT-PCR analysis was performed using SYBR^® Premix Ex Taq^TM II (Tli RNaseH Plus) (Takara, Dalian, China) and a Bio-Rad CFX96 system (Bio-Rad, Hercules, CA, USA) with the following reaction conditions: initial denaturation at 95°C for 30 s, 40 cycles of 95°C for 5 s and 60°C for 30 s, and a melting-curve program (65–95°C with a temperature increment of 0.5°C every 5 s). The melting curve was created to identify the amplicon specificity. All reactions were performed with three biological replicates and three technical replicates. The primer efficiency of each gene and regression coefficient (R²) were evaluated using a standard curve generated from a fivefold dilution series of cDNA (1/5, 1/25, 1/125, 1/625). A no-template control and a reverse transcription negative control were included to monitor the potential reagents and genomic DNA contamination, respectively.

To show the variation in the expression of each gene, boxplots of quantitative cycle (Cq) values for the 10 candidate reference genes were drawn using the boxplot R package. Expression stability of these 10 reference genes under different experimental conditions was analyzed using geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), BestKeeper (Pfaffl et al., 2004), and the comparative ΔCt method (Silver et al., 2006). All of these software tools were run in accordance with their manuals. Furthermore, comprehensive stability rankings were generated using the web-based tool RefFinder². It should be noted that the data of each biological replicate were analyzed separately.

Two lily anthocyanin biosynthesis-related genes, LhF3′H and LhUFGT, were selected as target genes to validate the reliability of the identified reference genes. Their gene expression profiles in the upper parts and bases of inner tepals at different developmental stages were normalized using the most stable candidate reference gene(s) as well as the least stable reference gene(s). Sample collections and experiments were performed as described above. The average Cq value was calculated from three biological replicates and three technical replicates and used for relative expression analyses. Relative gene expression levels of LhF3’H and LhUFGT were calculated using the 2^-ΔΔCq method (Livak and Schmittgen, 2001).

A total of 10 candidate reference genes were selected from the transcriptome of the Asiatic lily cultivar ‘Tiny Padhye’ for qRT-PCR analysis. The presence of a single PCR product with the expected size and single peak in the melting curve analysis for each gene confirmed the specific amplification of each gene (Figure 2). The amplification efficiency (E) of all PCR reactions ranged from 89.59% for ACT11 to 102.25% for 18S (Table 1), suggesting that these genes are suitable for further gene expression analysis. Meanwhile, the standard curves showed good linear relationships, with correlation coefficients (R²) above 0.99 (Table 1).

To assess the expression stability of the reference genes in different samples, the transcript abundances of 10 candidate reference genes were presented as their Cq values, which varied from 16.22 (GAPDH) to 32.76 (TUB), while the mean Cq values varied from 20.02 (GAPDH) to 28.39 (18S) (Figure 3).

The orders of candidate gene stability ranking under different experimental conditions were determined separately using geNorm, NormFinder, BestKeeper, the comparative ΔCt method, and RefFinder.

The geNorm program is used to rank the stability of expression of tested genes by calculating their expression stability values (M) based on the average pairwise expression ratio (Vandesompele et al., 2002). Candidate reference genes with lower average expression stability values (M) are considered to be more stably expressed reference genes (Vandesompele et al., 2002). The stability values (M) of tested genes evaluated by geNorm are shown in Figure 4. TIP41 and ACT were ranked as the two most stable genes for tepals at different developmental stages, and ACT11 and TIP41 for different tissues and all samples, while TUB and GAPDH were the two least stable genes. For tepals under shading treatment, CYP and F-box were the two most stable genes with the lowest M-value.

The pairwise variation (V_n/V_{n + 1}) value calculated by the geNorm algorithm determines the optimal number of reference genes for gene expression normalization. A value of V_n/V_{n + 1} < 0.15 suggests that the optimal number of reference genes equal to a value of n is sufficient (Vandesompele et al., 2002). In this study, the value of V₂/V₃ was below 0.15 for tepals at different developmental stages and under shading treatment (Figures 5A,B), indicating that two reference genes would be sufficient for gene normalization under these experimental conditions. However, when different tissues and total samples were considered, all V_n/V_{n + 1} values were still above 0.15 (Figures 5C,D), and thus the cut-off value of 0.15 was somewhat strict in these cases. Many studies have shown that the use of 0.15 as a cut-off value is just a recommendation, and whether 0.15 is actually used is dependent on the data (Silveira et al., 2009; Fernandez et al., 2011; Lu et al., 2013; De Lima et al., 2016). Generally, a reliable result could be obtained by using three reference genes in the majority of experiments (Toegel et al., 2007; Maroufi et al., 2010; Wang et al., 2014; Wang and Lu, 2015; Li M.Y. et al., 2016). Considering this, we opted to use three reference genes in these experimental conditions.

NormFinder ranks the stability of expression of candidate reference genes by calculating the average pairwise variation in a gene relative to that of other candidate genes (Andersen et al., 2004). Genes with the lowest expression stability value are characterized as the most stable reference genes (Andersen et al., 2004). The stability values of tested genes evaluated by NormFinder are listed in Table 2. For different tissues, the most stable gene was ACT, followed by ACT11 and EF1-α. For tepals under shading treatment, TIP41 and F-box were the two most stable genes. For tepals at different developmental stages, the three most stable genes were TIP41, ACT, and ACT11. Finally, ACT was the most stable reference gene while TUB was the least stable one for all samples.

BestKeeper is another software tool for determining the expression stability of reference genes by calculating the coefficient of variation (CV) and the standard deviation (SD) of Cq values (Pfaffl et al., 2004). The most stable reference genes are identified as those with the lowest CV and SD (Pfaffl et al., 2004). In this study, EF1-α and CYP were identified as the most stable genes for different tissues, ACT and ACT11 for tepals under shading treatment, ACT and TIP41 for tepals at different developmental stages, and EF1-α and ACT for all samples (Table 2).

The comparative ΔCt method identifies potential reference genes by comparing the relative expression of gene pairs within each sample (Silver et al., 2006). The ranking orders of the ten tested genes generated by the ΔCt method are shown in Table 2. ACT, TIP41, and ACT11 were the three most stable genes for tepals at different flower developmental stages, different tissues, and all samples. For tepals under shading treatment, the most stable gene was TIP41, followed by SAND and F-box (Table 2).

The ranking orders of the ten tested genes generated by geNorm, NormFinder, BestKeeper, and the comparative ΔCt method showed some differences (Table 2). To provide a comprehensive evaluation of candidate reference genes, further analysis was thus carried out using the web-based comprehensive tool RefFinder, which integrates geNorm, Normfinder, BestKeeper, and the comparative ΔCt method. The comprehensive ranking orders recommended by RefFinder are shown in Table 2. For all samples and tepals at different developmental stages, TIP41, ACT, and ACT11 were ranked as the top three most stable genes. TIP41 and F-box were found to be the most stable genes for tepals under shading treatment, while ACT, ACT11, and EF1-α were the most stable genes for different tissues.

The best combination of reference genes in different experimental sets was determined based on the optimal number of reference genes calculated by geNorm and the ranking list obtained using RefFinder. Specifically, TIP41 and ACT were found to be the best combination of reference genes for tepals at different developmental stages, TIP41 and F-box for tepals under shading treatment, ACT, ACT11, and EF1-α for different tissues, and ACT, TIP41, and ACT11 for all samples.

To validate the selected candidate reference genes, the expression levels of two anthocyanin biosynthesis-related genes (LhF3′H and LhUFGT) were investigated using different reference genes in inner tepals at different developmental stages. Each of the two most stable reference genes (ACT and TIP411), combination of stable genes (ACT + TIP41), and the least stable reference gene (TUB) were used as internal controls. When using ACT alone, TIP41 alone, and the combination of ACT + TIP41 for normalization, LhF3′H and LhUFGT showed the highest expression in tepal bases at S2 and an extremely low level of expression in upper tepals at all four stages (Figure 6). However, when the least stable gene TUB was used for normalization, the expression patterns were completely different. Namely, LhF3′H and LhUFGT had the highest expression levels in tepal bases at S4 (Figure 6).