Selection of appropriate reference genes for RT-qPCR analysis under abiotic stress and hormone treatment in celery

Celery is one of the most important vegetable crop and its yield and quality is influenced by many environmental factors. Researches on gene expression not only help to unravel the molecular regulatory mechanism but also identify the key genes in the biological response. RT-qPCR is a commonly used technology to quantify the gene expression. Selecting an appropriate reference gene is an effective approach to improve the accuracy of RT-qPCR assay. To our knowledge, the evaluation of reference genes under different treatments in celery has not been reported yet. In this study, the expression stabilities of eight candidate reference genes (ACTIN, eIF-4α, GAPDH, TBP, TUB-A, UBC, TUB-B, and EF-1α) under abiotic stresses (heat, cold, drought, and salt) and hormone treatments (SA, MeJA, GA, and ABA) were detected. The expression stabilities of candidate genes were compared and ranked by geNorm, NormFinder, BestKeeper, ΔCt, and RefFinder programs. The results calculated by different programs were not completely consistent. Considering the comprehensive analysis results, ACTIN was the most stable reference gene and TUB-B showed the worst expression stabilities under the selected abiotic stress and hormone treatments in celery. The reliability of reference genes was further confirmed by the normalization of CAT1 gene under drought stress. This study presented evidences and basis to select the appropriate reference genes under different treatments in celery.


INTRODUCTION
Celery (Apium graveolens L.), one plant of Apiaceae, is an important vegetable and its leaves are the mainly edible organs (Li et al., 2018). Nowadays, celery is commonly consumed for its abundant nutritional values (apigenin, vitamin C, and cellulose etc.) and low calorie contents (Dianat et al., 2015). The yield and quality of celery are influenced by many environmental factors (temperature, moisture, soil salinity, and hormone) (Golldack, Luking & Yang, 2011). During plant development, many environmental stresses disturb the physiological processes and affect the growth and development (Cattivelli et al., 2008;Kosova, Prasil & Vitamvas , 2008;Wang et al., 2012;Wang et al., 2014). Phytohormones are known to be plant growth regulators and play vital roles during plant development, such as gibberellic acid (GA), salicylic acid (SA), abscisic acid (ABA), and methyl jasmonate (MeJA) (Moons et al., 1997;Wang et al., 2015;Liu et al., 2016;Galimba et al., 2019). Under abiotic stress and hormone treatments, plants generate many responsive mechanisms to relieve environmental damages (Kosova, Vitamvas & Prasil , 2014). The molecular mechanisms including physical, physiological, and biochemical responses were associated with the expressions of certain genes (Wang, Vinocur & Altman, 2003). Researches on expressions of abiotic stress-related genes provided strategies to improve the stress resistance in molecular breeding (Brikis et al., 2018;Wang et al., 2019).
The gene expression analysis was commonly applied to understand the molecular regulatory mechanisms and identify the key genes in the current molecular biology (Bustin et al., 2005;Silva et al., 2019). Quantitative real time polymerase chain reaction (RT-qPCR) has become a recognized technology for quantifying the gene expression due to its advantages of high-throughput, high-sensitivity, high-veracity, and low-cost (Gachon, Mingam & Charrier , 2004;Bustin et al., 2005;Derveaux, Vandesompele & Hellemans , 2010;Miao et al., 2019). Nevertheless, many factors including enzymatic efficiency, RNA purity, and cDNA quality may affect the accuracy and credibility of RT-qPCR results (Vandesompele et al., 2002). Several strategies were applied to ensure the accuracy of RT-qPCR. Selection of one or more suitable internal control genes, also known as housekeeping genes, has become a frequently method to normalize the gene expression (Gutierrez et al., 2008).
The house-keeping genes have been identified in many species, including Arabidopsis (Czechowski et al., 2005), rice (Jain et al., 2006), and soybean (Libault et al., 2008). Housekeeping genes, e.g., glyceraldehyde-3-phosphate dehydrogenase (GAPDH ), ubiquitin C (UBC), actin (ACTIN ), eukaryotic translation initiation factor 4α (eIF -4α), elongation factor-1α (EF -1α), TATA-box binding protein (TBP), and tubulin (TUB) were widely used as reference gene to standardize the expressions of target genes (Dheda et al., 2004;Galli et al., 2013;Li et al., 2016b). However, some studies also indicated that the expressions of certain house-keeping genes under different tissues or treatments were fluctuant (Barber et al., 2005;Borowski et al., 2014;Tian et al., 2015). The unstable reference gene would significantly influence the accuracy and reliability of target gene expression quantification. The suitable reference gene for RT-qPCR analysis of celery under different tissues and developmental stages has been identified in previous study (Li et al., 2016b). To our knowledge, the appropriate reference genes for RT-qPCR analysis of celery under abiotic stress and hormone treatment have not been reported yet. Considering the roles of gene expression in the molecular biology of celery, the comparison and selection of reference genes under various experimental treatments is necessary.

