Evaluation of putative reference genes for quantitative real-time PCR normalization in Lilium regale during development and under stress

Plant materials and sampling

Bulbs of L. regale with a circumference of 120 mm were stored at 4 °C for 3 months, and then planted in pots (diameter 100 mm) filled with substrate (sterile turf: vermiculite: perlite = 1:1:1, v/v/v). They were grown in a growth chamber with 12 h light (200 μmol · m⁻² · s⁻¹)/12 h dark photoperiod and 25 °C day/22 °C night thermoperiod, at a relative humidity of 80%. The plants were daily watered to maintain the field capacity to the maximum. For quantification at different developmental stages, one leaf of each plant was sampled every month after germination. The sampling was conducted for six months from the fifth to tenth leaf below the apex, throughout the vegetative and reproductive stages of the plants before they entered senescence.

The plants generated from the bulbs were used for other treatments, which were conducted during the flower bud stage. For drought stress treatment, the water supply was withheld and leaf samples were collected at 2–4, 6 and 8 d during the treatment. For salinity stress treatment, the pots with plants were completely saturated with water containing 100 mM NaCl, followed by leaf sampling at 2, 6, 12, 24, 36 and 48 h during the treatment. For heat or cold stress treatment, the pots with plants were transferred to another growth chamber with a same photoperiod at 37 °C (heat stress) or 4 °C (cold stress). Leaves from the treated plants were collected at 2, 6, 12, 24, 36 and 48 h during the treatment. One leaf of each plant was sampled at every time point from the fifth to tenth leaf below the apex.

For biotic stress treatment, the fifth to tenth leaf below the apex of each plant were detached and surface-sterilized by 0.1% sodium hypochlorite for inoculation. Botrytis elliptica was cultured on Potato Dextrose Agar medium (PDA) under near-UV light for 5 days. The conidia were collected by successive gentle vortexing in Tween 20 solution (0.05% Tween 20 in sterile deionized water). After adjusted to 5 × 10⁴ conidia · mL⁻¹, the suspension was sprayed on the abaxial surface of the leaves. Mock inoculations with water were also performed as a control group. The petioles of inoculated leaves were covered with cotton infiltrated by fertilized water, and the leaves were incubated on wet filter papers in petri dishes, which were placed in a growth chamber with 100% relative humidity, 12 h light (200 μmol · m⁻² · s⁻¹)/12 h dark photoperiod and 25 °C day/22 °C night thermoperiod. Samples were collected at 2, 6, 12, 24, 36 and 48 hpi (hours post inoculation).

All the samples collected were immediately frozen in liquid nitrogen and stored at −80 °C. Two biological replications were designed for every developmental stage or every time point during the stress treatments, each replication consisting of three leaves from three individual plants as a sample pool. Twelve sample pools from each treatment (except 10 from drought treatment), namely 70 sample pools from 210 plants in total were used for study.

Total RNA extraction and cDNA synthesis

Total RNA from each sample was extracted using the EASYspin Plus Plant RNA kit (RN38, Aidlab Biotech, Beijing, China) according to the manufacturer’s protocol. Each RNA sample was treated with an RNase-free DNase I (Promega, Madison, WI, USA) at 37 °C for 30 min to eliminate the trace contaminants of genomic DNA. The concentration and integrity of the samples were examined by NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and 1.4% agarose gel electrophoresis. The first strand cDNA was synthesized based on 1 μg total RNA of each sample using the M-MLV reverse transcription system (Promega) according to the manufacturer’s protocol.

Primer design and PCR confirmation

Primers were designed using Beacon designer 7 (Premier Biosoft, Palo Alto, CA, USA) following the stringent criteria: PCR product size between 75–200 bp, primer length between 18–25 bp, optimal T_m between 55–65 °C, GC content between 40–60%. The theoretical T_m of the amplicon for each gene were predicted by Primer Premier 5 (Premier Biosoft), in which the salt concentration was set as 50 mM. In order to check the specificity and suitable reaction condition of the primers, they were initially tested in a general PCR reaction with gradient temperature (52–60 °C) using cDNA as a template. For further confirmation, PCR products were cloned into pGEM-T Easy vectors (Promega) respectively and then sequenced by Sangon Biotech (Shanghai, China).

