Validation and Comparison of Reference Genes for qPCR Normalization of Celery (Apium graveolens) at Different Development Stages

A suitable reference gene is an important prerequisite for guarantying accurate and reliable results in qPCR analysis. Celery is one of the representative vegetable in Apiaceae and is widely cultivated and consumed in the world. However, no reports have been previously published concerning reference genes in celery. In this study, the expression stabilities of nine candidate reference genes in leaf blade and petiole at different development stages were evaluated using three statistics algorithms geNorm, NormFinder, and BestKeeper. Our results showed that TUB-B, TUB-A, and UBC were the most reference genes among all tested samples. GAPDH represented the maximum stability for most individual sample, while the UBQ displayed the minimum stability. To further validate the stability of reference genes, the expression pattern of AgAP2-2 was calculated by using the selected genes for normalization. In addition, the expression patterns of several development-related genes were studied using the selected reference gene. Our results will be beneficial for further studies on gene transcription in celery.


INTRODUCTION
Celery is widely cultivated as a vegetable crop in the world, and is rich in flavonoids, carotenoids, carbohydrate, and fibrin. During the plant development process, the physiological and morphological characteristics of celery have significant changes, which will affect the eating quality of leaf blade and petiole. In the last few years, high throughput sequencing technology has been widely used in celery research, for exploration gene transcriptional mechanism and regulation networks (Jiang et al., 2014a,b;Li et al., 2014a,b;Jia et al., 2015). These studies have identified many key genes that involved in lignin biosynthetic pathway during celery development stages (Jia et al., 2015). Another study indicated that some genes were associated with apigenin biosynthesis during celery leaf development (Yan et al., 2014). However, more regulatory networks during plant development are still waiting to be studied.
Gene expression pattern reflects the tendency of gene expression regulation, and provides a novel insight into understand the biological functions of genes. qPCR is a reliable and rapid method to evaluate the expression level of target gene, especially has very sensitive detection ability for some low copy mRNAs (Heid et al., 1996;Bustin, 2000;Mackay, 2004). To exclude the errors in mRNA extraction quality, reverse amplification efficiency, and qPCR procedures, the reference genes are needed for data correction and standardization (Radonić et al., 2004;Dheda et al., 2005). An ideal reference gene should have a constant expression in various tissues and different experimental conditions (Dheda et al., 2004). However, to date, no absolute reference gene has been identified in plants or animals (Volkov, 2003;Derveaux et al., 2010). Some studies directly selected the common reference genes such as ACTIN, GAPDH, and TUB to normalize the target genes without evaluating the expression stability. However, these reference genes have significant differences under different experimental conditions (Kim et al., 2003;Yan et al., 2012). The unstable expression of reference genes may cause the deviation of final result. Other researches pointed that two or more reference genes should be needed to normalize (Vandesompele et al., 2002;Schmid et al., 2003). Some valid statistical software have been developed, such as geNorm (Vandesompele et al., 2002), Bestkeeper (Andersen, 2004), NormFinder (Pfaffl et al., 2004), to evaluate the stability of the candidate reference genes under specific experimental conditions. Currently, several reliable reference genes have been reported in plants, and the stability of reference genes in different plant species are not completely consistent (Czechowski et al., 2005;Jiang et al., 2014a;Tian et al., 2015). ACT7 and PP2A genes displayed the maximum stability under abiotic stress conditions in Oenanthe javanica (Bl.) , ACTIN and TUB were the most stable genes in carrot (Tian et al., 2015). Moreover, the reference genes under different experimental conditions are also not the same. In the study of rice, eIF-4α and ACT1 were the most suitable reference genes during seed development (Li et al., 2010), UBQ5 and eEF-1α were most stable across all the tissue samples (Jain et al., 2006), while the 18S rRNA was the most reliable reference gene under various growth stages of etiolated seedlings and different cultivars (Kim et al., 2003). However, none of reference gene in celery has been reported. Hence, identification of suitable reference genes in various tissues and at different development stages will be required, which will greatly contribute accurate and reliable analysis of gene expression.
To accurately normalize the target gene expression in celery tissues and development stages, nine candidate reference genes were selected and their expression stability was evaluated. The target gene AgAP2-2, a gene encoding an AP2/ERF transcription factor, was used to validate the selection of reference gene. In addition, the expression patterns of development-related genes were also analyzed using the selected reference gene. This study aims to identify the most suitable reference genes that will provide a more accurate and reliable expression analysis of other celery genes among various tissues and development stages.