Plant materials and experimental treatments
The seeds of celery cultivar 'Jinnan Shiqin' were germinated in a petri dish at room temperature. Celery seedlings were transferred into the plastic pots with 1:1 mixture of soil and vermiculite. Seedlings were grown in an artificial climatic chamber with the condition as previously described (Feng et al., 2018c). After 8 weeks of growth, the vigorous seedlings with consistent growth were selected for experimental treatments. As for heat and cold stresses, seedlings were placed in the light incubators with temperatures of 38 • C and 4 • C, respectively. For drought and salt stresses, seedlings were irrigated with 0.5 L of PEG 6000 (20%) and NaCl (0.2 M) solution, respectively (Tian et al., 2015). As for hormones treatments, celery leaves were sprayed with 0.5 L of GA (1.4 mM), SA (1.4 mM), ABA (0.1 mM), and MeJA (0.8 mM), respectively (Chan, 2012;Zhu et al., 2013;Li et al., 2016a). All of the treatments were performed with three biological replicates. Leaf blades were collected from untreated and treated celery plants after 2 h of treatments.

Preparation of RNA and cDNA
Total RNA was extracted from the celery samples by using Total RNA Kit (Tiangen, Beijing, China) based on manufacturer's protocol. The concentration and quality of total RNA was measured by using a One-Drop TM spectrophotometer. The qualified RNA (1 µg) were used to synthesize cDNA by using the Prime-Script RT reagent kit (TaKaRa, Dalian, China) with a 20 µL system.

RT-qPCR analysis
Eight candidate celery genes, ACTIN, eIF -4α, GAPDH, TBP, TUB-A, UBC, TUB-B, and EF -1α, were used to screen the appropriate reference genes under different abiotic stresses and hormone treatments based on the celery transcriptome and genome data (Li et al., 2014;Li et al., 2016b;Feng et al., 2018a). The primer sequences of ACTIN, GAPDH, TBP, TUB-A, UBC were consistent with previous study (Li et al., 2016b). The gene sequences of eIF -4α, TUB-B, and EF -1α cloned from 'Jinnan Shiqin' were different from the previous study, which were listed in Table S1. The RT-qPCR primer sequences of TUB-B, eIF -4α, and EF -1α genes were re-designed by using Primer Premier 6.0 software. The gene information and primer sequences were listed in Table 1. The specificity and accuracy of primers were determined by the PCR assay and single peak in the melting curve of RT-qPCR assay.

