Selection of Suitable Reference Genes for qPCR Normalization under Abiotic Stresses and Hormone Stimuli in Carrot Leaves

Carrot, a biennial herb of the Apiaceae family, is among the most important vegetable crops in the world. In this study, nine candidate reference genes (GAPDH, ACTIN, eIF-4α, PP2A, SAND, TIP41, UBQ, EF-1α, and TUB) were cloned from carrot. Carrot plants were subjected to abiotic stresses (heat, cold, salt, and drought) and hormone stimuli (gibberellin, salicylic acid, methyl jasmonate, and abscisic acid). The expression profiles of the candidate reference genes were evaluated in three technical and biological replicates. Real-time qPCR data analyses were performed using three commonly used Excel-based applets namely, BestKeeper, geNorm, and NormFinder. ACTIN and TUB were the most stable genes identified among all sample groups, but individual analysis revealed changes in their expression profiles. GAPDH displayed the maximum stability for most of single stresses. To further validate the suitability of the reference genes identified in this study, the expression profile of DcDREB-A1 gene (homolog of AtDREB-A1 gene of Arabidophsis) was studied in carrot. The appropriate reference genes were selected that showed stable expression under the different experimental conditions.


Introduction
Quantitative real-time reverse transcription polymerase chain reaction (qPCR) allows accurate high-throughput RNA quantification over a wide dynamic range at a relatively low cost; this technique has high sensitivity and has been widely used for gene expression analysis [1][2][3][4]. Appropriate reference genes could eliminate the discrepancy that may exist in different samples and ensure the accuracy and reliability of the experimental results. Discrepancies may be due to variations in RNA expression levels and the quality and efficiency of reverse transcription. The use of reference genes to measure the temporal and spatial expressions of the target gene is widely acknowledged as a standardized method. In higher plants, suitable internal controls for gene expression studies have been recognized for pepper [5], rice [6], Arabidopsis thaliana [7], Brachypodium distachyon [8], chicory [9], poplar [10], coffee [11], Oenanthe javanica (BI.) genes. The experimental data of the genes were determined by qPCR in carrot leaves under different hormone stimuli treatments (GA, SA, ABA, and MeJA, respectively) and abiotic stresses treatments (heat, cold, salt, and drought). All nine reference genes displayed a wide range of quantification cycle (Cq) values across experimental samples, indicating variable expression. Furthermore, the expression level of DcDREB-A1, the homolog of AtDREB-A1 (DREB, Dehydration responsive element binding factor) gene of Arabidopsis, was assessed using different reference genes to validate the selection of candidate reference genes. We assumed that the reference genes identified in current study would enable better normalization and quantification of transcript levels in future expression studies on carrot plants.

Plant materials and treatments
Seeds of D. carota variety of Kurodagosun were sown in plastic pots containing a soil/vermiculite mixture (1:1) [45][46][47] and grown in an artificial climate chamber programmed for 16 h/8 h at 25°C/16°C for day/night conditions at a light intensity of~300 μmol•m -2 •s -1 and relative humidity 60%. Healthy and vigorous eight-week-old seedlings were used for treatments. In drought experiment, soil were irrigated with 500 mL of 20% PEG 6000 for 2 h in each pots. In salt experiment, leaves were sprayed with 500 mL of 0.2 M NaCl for 2 h. Cold and heat treatments were performed by exposing eight-week-old seedlings to 4 and 40°C temperatures in light incubators for 2 h, respectively. For hormone treatments, leaves were sprayed with 500 mL of SA (1.4 mM) [37,48], MeJA (0.8 mM) [38], GA (1.4 mM) [39], and ABA (0.1 mM) [40] for 2 h, respectively. Plants were sprayed or irrigated only once. GA, SA, MeJA (containing 0.02% (v/v) absolute ethanol and 0.02% (v/v) Tween-20), and ABA were dissolved in distilled water [49][50][51][52][53]. The pH of GA, SA, MeJA, and ABA dilutions were 2.8, 2.8, 6.7, and 5.3, respectively. In all cases, pots were placed in light incubators under optimal conditions with constant light intensity, being processed at the same time as plants subjected to the different stress conditions. Three biological experimental replicates were collected from three seedling samples performed in different pots for each treatment. Leaves were collected from the eight-week-old seedlings subjected to all treatments. The samples were frozen in liquid nitrogen and stored at −80°C until further use.