Plant Materials and Growth Conditions
The seeds of celery (Apium graveolens L. cv. Ventura) were cultivated in a controlled-environment growth chamber in Nanjing Agricultural University, China (32 • 02 ′ N, 118 • 50 ′ E). All plants were grown under a photoperiod of 16 h with 300 µ mol m −2 s −1 light intensity at 25 • C and 8 h dark condition at 16 • C. The relative humidity varied from 60 to 65%. Three development stages of celery were evaluated, and the height of the plant at Stage 1 was 10 cm (35 d), the height of the plant at Stages 2 was 20 cm (50 d), and the height of the plant at Stages 3 was 30 cm (65 d; Figure 1). Three biological replicate samples of celery leaf blade and petiole at each developmental stage were collected, then immediately frozen in liquid nitrogen and stored at −80 • C until use.

RNA Isolation and cDNA Synthesis
Total RNA was extracted using the total RNA kit (Tiangen, Beijing, China) and then treated with RNase-free DNase I (Takara, Dalian, China) to eliminate genomic DNA contamination. The quantity and quality of RNA samples were measured by agarose gel electrophoresis and the use of a Nanodrop ND 1000 spectrophotometer (Nanodrop Technologies Inc., Delaware, USA). Only the samples with an A 260 /A 280 ratio of 1.8-2.2 and an A 260 /A 230 ratio >1.8 were used for further analysis. Total RNA (1.0 µg) was reverse-transcribed into cDNA using a PrimeScript RT reagent kit (Takara, Dalian, China). The cDNA was efficacy dilutions (10X, 10 2 X, 10 3 X, FIGURE 1 | Growth status of celery at three developmental stages. The leaf blades and petioles at different developmental stages were presented, respectively. Stage 1, 35 days after sowing; Stage 2, 50 days after sowing; Stage 3, 60 days after sowing. Frontiers in Plant Science | www.frontiersin.org 10 4 X, 10 5 X dilution) for detecting the amplification efficiency (E) and correlation coefficient (R 2 ), and 18-fold dilution for qPCR experiments.