Data analysis
As for the primers of eIF -4 α, TUB-B, and EF -1α genes, the standard curves were established based on the Cq values of different dilution gradient with their corresponding logarithm of dilution multiples. The slope of standard curve was used to calculate the amplification efficiency (E) of primers, according to the formula: E% = (10 [−1/slope] -1) ×100% (Radonic et al., 2004). The amplification efficiency of other primers was reported in previous study (Li et al., 2016b). The stabilities of candidate reference genes under different treatments were ranked based on the analysis results of geNorm (Vandesompele et al., 2002), NormFinder (Andersen, Jensen & Orntoft , 2004), BestKeeper (Pfaffl et al., 2004), Ct (Silver et al., 2006), andRefFinder (Xie et al., 2012). Nine Cq values were obtained from each sample, including three biological and three technical replicates. Before geNorm and NormFinder analysis, the raw Cq values should be calculated with the 2 − Ct formula ( Ct indicated the Cq value of the sample minus the minimum Cq value). The geNorm program calculated the M value of each gene, and two genes with the lowest M value were the most stable reference genes. The pairwise variation (Vn/n+1) in geNorm program indicated the optimal number of reference gene for the normalization of RT-qPCR. If the Vn/n+1 value <0.15, the optimal number of optimal internal reference genes is n; adversely, the optimal number of reference genes is n + 1. NormFinder ranked the candidate genes according to the stability value calculated by the expression variations of intragroup and the intergroup of each gene. In BestKeeper, the standard deviation (SD) and coefficient of variation (CV) were calculated based on the Cq values. BestKeeper ranked the stability of candidate genes according to the SD and CV values. The lowest SD and CV values indicated the most stable reference gene. The expression stabilities of 'pairs of genes' were compared by Ct program. Based on the analysis results of various programs, the stability of candidate reference genes were comprehensively evaluated and ranked by using RefFinder.

Validation of reference genes
Catalase was involved in plant regulatory mechanism under abiotic stress (Hu et al., 2010). Based on the protein sequence of AtCAT2 (GenBank accession number NP_195235.1) and our transcriptome and genome data, the celery CAT1 gene was identified and cloned from 'Jinnan Shiqin'. The sequence of CAT1 has been submitted to GenBank (accession number: MN365877). The RT-qPCR primer of CAT1 was designed by using Premier 6.0 software according to the gene sequence (forward: 5 -TTCACCTTCCTCTTGGATGACATTGG-3 and reverse: 5 -GCTCCTCCGATCTTGATGGCTTC-3 ). The RT-qPCR assay of CAT1 gene under drought treatment was conducted to validate the candidate reference genes. The relative expression level of CAT1 gene was normalized using various reference genes according to the 2 − Ct method (Schmittgen & Livak, 2008).

Selection of candidate reference genes
Eight genes, ACTIN, eIF -4α, GAPDH, TBP, TUB-A, UBC, TUB-B, and EF -1α were selected as candidate reference genes from celery. The specificity and efficiency of primers were confirmed by PCR amplification assay and melting curve of RT-qPCR. The single band corresponding to various reference genes was detected in the electrophoretogram with 1.5% agarose gel, respectively (Fig. S2). The melting curves showed that candidate reference genes had a single peak in RT-qPCR reaction (Fig. S3). The amplification efficiency (E) and correlation coefficient (R 2 ) of the candidate genes meets the standard of RT-qPCR assay (90% <E % <110%; R 2 > 0.99; Table 1) (Li et al., 2016b).

RT-qPCR results of candidate reference genes
The gene expression levels were indicated with the Cq values in RT-qPCR assay. The raw Cq values of candidate reference genes were listed in Table S2 and their statistics were demonstrated in Fig. 1 (Fig. 2).