qPCR analysis and efficiency

Real-time quantitative PCR was conducted with MiniOpticon Real-time PCR System (Bio-Rad, Hercules, CA, USA), using SYBR Premix Ex Taq (Takara Bio, Kusatsu, Shiga, Japan) to monitor dsDNA synthesis. It was carried out in a 20 μl mixture containing 2 μl cDNA template, 10 μl 2× SYBR Premix Ex Taq, and each primer in the optimized concentration shown in Table 1, by the following thermal cycles: 95 °C for 2 min, followed by 40 cycles of 95 °C for 5 s, the optimized annealing temperature shown in Table 1 for 10 s and 72 °C for 15 s. In order to ensure the specificity of each primer pair, melting curves were recorded after Cycle 40 by heating from 65–95 °C stepwise by 0.5 °C every 5 s. Each qPCR reaction was performed in technical triplicate.

The reaction condition including the best primer concentration and annealing temperature was optimized using the cDNA sample pool of the total samples used in this study as templates, which were analyzed as five serial ten-fold dilutions. The amplification efficiency of each primer pair was derived from a standard curve generated after the optimization. Mean quantification Cycle (Cq) values of each ten-fold dilution were plotted against the logarithm of the cDNA dilution factor (Ginzinger, 2002).

Data analysis

The Cq values for each tested reference gene were analyzed by a specific threshold. The raw Cq values were converted into relative quantities using the formula Q = E^−ΔCq, where E is the efficiency of the gene amplification for each primer pair and ΔCq is the lowest Cq value as calibrator (which corresponds to the sample with the highest expression) minus the Cq value of the sample tested. The data obtained for each primer were analyzed using geNorm V3.5 (Vandesompele et al., 2002), NormFinder (Andersen, Jensen & Orntoft, 2004) and BestKeeper (Pfaffl et al., 2004) for reference gene selection.

Normalization of LrLOX

The gene LrLOX (GenBank accession number KM051414), a putative homolog of Lipoxygenases (LOXs) in L. regale, was used to validate the selected reference genes. The primers were designed using Beacon designer 7 (Forward: 5′-TTCCAGCGACAACAGGAGCAC-3′; Reverse: 5′-CGTCGTCCACCAAATCCACTT-3′). The leaves of L. regale were inoculated with B. elliptica using the method mentioned above. The expression profiling of LrLOX in the samples collected at 2, 6, 12, 24, 36 and 48 hpi was analyzed under the same experimental procedure as above. The samples at 0 hpi were taken as control group. For comparative purposes, the relative expression of the target gene was calculated with different normalization factors based on the geometric mean of the three most stable genes and the most unstable genes under biotic stress, as well as the most stable genes in all conditions concluded from the analysis. The relative expression was calculated by standard E^−ΔΔCq method. The data collected were subjected to one-way ANOVA and the means were analyzed by Tukey’s range test. Data were represented as mean ± SD.

Selection of candidate reference genes and amplification specificity

In our previous study, a cDNA library of L. regale was constructed to analyze the gene expression pattern after inoculation with B. elliptica (Q Cui et al., 2015, unpublished data). Based on the RPKM (Reads per Kb per Million reads) values of the unigenes in the library (Table S1), as well as the previous studies on reference gene selection in other plants (Bin et al., 2012; Fu et al., 2012), a total of nine genes associated with a wide variety of biological functions were selected as candidate reference genes during development or under biotic and abiotic stress: 18S rRNA, ACT, BHLH, CLA, CYP, EF1, GAPDH, SAND and TIP41. Their full names, cellular functions and EST GenBank accession numbers are listed in Table 1. From at least three primer pairs per gene designed, the final primer pairs for the nine genes were selected on the basis of their single PCR products, amplification efficiencies and regression coefficients (R²), as shown in Table 2. The primers for each gene amplified a single PCR product at an expected size as confirmed in a 2% agarose gel electrophoresis (Fig. S1) and later sequencing (Seq. S1). Melting curve analysis also confirmed the single PCR products with no primer-dimers (Fig. S2), and their melting temperatures were consistent with the prediction (Table 2). Under optimized reaction conditions (Table 2), their R² values ranged from 0.992–0.999, and amplification efficiencies ranged from 93.0–105.4%.