Selection of Candidate Reference Genes and Primer Design
Nine typical candidate reference genes of celery, ACTIN, EF-1α, GAPDH, RAP2, TBP, TUB-A, UBC, TUB-B, and UBQ, were selected for qPCR analysis in this study. Arabidopsis reference genes were downloaded from the TAIR database (http://www.arabidopsis.org) and used as query sequences to retrieve the homologs genes in celery. Based on the A. graveolens transcriptome sequencing data built by our group Jia et al., 2015), nine potential genes were cloned. We have submitted all the nucleotide sequences to GenBank, and the corresponding accession numbers were KU234487 (ACTIN), KU234488 (EF-1α), KU234489 (GAPDH), KU234490 (RAP2), KU234491 (TBP), KU234492 (TUB-A), KU234493 (UBC), KU234494 (TUB-B), and KU234495 (UBQ). The primers used for qPCR were designed by Primer Premier 6 using a standard set of design criteria (annealing temperatures 58-60 • C, primer lengths 18-26 bp, GC content between 40 and 60%, and the PCR product between 60 and 150 bp; Udvardi et al., 2008). Primer sequences and amplicon characteristics are listed in Table 1.

qPCR and Statistical Analysis
The qPCR reactions were performed in a 96 well plate using the MyiQ Single Color Real-Time PCR Detection System (Biorad, Hercules, CA, USA). Each 20 µL PCR reactions contained 2.0 µL of diluted cDNA, 0.4 µL of each primer (10 mM), 10 µL of SYBR Green I mix (Takara, Dalian, China), and 7.2 µL of ddH 2 O. The PCR conditions were as follows: at 95 • C for 30 s for predenaturation, 40 cycles of 95 • C for 5 s for denaturation, and 60 • C for 30 s for annealing and extension. A melting curve (65-95 • C, at increments of 0.5 • C) was generated to verify the specificity of primer amplification. Each PCR reaction was repeated three times, and three biological replicates were analyzed. In Addition, each assay contained a standard curve of different dilutions of the template and a no-template control. Amplification efficiency of each primer pair was calculated (%E = (−1 + 10 [−1/slope] ) × 100%) and correlation coefficient (R 2 ) was tested.
Expression level of nine genes in each reaction was determined by the cycle threshold Cq (the cycle at which a threshold fluorescence was obtained). The original data was presented in the Table S1. Three different Microsoft Excel-based softwares, geNorm (Vandesompele et al., 2002), NormFinder (Andersen, 2004), and BestKeeper (Pfaffl et al., 2004), were used to determine the best reference genes. These raw data can be directly used for BestKeeper program, but for geNorm and NormFinder, Cqvalues were converted into relative quantity values by the formula 2 − Cq , Cq = the corresponding Cq-value -minimum Cq.
(1) geNorm. In geNorm, the calculation principle relies on the expression ratio of two ideal internal control genes is identical in all samples. The gene expression stability measure M is calculated as the level of pairwise variation for that gene compared with all other tested reference genes. geNorm identify these genes by progressively eliminating less stable genes from the analysis, and the reference gene with the lowest pairwise variation is the most stable. Besides, optimal number of multiple reference genes was determined by pairwise variation (V n/n+1 ) between normalization factors (NF n and NF n+1 , n ≥ 2) in geNorm.
(2) NormFinder. The NormFinder program ranks all candidate reference genes on the basis of intragroup and intergroup expression variations, and then combines them into a stability value for each candidate reference gene. This program can avoid the misinterpretations which caused by artificial selection of co-regulated genes. The reference gene with the lowest stability value is the most stable.

Selection and Expression Analysis of Development Related Genes
The transcriptome sequencing of three celery development stages were finished by our group Jia et al., 2015). Base on the annotation of celery genes, some of the genes were related to plant growth and development. Transcript abundances were estimated by calculating read density as "reads per kilobase of exon model per million mapped reads" (RPKM; Mortazavi et al., 2008). The expression clusters of genes were analyzed by using Cluster (http://bonsai.hgc.jp/~mdehoon/software/cluster/ software.htm), and the heatmap was drawn using Tree View (http://jtreeview.sourceforge.net/). qPCR was also performed to analyze the gene expression. Each reaction had three technical and biological repeats. The relative gene expression levels were calculated with the 2 − CT method (Pfaffl, 2001). The gene-specific primers are shown in Table S2.

Amplification Specificity and Efficiency of Candidate Reference Genes
The specific primers of nine candidate reference genes were designed for qPCR, with the amplicon length ranging from 105 to 159 bp. A single peak in the melting curve showed the expected amplification effect (Figure 2). The correlation coefficients (R 2 ) and PCR amplification efficiencies of nine genes in leaf blade and petiole were calculated, respectively, the results met the standard (R 2 > 0.99, 90 < E% < 110; Ramakers et al., 2003; Figures S1, S2).
The Cq-values in qPCR provided an overview of the gene expression levels of nine candidate reference genes in test samples. The transcript abundant for each gene has significant difference and the raw data were listed in Table S2. In general, the Cq-values between 18 and 30 are considered to be appropriate and effective data (Czechowski et al., 2005;Derveaux et al., 2010;Jiang et al., 2014a). In this study, all the Cq-values were between 17.03 and 29.84, and the mean Cq-values of the genes ranged from 21.57 for TUB-A to 26.18 for TBP (Figure 3 and Table S3).
High Cq-value indicates the low expression levels, conversely mean the high expression. Among nine genes, GAPHD and TUB-A showed high expression with low Cq-values, whereas UBC and TBP showed low expression. The significant variation in gene expression indicated screening the appropriate reference gene should measure the stability.

Stability of Candidate Reference Genes Under Different Development Stages and Tissues
To select the most suitable reference gene, three methods (geNorm, NormFinder, and BestKeeper) were used to analyze the stability of each reference gene. The stability ranking of candidate reference genes in six individual samples were FIGURE 2 | Melting curves generated for nine candidate reference genes by qPCR.
Frontiers in Plant Science | www.frontiersin.org calculated, respectively ( Table 2). In addition, we divided these six samples into three groups: "Leaf blade, " "Petiole, " and "Total." For complexity of the groups, the ranks of the nine genes were calculated again and the results were shown in Table 3.
In geNorm analysis, the M-value was used to represent the average expression stability, the lower the M-value indicates a higher stability. In all six samples, the M-value of candidate reference genes were less than the default limit of 1.5, but the most suitable reference gene was different in different tissues and development stages. GAPDH was the most stable gene at Stages 1 and 2 in leaf blade, and was the most stable gene at Stages 2 and 3 in petiole. UBC exhibited relatively stable expression at Stage 3 in leaf blade and Stage 1 in petiole. Under all tissue sets, TUB-B was the most stable gene among the nine candidate reference genes in "Leaf blade, " "Petiole, " and "Total, " whereas UBQ was the least stable gene with the largest M-value in "Leaf blade" and "Total." Another method, NormFinder, also classified TUB-B as the most stable reference gene with the minimum value of 0.005 at Stage 3 in leaf blade. GAPDH showed good stability at Stage 1 in leaf blade and at Stage 3 in petiole. However, GAPDH showed the worst stability at Stage 1 in petiole. In three groups, TUB-B ranked first in "Leaf blade" and "Total" with the value of 0.051 and 0.117, and ranked second in "Petiole" with the value of 0.163. Moreover, UBQ was the most stable gene in "Petiole, " but was the least stable gene in "Leaf blade" and "Total." Based on calculations by BestKeeper software, the smaller SDand CV-value means the gene is more stable. BestKeeper ranked TUB-A, RAP2, UBC, UBQ, EF-1α, GAPDH, respectively, as the best reference gene under three development stages in leaf blade and petiole. Although the best reference gene in six individual samples was not the same, UBC had the lowest SD-and CVvalues in three groups. That meant UBC was more stable than the other genes in three groups "Leaf blade, " "Petiole, " and "Total." At the same time, we also found that UBQ had poor performance according to the ranking by BestKeeper.

The Optimal Number of Reference Genes for Normalization in Celery
The geNorm algorithm was used to determine the optimal number of reference genes by evaluating pairwise variation (V n/n+1 ) between normalization factors (NF n and NF n+1 , n ≥ 2). The V n/n+1 -value below 0.15 suggested that an additional reference gene is not necessary for normalization (Vandesompele et al., 2002). Among the nine reference genes, the most stable reference genes varied in different samples and groups (Figure 4). For leaf blade development, two reference genes were enough for normalization at Stages 1 and 2, while four reference genes were needed at Stage 3 with the V 4/5 -value dropping to 0.15. At three petiole development stages, four, seven, and two reference genes were needed, respectively. For group "Leaf blade, " five stable reference genes were proposed to be used. When the samples were analyzed as "Petiole" or "Total, " it seems that all the V-values were higher than 0.15.

Validation of the Selected Reference Genes
To validate the suitability of reference genes, five reference genes including two most stable reference genes TUB-B and TUB-A, two less stable reference genes UBC and RAP2, and the least stable reference gene UBQ were selected as calibrator. The relative expression levels of AgAP2-2 were, respectively, calculated by using the selected reference genes. The homologous gene in Arabidopsis, AtAP2 (AT4G36920), which have been confirmed to play an important role in plant development (Jofuku et al., 1994;Liu et al., 2014). In our transcriptome data (Jia et al., 2015), the transcript abundances of AgAP2-2 were significantly different at three developmental stages. As illustrated in Figure 5, the expression patterns of AgAP2-2 had a greater difference using different reference genes. When using the most stable reference genes TUB-B and TUB-A, the expression patterns of AgAP2-2 were consistent, and the expression level was higher in petiole, especially at Stage 3. Similar expression patterns were generated by using the less stable reference genes UBC and RAP2. However, when the least stable reference gene UBQ was used, the expression level of AgAP2-2 had a strong bias compared with other genes. This result demonstrated that the reference gene with stable expression was necessary to accurately normalize the expression of target gene.

Expression Abundances of Development-Related Genes in Three Development Stages
Many genes have been identified to play key roles in plant growth and development, such as transcription factor genes, hormonerelated genes, and genes encoding ribosomal proteins (Katagiri, 1992;Horiguchi et al., 2012;Saini et al., 2013). In Arabidopsis, lots of genes were identified to involve in development and multiple biological processes (Jofuku et al., 1994;Horiguchi et al., 2011;Köllmer et al., 2011). The high-throughput transcriptome sequencing of three celery development stages were finished by our lab and the annotation of celery genes showed many genes were involved in plant development (Jia et al., 2015). Basing on the transcriptome data in celery and   the orthologous genes between celery and Arabidopsis, we selected 15 genes which, respectively, belonged to AP2/ERF transcription factor family, auxin-related gene, and ribosomal protein gene for further study. The expression abundances of these genes were analyzed by calculating the RPKMvalues. As showed in Figure 6, the auxin-related genes (AgWAT1, AgABCB19, AgSNX1, AgAILP1, AgARF1) and ribosomal protein FIGURE 4 | Determination of the optimal number of reference genes. Pairwise variation (V n/n+1 ) were calculated between the normalization factors (NF n and NF n+1 ) in all tested samples. The V n/n+1 -values below 0.15 suggested that an additional reference gene is not necessary for normalization. genes (AgRPL24, AgRPS5, AgRPL15, AgRPL28, AgRPL27) have a relatively high abundance, whereas the AP2/ERF family genes (AgAP2-1, AgAP2-2, and AgANT belong to AP2 subfamily, AgRAP2 belongs to RAV subfamily, AgERF-4 belongs to ERF subfamily) were relatively low. Several genes which belonged to the same group showed a similar expression pattern, such as AgWAT1/AgAILP1, AgABCD19/AgARF1, and AgRPL24/AgRPL27. In addition, we also found that many genes showed differential expression at three development stages. Most AP2/ERF transcription factor genes and auxin-related genes showed low expression abundance at Stage1 and high abundance at Stages 2 and 3. It seemed that the expression of these genes  increased over the process of plant development. In contract, the expression levels of all the ribosomal protein genes at Stage 1 were significantly higher than those at Stages 2 and 3.

qPCR Analysis of Different Expression Genes in Three Development Stages
Among the selected genes, nine different expression genes were further selected to examine their expression levels in leaf blade and petiole at three development stages by qPCR (Figure 7). TUB-B was used as control gene to normalize the expression. The expression patterns of AgAP2-1 and AgAP2-2 were similar: the expression level was higher in petiole, especially at Stage 3. Overall, the expression of these two genes increased with the plant development. AgRAP2 had a little change in leave blade and petiole at three stages, yet the expression abundance kept a high level in Figure 5. The two auxin-related genes, AgARF1 and AgAILP1, also showed the similar expression pattern: the expression differences were not significant between Stages 1 and 2, but at Stage 3 the expression level increased about five times both in leaf and petiole. The other gene AgSNX1 seemed to have a stable expression. For the three ribosomal protein genes, the expression level of AgRPS5 in leaf blade had a slight increase, while in petiole showed no remarkable difference. Transcriptome data showed that this gene was detected with high abundance in three development stages (Figure 5). Interesting, the expression patterns of AgRPL15 and AgRPL28 in leaf blade were significantly decreased during the development of the leave blades, while the expressions in petiole were gradually increased.

DISCUSSION
qPCR has become an important tool for molecular biology research. Using a suitable reference gene can efficiently correct the errors of RNA quantity and reverse transcription efficiency, which can help to obtain the real differential expression of target gene (Udvardi et al., 2008). The most commonly used reference genes, which used as basic component of organelles skeleton (ACTIN, TUB-A, and TUB-B) or involved in biochemical metabolic processes of organisms (GAPDH, EF-1α, and UBQ), can stably express in tissues and organs Gutierrez et al., 2008). However, recent studies have found that these reference genes may show instability under various plant species or genotypes (Wang et al., 2014. For example, GAPDH showed the most stability in grapevine but ranked worst in wheat (Reid et al., 2006;Long et al., 2010); UBI and ACT showed good stability in wheat, yet performed unsatisfactory in tomato (Long et al., 2010;Mascia et al., 2010). Moreover, the expression stabilities of reference genes have been demonstrated to vary under different tissues and environmental stresses (Czechowski et al., 2005;Libault et al., 2008).
Since a large sequence data have been obtained in celery, the expression patterns and function analysis of many genes will be more convenient. As the vegetative organs, leaf blades and petioles are the product of specific development stage of celery. Previous study have revealed several key genes contribute to the complex network of celery development (Jia et al., 2015), further study involving the expression patterns and function analysis of these genes in various tissues and tissues developmental stages will be needed. To ensure accurate and reliable results, reference genes should be evaluated in target conditions. To date, there is no report on the selection of the most suitable reference genes in celery.
Three commonly used algorithms (geNorm, NormFinder, and BestKeeper) were used to evaluate and identify suitable reference genes (Vandesompele et al., 2002;Andersen, 2004;Pfaffl et al., 2004). In our study, geNorm ranked TUB-B, TUB-A, UBC as the best reference genes. NormFinder generated a similar ranking to the geNorm analysis. BestKeeper recommended UBC, EF-1α, and UBQ as the most stable reference genes under all samples. Taken all the results into consideration, TUB-B, TUB-A, and UBC can be used as reference for normalization in celery development. ACTIN is always considered as a suitable reference gene, but the results of the current study indicated that ACTIN is not the best suitable reference gene in celery. The tubulin gene family members, including TUB-A and TUB-B, also often serve as the suitable candidate reference genes. The expression stability of TUB-A is generally higher than the TUB-B (Brunner et al., 2004;Jian et al., 2008). But in this study, TUB-B appears to be more stable. Using a single reference gene for calibration and standardization is deemed to affect the accuracy of the result (Zhu et al., 2008). Schmid et al. (2003) suggested that two or more reference genes can contribute to the calibration of system deviation under a set of samples or experimental conditions. The geNorm programmer determined the optimal number of reference genes necessary for normalization under different samples in our study. With a threshold of 0.15, two reference genes were enough for normalization at Stages 1 and 2 in leaf blade and at Stage 3 in petiole, while more genes were needed for other tissue or conditions. In "Petiole" and "Total" groups, there was no suitable number of reference genes. The complexity of the samples may result in higher variability. However, using more reference genes can help to reach possible optimum results, but not a necessary criterion (Vandesompele et al., 2002).
A large number of genes are involved in the process of plant development. Some transcription factors, such as AP2/ERF, NAC, MADS-box, participated in cell differentiation, organ development, and construction of plant morphology (Rounsley et al., 1995;Köllmer et al., 2011;Le et al., 2011). In this study, genes encoding AP2/ERF family transcription factors expressed differently at different stages of celery development, suggested that these genes might involve in activity developmental regulation in celery. The expression levels of auxin-related genes and ribosomal protein genes were relatively high, especially the expressions of several ribosomal protein genes in early development stage. With the exuberant cell division activity, the young leaves require a lot of protein synthesis to meet the growing needs. The high expression of ribosomal protein genes in the initial stage of plant growth provides a prerequisite for translating other development-related genes. Some of the ribosomal protein genes associated with leaf development in Arabidopsis were also confirmed (Schippers and Mueller-Roeber, 2010). Overall, the coordinated regulation of a large number of development related genes has realized normal development of plants. Selection of the suitable reference genes provides a favorable basis for the further research on celery development.

AUTHOR CONTRIBUTIONS
AX and ML conceived and designed the experiments. ML, GW, QJ, CT, and FW performed the experiments. ML, QJ, CT, and AX analyzed the data. AX contributed reagents/materials/analysis tools. ML wrote the paper. ML and AX revised the paper.

ACKNOWLEDGMENTS
The research was supported by the National Natural Science Foundation of China (31272175); Jiangsu Natural Science Foundation (BK20130027); and Priority Academic Program Development of Jiangsu Higher Education Institutions.