Total RNA extraction and cDNA synthesis
Frozen carrot tissues were disrupted under liquid nitrogen conditions using mortar and pestle. Total RNA extraction was performed according to the manufacturer's protocol (Tiangen, Beijing, China). The concentration and purity of RNA samples were measured by NanoDrop ND1000 spectrophotometer, and cDNA synthesis was performed using an A 260 /A 280 ratio of 1.8 to 2.0 samples. The genetic integrity was evaluated by 1.5% agarose gel electrophoresis. cDNA was synthesized from approximately 1,000 ng total RNA using the PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa, Dalian, China). The cDNA was ten-fold diluted series (10×, 10 2 ×, 10 3 ×, 10 4 ×,10 5 ×, and 10 6 × dilutions) for determining the amplification efficiency (E) and correlation coefficient (R 2 ) analysis; and eighteen-fold diluted for conducting the qPCR analysis of elicitor treatments.

Selection of candidate reference genes and primer design
Nine genes, TIP41, TUB, eIF-4α, UBQ, SAND, GAPDH, EF-1α, PP2A, and ACTIN, were used to identify the most stable reference genes for qPCR expression analyses of target carrot genes. These genes have already been identified and have been commonly used as internal controls in previous studies [12,21,28,44]. For this study, the Arabidopsis genes were selected from the TAIR database (http://www.arabidopsis.org). Potential homologs of the nine reference genes were identified from the genome and transcriptome data sequences of carrot, which were sequenced and analysized by our group (CarrotDB: http://apiaceae.njau.edu.cn/carrot/) [54]. The potential homologs sequences were aligned and edited by using BioEdit Sequence Alignment v 7.0.9 software. Primers were designed using Primer 6.0 (Premier Biosoft International, Palo Alto, CA) and DNAMAN 6.0 (Lynnon Biosoft, USA) according to the manufacturer's instructions. The primers used in qPCR, as well as their melting temperatures (80°C to 85°C), primer lengths (22 bp to 25 bp), GC content (44% to 60%), and amplicon lengths (80 bp to 180 bp) are provided in Table 1. Cloning information is presented in Table A in S1 File. The specificity of the amplicons was verified by using a single band of expected size in 1.5% agarose gel following electrophoresis and by the presence of a single peak in the qPCR melting curve. The target amplicons were sequenced to confirm specificity of the PCR products.

Quantitative real-time PCR assay
The qPCR was designed according to the minimum information forpublication of quantitative real-time PCR experiment guidelines [55]. Reactions used SYBR Green I Mix (TaKaRa, Dalian, China) in a 20 μL reaction volume and were performed in a 96-well plate on MyiQ single color real-time PCR detection system (Bio-Rad, Hercules, USA). Reaction mixtures contained 10 μL SYBR Green I Mix, 2 μL diluted cDNA, ddH 2 O, and a final primer concentration of 0.4 μM.
The following amplification conditions were applied: an initial denaturation step of 95°C for 30 s; 40 cycles at 95°C for 5 s; and 60°C for 20 s. The final dissociation curve was obtained from 65°C to 95°C to verify primer specificity. Each assay included three technical and biological replicates, and a standard curve of six serial dilution points. The general quality assessment of the PCR results was based on the amplification and melting curve profiles of the samples in relation to the assay controls (non-template controls). Mean Cq values of the ten-fold dilution series were plotted against the logarithm of the pooled cDNA dilution factors. The Cq values and the following equation were used to determine efficiency (E) of each gene with the slope of a linear regression model: % E = (10 [−1/slope] -1) × 100% [56]. Amplification efficiencies were calculated from standard curves with satisfactory linear relationships (R 2 > 0.99). All PCR processes displayed efficiency between 90% and 110%.

Data analysis
Three different types of Microsoft Excel-based software, namely, geNorm [29], NormFinder [57], and BestKeeper [58], were used to rank the expression stability of reference genes across all experimental sets. These data were either used directly for stability calculations (BestKeeper analysis) or were converted into relative quantities and imported into the geNorm and Norm-Finder using the formula 2 -ΔCq , in which ΔCq = the corresponding Cq value-minimum Cq.
The raw data are listed in Table B in S1 File.
In geNorm, the reference gene expression stability measurement (M) value is calculated as the level of pairwise variation for each reference gene with all other control genes and as the standard deviation (SD) of the logarithmically transformed expression ratios [29]. The reference gene with the lowest M value is considered the most stable gene [59]. Similar to geNorm, the NormFinder program is another Visual Basic application tool for Microsoft Excel that is used to determine the expression stabilities of reference genes [12]. Misinterpretations caused by artificial selection of co-regulated genes are avoided with this program [57]. BestKeeper determines the most stably expressed genes based on the coefficient of correlation to the candidate reference gene's Cq values [58]. Genes with the lowest SD and CV values are the most stable [60].

Cq values of candidate reference genes in carrot
Based on primer sequences from Table A in S1 File, cDNA of nine genes were cloned and identified in carrot leaves based on the data of carrot genome and transcriptome sequences. The gene expression levels were determined as Cq values (Table B in S1 File), and the transcripts of the reference genes showed different levels of abundance (Fig. 1). Mean Cq values of the genes ranged from 24.49 (EF-1α) to 32.96 (TIP41), and the Cq values of all the tested samples were between 18.62 (EF-1α) and 38.01 (TIP41). Low Cq values corresponded to high levels of expression. EF-1α showed high expression level with low Cq value. TIP41 and SAND showed low expression levels with high Cq values (Fig. 2).

Determination of the optimal number of reference genes in carrot
The optimal number of reference genes required for normalization was determined with geNorm using pairwise variations (Vn/n + 1) between the sequential normalization factors (NFn and NFn + 1, n 2). A large variation between the sequential normalization factors indicates that the added gene has a significant effect and is preferred for inclusion and calculation of a reliable normalization factor [29]. As shown in Fig. 3, the third gene had no significant effect (V 2/3 , low value) in cold and drought conditions. Thus, two reference genes were sufficient for normalizing gene expression under the cold and drought conditions. With a threshold of 0.15, three genes were sufficient for normalizing gene expression under heat and GA stress conditions, five for SA stress and six for MeJA stress. None of the gene selected was found to be appropriate in salt stress condition in the current study.

Expression stability of candidate reference genes in carrot
Three different software programs were used to calculate the expression stability of the candidate reference genes: geNorm, NormFinder, and BestKeeper. Eight different treatment sets were sorted into three groups: "abiotic stress" (heat, cold, salt, and drought), "hormone stimuli" (SA, GA, ABA, and MeJA), and "total" (samples in all treatments). Accordingly, 11 evaluation patterns were generated for both single stress treatments and groups.
According to geNorm, in which the default limit comprised M values less than 1.5, and except for TIP41 under ABA stress, all the other reference genes performed well under individual stress conditions (Table C in S1 File). EF-1α and ACTIN were the two best genes among the nine reference genes in SA and salt stress treatments. However, in the MeJA treatment, EF-1α and UBQ were the two best reference genes. In NormFinder, TIP41 was the most stable gene among the nine candidate genes under salt and SA stress conditions. UBQ was the most stable gene under GA and cold stress conditions. The BestKeeper analysis showed that most of the nine candidate genes had satisfactory stability. TUB, GAPDH, and UBQ were ranked at the top positions in most of the single stress treatments by BestKeeper analysis. All candidate genes were confirmed to be stable under GA treatment in BestKeeper.
Recognizing the best reference gene was difficult because of the complexity of the groups. The results of the analysis of the three groups of samples are shown in Table 2. The nine candidate genes performed well by geNorm analysis. In the "abiotic stress" group, ACTIN and UBQ were the two most stable genes, and eIF-4α and GAPDH ranked top two in the "hormone stimuli" group. ACTIN and EF-1α were the two most stable genes in the "total" group. ACTIN, EF-1α, and GAPDH performed well in all three groups by geNorm analysis. eIF-4α was the most stable reference gene with the minimum value of 0.005 obtained by NormFinder in the "hormone stimuli" group. ACTIN was the most stable reference gene with the value of 0.012 and 0.015 obtained by NormFinder in "total" and "abiotic stress" groups, whereas it ranked fourth in "hormone stimuli" group. GAPDH performed well in "hormone stimuli" group by NormFinder analysis, while it ranked the last one in "abiotic stress" and "total" group. In all three groups, ACTIN, eIF-4α, and TIP41 performed well in terms of stability according to NormFinder analysis. In BestKeeper, EF-1α was the most stable reference gene in "abiotic stress" group, whereas it was the least stable in "hormone stimuli" and "total" groups; PP2A ranked first in "hormone stimuli" and "total" group and it ranked fifth in "abiotic stress" group. TUB and GAPDH were more stable than the other genes in all three groups according to BestKeeper. In both "abiotic stress" and "total" groups, ACTIN was ranked first according to geNorm and NormFinder; whereas ACTIN was ranked sixth and eighth by BestKeeper, respectively.

Reference gene validation
To validate the selection of candidate reference genes, the relative expression of DcDREB-A1 was calculated by using the selected reference genes (Fig. 4). In A. thaliana, the expression of the AtDREB-A1 gene was induced by abiotic stress treatments, including cold, drought, and heat [61][62][63]. In this study, the expression profiles of DcDREB-A1 in carrot under heat stress condition were assessed by using six candidate reference genes. When the two most stable reference genes ACTIN and TUB were used for normalization, the expression levels of DcDREB-A1 peaked at 1 h and subsequently decreased at 2 and 4 h (Fig. 4). When the less stable reference gene UBQ and EF-1α were used for normalization, similar expression patterns were generated. By contrast, when PP2A were used for normalization, the transcript levels and expression patterns differed from those obtained using ACTIN and other suitable reference genes. When normalization was conducted based on the least stable reference gene PP2A, the expression patterns of DcDREB-A1 peaked at 2 h and decreased at 4 h.

Discussion
qPCR is broadly accepted as a method with high sensitivity and specificity. Such method is used because of its repeated quantitative dynamic range and the high-throughput analysis of gene transcript levels. Accurate normalization of gene expression against an appropriate internal control is required for a valid qPCR analysis. Gene transcripts with invariant abundance under various environmental stimuli are essential reference points for accurate data analysis [7]. Thus, reference genes should be validated under certain experimental conditions and in different species [60,64]. Leaves serve important functions in process of photosynthesis. In the growth and development of carrot, the accumulated photosynthetic products were transported to tuberous roots. Furthermore, leaf is a vital organ of the response to the abiotic stress and hormone signal. In this study, we tested suitable reference genes for the expression of target genes in carrot leaves. Nine genes (GAPDH, ACTIN, eIF-4α, PP2A, SAND, TIP41, UBQ, EF-1α, and TUB) were selected as candidate reference genes for stable expression assessment tests in carrot leaves. All nine candidate genes were cloned from carrot in this study based on our transcriptome and genome database, CarrotDB [54]. A single peak in the melting curve analyses confirmed the primer pair that showed specificity. Plants were subjected to different hormone stimuli (GA, SA, MeJA, and ABA), abiotic stresses (heat, cold, salt, and drought), and efficacy dilutions (10×, 10 2 ×, 10 3 ×, 10 4 ×,10 5 ×, and10 6 × dilutions). The expression data were collected following qPCR amplification and detection. The curves showed a good linear relationship default limit with R 2 >0.99, and their amplification efficiencies ranged from 93.8% to 108.8% (Table 1). The primer pairs and amplification conditions were acceptable in qPCR-based quantification [55].
Most of the Cq values were lying between 18 and 35 across all tested samples, and the mean Cq values ranged from 24 to 33 [12,60,65,66]. Here, the Cq values of PP2A (Cq max -Cq min <10 cycles; SDs = 1.24) were distributed more centrally than those of the other candidate genes, whereas the Cq value of EF-1α showed the highest variation (Cq max -Cq min <13 cycles; SD = 1.99). TUB, eIF-4α, ACTIN, PP2A, GAPDH, and UBQ showed moderate expression levels. The mean values of all reference genes except EF-1α were close to the median values of these candidate genes, indicating that Cq values are evenly distributed.
We compared three different approaches implemented in the software programs, namely, geNorm, NormFinder, and BestKeeper. Each program differed in terms of the composition and ranking of the most stably expressed reference genes candidates; such differences may be caused by the variations between the approaches [3]. Comparison of the results obtained from the three software programs could reveal the most stable reference genes under specific experimental conditions. We also have determined the optimal number of reference genes required for accurate normalization. However, setting a cut-off value was used in some references but not a necessary criterion [29].
In heat stress, three genes were sufficient for normalizing gene expression under heat stress conditions calculated by geNorm (V 3/4 value = 0.15). The software program suggested the use of eIF-4α, SAND and TUB for normalization. SAND was ranked the best reference gene in NormFinder, while eighth in BestKeeper. EIF-4α was ranked fourth in NormFinder, while seventh in BestKeeper. TUB was ranked fifth in BestKeeper, sixth in NormFinder. EF-1a was ranked the best reference gene in BestKeeper, sixth in geNorm and the last one in NormFinder. ACTIN was ranked the third in BestKeeper, fourth in geNorm and fifth in NormFinder. We considered the rankings of three algorithms together, and recommended ACTIN and TUB combined with eIF-4α or SAND, as the best combination of stable reference genes for qPCR in the heat treatment.
In cold stress, the pairwise variation V 2/3 = 0.13, indicated that the addition of third gene had no significant effect for normalization. Two most stable genes ACTIN and UBQ can be used by geNorm analysis. ACTIN could be the preferred reference gene with the ranking of fifth and sixth by NormFinder and BestKeeper. UBQ was identified by NormFinder as the most stable reference gene and showed a variation in BestKeeper (seventh-ranked). BestKeeper ranked SAND as most stable, while ranked seventh by NormFinder and eighth by geNorm. BestKeeper ranked GAPDH at the fifth position and it was ranked fourth by geNorm and NormFinder. Based on these results, UBQ combined with ACTIN or GAPDH were recommended as the best combination of stable reference genes for normalization in cold treatment.
In drought stress, the pairwise variation V 2/3 = 0.12, indicated that two genes were sufficient for normalizing gene expression according to geNorm. Two most stable genes GAPDH and ACTIN could be used in qPCR by geNorm analysis. GAPDH could be the preferred reference gene with the ranking of third and ninth by BestKeeper and NormFinder, respectively. ACTIN was identified by NormFinder as the most stable reference gene and ranked the second place in BestKeeper. BestKeeper ranked TUB as most stable, and third-ranked by NormFinder, while ranked eighth by geNorm. Based on these results, ACTIN combined with TUB or GAPDH were recommended as the suitable combination of stable reference genes for normalization in drought treatment. Similarly, GAPDH combined with eIF-4α and UBQ would be sufficient for the GA treatment; the suitable combination of GAPDH, ACTIN, eIF-4α, TIP41, and EF-1α would be sufficient for the SA treatments, and a suitable combination of GAPDH, eIF-4α, PP2A, SAND, UBQ, and EF-1α would be sufficient for the MeJA treatment.
The results of BestKeeper analysis in the three groups, namely, "abiotic stress", "hormone stimuli", and "total", indicated that they did not perform well. In "abiotic stress" group, two most stable genes ACTIN and UBQ can be used in qPCR by geNorm analysis. ACTIN was identified by NormFinder as the most stable reference gene and ranked sixth in BestKeeper; UBQ could be the preferred reference gene with the ranking of second by geNorm and BestKeeper and eighth calculated by NormFinder; EF-1a was identified by BestKeeper as the most stable reference gene and ranked third and seventh in geNorm and NormFinder, respectively; Best-Keeper ranked TUB at the second place (SD of TUB, UBQ and GAPDH = 1.33), and sixthranked by NormFinder and geNorm. ACTIN, UBQ, EF-1a and TUB were chosen as the stable reference gene combination in "abiotic stress" group. Similarly, eIF-4α, GAPDH, ACTIN, and TUB were selected for "hormone stimuli" group; ACTIN and TUB could be chosen as reference genes for the "total" group.
In recent studies, UBQ showed stability in tomato [67] and A. thaliana [7], however, failed to perform satisfactorily in rice [59] and soybean [68]. GAPDH is among the best reference genes for measuring gene expression in many tissues [11,19,69]. ACTIN showed instability under numerous experimental conditions [70], but this gene is shown to be a suitable reference gene in developmental studies [68]. TUB also displayed a acceptably variable expression pattern and could be regarded as a commonly used reference gene in recent studies [7,60,66].
Plants respond to abiotic stress in their environments in developmental, physiological, and biochemical ways using a network of transcription factors [71,72]. AP2/ERF transcription factor (APETALA2/ethylene-responsive factor) is a large family of plant-specific transcription factors that activates the expression of abiotic stress-responsive genes via specific binding to the dehydration-responsive element and cis-acting element in their promoters [73][74][75]. These DREB homolog genes were induced by heat in many plants, for example Zea mays [76], Chinese cabbage [77], Arabidopsis thaliana [78], and so on. Previous studies have shown that heat shock transcription factors could be transcriptionally controlled by DREB and important for the establishment of thermotolerance [78,79]. Over-expression of several DREB transcription factors in transgenic plants could enhance tolerance to heat stress in plants [42]. We collected six reference genes to normalize the relative expression of DcDREB-A1 under heat stress condition. The expression level of the DcDREB-A1 gene was normalized by the most stable reference genes (ACTIN and TUB). The less stable reference genes (UBQ and EF-1α) showed similar expression patterns, but expression levels varied for these reference genes. When PP2A was used for normalization, the expression patterns and transcript levels obviously differed from those obtained by normalization against ACTIN and other suitable reference genes. Thus, use of an untested reference gene may reduce accuracy or produce misleading results.
To our knowledge, this study is the first systematic analysis for the selection of superior reference genes for qPCR in carrot leaves under different 'abiotic stress' (osmotic, salt, cold and heat), 'Hormone stimuli' (SA, GA, ABA, and MeJA), and 'Total' (samples in all treatments) conditions. The most stable reference gene were not the same ones depending on the stress. This study also proved that no single gene could express stably in all cell types and under all experimental conditions. Our shortlist may provide further supports to find putative candidate genes for future experiments that address other environmental variables as treatment factors in carrot.

Author Contributions
Conceived and designed the experiments: ASX CT. Performed the experiments: CT QJ GLW FW ZSX. Analyzed the data: CT QJ. Contributed reagents/materials/analysis tools: ASX. Wrote the paper: CT. Revised the paper: ASX CT QJ.