Transcription profiling of the reference genes

The expression levels of the nine candidate reference genes were assessed for expression stability under diverse conditions including different developmental stages and biotic or abiotic stresses. The Cq (quantification cycle) values for the candidate reference genes in all the samples were used to compare the expression levels among these genes (Fig. 1). The data analysis showed that these reference genes were moderately abundant in all tested samples, exhibiting a wide range of expression levels. The mean Cq values of the reference genes varied from 18.88–24.40, with most lying between 20 and 23. GAPDH showed the highest expression level in all the samples with the lowest Cq value (18.89), while SAND showed the lowest expression level with the highest Cq value (24.41). CYP, ACT and GAPDH showed great variability with the Cq value ranging from 18.40–27.65, 19.90–26.85 and 16.37–23.29 respectively. Since the expression levels of the nine candidate reference genes showed a variation across all samples, further evaluation would be necessary to identify the best selection(s) from these candidate genes for normalizing gene expression under certain experimental conditions.

Gene expression stability analysis

Expression stabilities of the nine candidate reference genes were analyzed by geNorm and NormFinder. The geNorm program calculated average expression stability (M) for a reference gene based on the average pairwise variation between all genes tested (Vandesompele et al., 2002). The rankings of the nine candidate reference genes under different conditions according to the M values are depicted in Fig. 2. When all the samples were taken together, BHLH, SAND, CLA and EF1 (in order, similarly hereinafter) were the most stably expressed genes, while CYP and 18S rRNA were the least stable ones. In leaves across various developmental stages, EF1, SAND, TIP41 and GAPDH performed well, while CYP and BHLH were the least stable genes. CLA, SAND and EF1 had the most stable expression under biotic stress, whereas CYP and 18S rRNA were the most variably expressed of the nine candidate genes. Under drought stress, CYP, SAND, GAPDH and TIP41 were the most stable genes. By contrast, the most stable ones for salinity stress were CLA, TIP41, BHLH and GAPDH. In response to cold and heat stress, GAPDH, CYP and SAND were identified as the best two reference genes for both treatments.

The geNorm program was also used to determine the optimal number of reference genes required for the normalization across different sets of experiment by calculating pairwise variation (V_n/V_{n + 1}) value. According to the instruction, 0.15 was suggested to be the cut-off value for determining the optimal number of reference genes, below which the inclusion of additional reference genes is not required (Vandesompele et al., 2002). For various developmental stages, the V4/5 value was 0.145 as shown in Fig. 3, which means four reference genes, namely EF1, SAND, TIP41 and GAPDH should be considered. The experimental data for biotic stress (V2/3 value = 0.098) showed that CYP and SAND were sufficient for normalization. For drought stress (V2/3 value = 0.088) and salinity stress (V2/3 value = 0.051), the suitable reference genes were CYP/SAND and CLA/TIP41, respectively. In addition, the V2/3 value of cold stress was 0.059, indicating that these two reference genes, GAPDH and TIP41, were sufficient for normalization of gene expression. EF1 and GAPDH produced a V2/3 value of 0.085 in heat stress treatment, suggesting these two candidate genes could be selected as appropriate reference genes. However, the V_n/V_{n + 1} values were higher than 0.15 when taken all the samples together, and thus the cut-off value of 0.15 was a little strict in this case (Fernandez et al., 2011). The V5/6 value was 0.158, indicating that the optimal number of reference targets would be five: BHLH, SAND, CLA, EF1 and TIP41.