Stability analysis of candidate reference genes
Five programs, geNorm, NormFinder, BestKeeper, Ct, and RefFinder were used to determine the expression stability of the celery candidate reference genes. To evaluate the gene stabilities under different treatments, celery plants were subjected with 8 treatments, abiotic stresses (heat, cold, drought, and salt) and hormone treatments (SA, MeJA, GA, and ABA). In addition to the stability analysis of single-treatment, these 8 treatments were also divided into three groups, namely abiotic stress (heat, cold, drought, and salt), hormone stimuli (GA, SA, ABA, and MeJA), and total (all treatments), for stability analysis.
geNorm analysis geNorm program calculated the M values of candidate reference genes and ranked their stabilities based on the M values. The lowest M value represents the most stable expression stability. As shown in Table 2, the M value of reference genes under all treatments were less than the default limit (1.5), which indicated that their expression stability were satisfactory. EF -1α and TBP genes were the most stable reference genes with the lowest M value under cold, drought, SA, and GA treatments. Meanwhile, TBP was also showed the highest expression stability under heat and MeJA treatments. Under salt treatment, ACTIN and TUB-A showed the lowest M value and they were the most stable reference genes. As for the abiotic stress and total groups, the M values of ACTIN and EF -1α were the lowest, which indicated that they were the most stable reference genes ( Table 3). As for the hormone stimuli group, TBP and GAPDH showed higher high expression stability than others. In the stability analysis of all three groups, TUB-B was the worst stable reference gene with the lowest M value.
In RT-qPCR analysis, multiple reference genes can be selected to quantify the expression of the target gene with more accuracy. The pairwise variations (Vn/n+1) calculated by geNorm program were used to determine the optimal number of reference genes in the normalization. As shown in Fig. 3, the V2/3 values of eight single-treatment and the three treatment groups were less than 0.15. Therefore, two suitable reference genes were adequate for gene expression normalization under the above treatments.

NormFinder analysis
NormFinder program ranked the reference genes based on the calculated stability values in different experimental designs. As shown in Tables 2 and 3, ACTIN was the most stable reference gene under heat and salt treatments, and EF -1α was the most stable gene under cold, drought and all tested hormone treatments. When performing stability analysis among multiple treatments, ACTIN was the most stable reference gene with the lowest stability value in both the abiotic stress group and the total group. As for the hormone stimuli group, the expression stability of EF -1α gene was the highest, followed by ACTIN gene. Similar to the geNorm analysis, eIF -4α, UBC, and TUB-B were the worst stable reference genes in the NormFinder analysis.

BestKeeper analysis
The BestKeeper program ranked the reference genes based on the SD and CV of Cq values in the RT-qPCR assay. Low SD and CV values represent the high expression stability. UBC was the most stable reference gene under drought, SA, and ABA treatments. EF -1α and ACTIN were the most stable reference genes under salt and MeJA treatments, respectively. As for the abiotic stress group and the total group, TUB-A showed the highest expression  stability, and eIF -4α showed the lowest expression stability. UBC was the most stable reference gene and GAPDH was the least stable reference gene in hormone stimuli group.

Ct analysis
The expression stability of the eight candidate gene was calculated and ranked based on the Ct method. In single treatment, ACTIN was the most stable expressed reference gene under heat, salt, SA, and MEJA treatments, respectively. As for drought, GA, and ABA treatments, EF -1α gene was the best reference gene for gene normalization. The expression stability under different groups was also investigated. ACTIN gene showed the most expression stability under abiotic stress treatments. In hormone stimuli and total groups, the expressions of ACTIN and EF -1α genes showed the highest stabilities. In addition, Ct analysis results indicated that TUB-B was the worst reference gene under most treatments.

RefFinder analysis
Considering the results of all statistic methods, the stability of those celery candidate reference genes was comprehensively evaluated by RefFinder. As shown in Table 2, ACTIN gene was the most stable reference gene under heat and salt treatments. EF -1α showed the  most stable expression under cold, drought, MEJA, GA, and ABA treatments. As for the stability analysis in different groups, ACTIN gene was the most stable reference gene and TUB-B was the worst stable reference gene in the abiotic stress, hormone stimuli, and total groups ( Table 3).