In general, similar results were obtained when the candidate reference genes were analyzed through NormFinder, which is based on variance estimation to calculate the Normalization Factor (NF) (Andersen, Jensen & Orntoft, 2004). These genes were ranked according to their stabilities under a given set of experimental treatments (Table 3). The ranking under biotic stress was similar as that of geNorm analysis. CYP and TIP41 were also recognized as the most stable genes under drought and salinity stress, respectively, while differences were observed in the analysis of genes under cold and heat stress, where SAND and CYP were ranked as the most stable one respectively instead. TIP41 and ACT were the most stable genes for various developmental stages. When all the samples were considered, TIP41, EF1, BHLH and SAND were the most stably expressed genes suggested by NormFinder, while CYP and 18S rRNA were suggested as the least stable genes by both the algorithms.

Their expression stability of expression was re-analyzed using the BestKeeper algorithm, which provides the Standard Deviation (SD) based on the Cq values of all candidate reference genes, as well as the correlation coefficients (r) of the genes with the BestKeeper Index calculated from the geometric mean of the remaining reference genes (Pfaffl et al., 2004). As shown in Table 4, the genes with a SD [± CP] value below 1 and a SD [± x-fold] value below 2 are considered to be qualified for the gene expression normalization, and their r values (correlation coefficient) are used for ranking. According to these criteria, EF1, BHLH and CLA were suggested as the three most stable genes when all the experimental conditions were taken into consideration. During different developmental stages, EF1, SAND and CLA were selected by the program for the most stable normalization. Almost all the reference genes were within acceptable criteria under various stresses, and the results were partially similar with those obtained by NormFinder.

Because of the diverse suggestions made by the three algorithms regarding the most stable reference genes, we used Venn charts to identify the most suitable ones under certain conditions by combining the top three genes for each algorithm together (Fig. 4). All the algorithms agreed that BHLH was the most stably expressed gene across all the experimental conditions, while CLA and EF1 were recommended by two of the three programs respectively. As a synthesis of the three programs’ output, EF1 under various developmental stages, SAND under biotic stress, CYP/GAPDH under drought stress, and TIP41 under salinity stress were generally considered suitable. All the three algorithms were consistently in agreement on the stability of SAND and GAPDH under cold stress, while only CYP was selected by all the three programs as for heat stress. Notably, ACT and 18S rRNA, the most commonly used reference genes, were hardly found to be among the most stable transcripts as suggested by our results.

Normalization of LrLOX

The expression level of LrLOX after inoculation with B. elliptica was calculated with the reference genes suggested as above. According to the result of the evaluation for biotic stress, SAND was the most stable reference genes, while CYP was the least stable one among the nine candidate genes. ACT and CLA were also recommended by two of the three algorithms, but ACT was considered to have a relatively lower stability than CLA and SAND. As shown in Fig. 5, when normalized with SAND, the transcription level of LrLOX was down-regulated (Fold Change (FC) = 0.65) at 2 hpi, and a maximum level (FC = 2.23) was obtained at 12 hpi (hours post inoculation) relative to uninoculated controls. Upregulation was still observed at 24 hpi, but the mRNA level was then steadily decreased in the next 24 h. Similar expression patterns were generated when normalized with CLA separately or in combination with SAND together. When normalized using ACT, the expression level of LrLOX at 24 hpi was shown as down-regulated (FC = 0.87), but it showed no significant difference if compared with the value normalized with SAND (Fig. 5). Normalization with a combination of the three most stable genes also made a consistent result. When normalized with BHLH, which was considered as the most stable gene overall but was ranked in the middle under biotic stress, the expression pattern of LrLOX was still roughly recognizable. However, significant difference could be observed at 12 hpi, as the maximum mRNA level was reduced to 1.70 fold. As expected, the expression profile showed more variation when normalized with CYP, the least stable genes calculated by the geNorm. Apart from the incorrect down regulation at 24 hpi here, the upregulation at 12 hpi was even as high as 6.66 fold.