Validation of reference genes
The CAT gene encoded the catalase, which is involved in the regulation of stress defense in plants (Willekens et al., 1997). The expression of the CAT gene could be induced by many abiotic stresses, including chilling, drought, and salt (Fadzillah et al., 1996;Kim et al., 2007). The celery CAT1 gene was cloned from the cDNA of 'Jinnan Shiqin' and sequenced. In this study, the relative expression level of CAT1 gene under drought stress was detected to validate the reference genes. As shown in Fig. 4, the expression levels of CAT1 gene normalized by various reference genes were different. The expression levels of CAT1 gene under drought stress were increased using ACTIN, GAPDH, TBP, TUB-A, UBC, TUB-B, and EF -1α as reference gene, respectively. When used the unstable reference gene eIF -4α for normalization, the expression of the CAT1 gene was decreased after 24 h of drought treatment.

DISCUSSION
Gene expression plays important roles in plant development and environmental stimuli defense. Research on gene expression contributed to unravel the complex regulatory mechanisms in life cycle of plant (Feng et al., 2018b;Silva et al., 2019). Nowadays, RT-qPCR is a general technique to determine the expression level of target gene (Gachon, Mingam & Charrier , 2004). However, the accuracy of RT-qPCR results was affected by many factors. Using proper reference gene was an effective approach to improve the accuracy of gene normalization during RT-qPCR assay (Nicot et al., 2005;Gutierrez et al., 2008). A previous study has investigated the suitable reference genes among various tissues and development stages in celery (Li et al., 2016b). In the process of growth and development, celery also encounters many environmental stimuli, including biotic and abiotic stresses. Selection of suitable reference genes is crucial to normalize the gene expression under different conditions in celery. The current study evaluated the expression stabilities of various candidate reference genes under abiotic stress and hormone stimuli in celery.
In this work, the expression profiles of 8 candidate reference genes of celery (ACTIN, eIF -4 α, GAPDH, TBP, TUB-A, UBC, TUB-B, and EF -1 α) under different abiotic stress and hormone stimuli were determined and compared. These candidate genes were commonly different treatments in celery. As the most stable reference gene in celery, ACTIN was also investigated to be the suitable reference gene in carrot and soybean under abiotic stress treatments (Jian et al., 2008;Tian et al., 2015). The ACTIN was the most stable gene at different development stages in carrot (Wang et al., 2016). It should be noted that the TUB-B investigated to be the most stable gene among tissues and developmental stages of celery (Li et al., 2016b), whereas our study indicated that TUB-B was the least stable gene under different treatments. This indicated that the stability of the same house-keeping gene was various under different conditions. Celery generated physiological regulation through the expressions of specific genes during development and environmental stress. Plant accumulated amounts of hydrogen peroxide (H 2 O 2 ) under stress conditions (Willekens et al., 1997). The CAT gene encoded the catalase, which is involved in the regulation of H 2 O 2 level in plants. Previous study indicated that the expression of CAT gene was up-regulated under drought stress (Nie et al., 2015). To validate the reliability of reference genes, the relative expression level of CAT1 gene was normalized by using different reference genes. Except when using the unstable reference gene eIF -4 α, the expressions of CAT1 normalized by other reference genes were increased under drought treatment.

CONCLUSIONS
This work aims to select the appropriate reference gene under abiotic stress and hormone stimuli in celery. The stability of eight candidate reference genes under different treatments was evaluated and ranked by geNorm, NormFinder, BestKeeper, Ct, and RefFinder programs. The analysis results indicated that ACTIN was the most recommended reference gene under abiotic stress and hormone treatments in celery, whereas the TUB-B was the worst stable gene. The reliability of celery reference gene was verified by expression normalization of CAT1 gene under drought stress. In conclusion, the results in this study provided reference and basis for the selection of suitable reference genes under abiotic stress and hormone treatment in celery.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This study was supported by the Jiangsu Agricultural Science and Technology Innovation Fund [CX (2018)2007], New Century Excellent Talents in University (NCET-11-0670); the National Natural Science Foundation of China (31272175); and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). There was no additional external funding received for this study. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: Jiangsu Agricultural Science and Technology Innovation Fund: CX(2018)2007.