Evaluation of housekeeping gene expression stability in carnation (Dianthus caryophyllus)

ABSTRACT The real-time quantitative reverse transcription PCR (RT-qPCR) is widely used for gene expression analysis, owing to its advantages of high specificity, sensitivity and repeatability. Suitable reference gene is an absolute prerequisite for accurate normalisation. In this study, three computational statistical methods before performing RT-qPCR, including geNorm, NormFinder and BestKeeper was used to integrated expression stability evaluations of 12 frequently-used reference genes in Dianthus caryophyllus across different experimental conditions. The results showed that the expression stability of candidate genes varies greatly in different sample pools. The expression of TIP41 (TIP41-like family protein) and UBQ10 (ubiquitin10) was relatively stable under different experimental conditions, while CYP (cytochrome P450) and TUA (a tubulin) could act as reliable internal controls in different tissues. The reliable reference genes under treatment of hormone, heavy metal, salt, heat, cold, flooding and drought were EF1α (elongation factor 1 alpha)/TIF5A (translation initiation factor 5A), UBQ10/18S (18S Ribosome RNA), UBQ10/CYP, TUB (β-tubulin gene)/TIP41 (TIP41-like family protein), TIF5A/EF1α, TUB/TIP41 and TIP41/PP2A (protein phosphatase 2A). Some common housekeeping genes, such as ACT (actin gene) and GAPDH (glyceraldehyde-3-phosphate dehydrogenase), exhibited unstable expression patterns. This was the first systematic study on the stability of internal reference genes of selection of reference genes in Carnation.


Introduction
Real-time quantitative reverse transcription PCR (RT-qPCR) is a nucleic acid quantification technology developed on the basis of traditional PCR technology (Ginzinger 2002;Bustin et al. 2005). It has strong specificity, sensitivity and high throughput. It has been widely used in molecular biology, medicine and other basic research fields, especially in gene function study, identification of strain and inspection and quarantine of animal and plant pests and diseases (Svingen et al. 2015;Oneill S et al. 2017;Umadevi et al. 2019). Because its sensitivity is too high, many factors could affect the correctness of expose the carnation to 4 or 42°C for 24 h. For other treatments, such as salt, drought, heavy metals, and hormone treatment, the seedlings are placed in NaCl, mannitol, CuSO 4 , CrCl 3 , ABA, NAA, 6-BA, GA solution for 24 h, and the solution concentration is 100 mM (Su et al. 2019). Both leaves and roots were sampled at the end of each treatment. To select reliable reference gene(s) specific to particular tissues, samples of the leaf, root, stem, sepal, peal, stamen and pistil were taken from non-treated plants. Each assay was represented by two biological replicates. All samples were frozen by liquid nitrogen, and stored in a −70°C refrigerator for use.

Primer design
12 candidate internal control genes, including TIF5A (translation initiation factor 5A), UBQ10, TIP41 (TIP41-like family protein), CYP, glyceraldehyde-3-phosphate dehydrogenase gene (GAPDH), elongation factor 1 alpha (EF1α), actin gene (ACT), 18S Ribosome RNA (18S), a tubulin and β-tubulin gene (TUA, TUB), S-adenosylmethionine decarboxylase gene (SamDC), protein phosphatase 2A (PP2A) were selected to screen the reliable reference genes with stable expression pattern, because the 12 genes were commonly used internal housekeeping genes, and their sequences were available in NCBI. The sequence was downloaded from NCBI and primers for RT-qPCR were selected using Primer3 software (v0.4.0) with annealing temperatures of 57-68°C, primer lengths of 18-20 bp, and approximately 50% GC content. The amplicon length is limited to 60-200 bp to ensure adequate amplification efficiency and minimise the impact of RNA integrity on the quantitation process. The primer sequences, amplicon sizes and melting points of all PCR products are shown in Table 1.

Total RNA extraction and cDNA synthesis
According to the manufacturer's protocol, the total RNA was isolated in samples using the Eastep® Super Total RNA Extraction Kit (TaKaRa, Japan), and the genomic DNA was removed by RNase-free DNase I (TaKaRa, Japan). The purity and concentration of RNA were determined by spectrophotometry and their integrity was observed by agarose electrophoresis. cDNA synthesis was carried out according to the PrimeScriptTM RT reagent Kit with gDNA Eraser [TaKaRa] operating instructions, in the presence of 500 ng RNA, 2.5 M oligo-(dT). A negative control was set up in the synthesis of cDNA, and all of them were added except reverse transcriptase.

Real-time PCR
RT-qPCR was performed using a CFX Connect™ Real-Time PCR Detection System and a TB Green® Premix Ex TaqTM II Kit. The RT-qPCR reaction was carried out in a 20 uL reaction, including 1.6 μl of cDNA, 0.8 μl of each primer (10 μM), and 10.0 μl of 2×TB Green Premix Ex Taq II (Tli RNaseH Plus). In addition, reactions without template or primer were usedas a negative control to exclude possible contamination. The thermal cycle program was set as the initial polymerase activation step, which was carried out for 30 s at 95°C, and then 40 cycles, including 15 s at 94°C for template denaturation, 15 s at 55-65°C for annealing, and 45 s at 72°C for extension and fluorescence detection. All samples were amplified in triplicate. The specificity of RT-qPCR was confirmed by agarose electrophoresis and dissociation curve analysis. The cycle threshold for each reaction (Cq, the first cycle of the signal over the background) was automatically determined, by default parameters of the CFX Connect™ Real-Time PCR Detection System device.

Data analysis
Based on the original fluorescence data obtained from Applied Biosystems 7300 Real-Time PCR System equipment, the RT-qPCR efficiency and all correlation coefficients of each gene were calculated by Lin-RegPCR program. The Cq value (Livak and Schmittgen 2001) used in the analysis represents the average of 6 values (3 biological repetition and 3 technical repetition), and were transformed to relative quantities using the comparative Cq method, according to the instructions of geNorm v3.5, NormFinder v0.953, BestKeeeper v1.0 (www.gene-quantifification.info) and RefFinder (https:// heartcure.com.au/) (Vandesompele et al. 2002).

Expression profiles of reference genes
The Cq value was negatively correlated with the initial template quantity of RT-qPCR (Gutierrez et al. 2008). The higher the Cq value, the less the gene content. As shown in Figure 1, the gene with the highest expression level was 18S (the average Cq value of the gene was the lowest), while the lowest expressed gene was ACT. Although the amount of RNA that initially added in RT-qPCR reaction was similar, the Cq value of each gene varied considerably across different experimental conditions, indicating that there were no genes with absolutely stable expression pattern ( Figure 1). Among all the genes, the largest change in Cq value was ACT, which seemed to suggest that the gene expression was the most unstable. The TIP41 has the smallest change in Cq value, indicating that the expression of this gene may be relatively stable. However, the stability of gene expression ultimately required further evaluation, and reliable reference genes were evaluated according to the specific situation.
geNorm analysis geNorm software is the most commonly used tool for evaluating the stability of internal reference genes. It is mainly based on the Microsoft Excel platform and is a VBA macro program. The algorithm is mainly evaluated by introducing the stability value M of the internal reference gene expression, and the magnitude of the M value is negatively correlated with the stability of the gene, that is the smaller the M value, the more stable the gene (Hellemans et al. 2007). Genes with a threshold M value lower than 1.5 are considered to be stable (Song et al. 2020). As shown in Table 2, when looking at the expression data set as a whole, EF1α and TIP41 showed the lowest M value therefore, the two genes could be the reliable internal control genes in this sample pool, while across different tissues, TUA and CYP exhibited the most stable expression pattern. In the treatment of hormone, heavy metal, salt, heat, cold, flooding and drought, the suitable reference genes could be EF1α/UBQ10, TIP41/18S, CYP/UBQ10, TUB/TIP41, TIP41/ EF1α, TIP41/TUB, TIP41/TIF5A. Generally speaking, when performing analysis in RT-qPCR, increasing the number of internal reference genes may improve the accuracy of the experiment. In order to determine the appropriate number of internal reference genes under specific experimental conditions, geNorm software was used to analyse the variation of pairs of internal reference genes in this experiment in order to calculate the appropriate number of genes. Researchers usually used 0.15 as the critical value for selecting the best internal parameter genes (Vandesompele et al. 2002). For example, in all samples, the value of V 2/3 was higher than 0.15, indicating that the first two reference genes were not enough to accurately quantify the results; while V 3/4 was lower than 0.15, indicating that the first three reference genes could be used as internal reference genes to obtain more accurate results. In addition, if you added the fourth-ranked gene as an internal reference gene, it would not only increased the cost of the test but would also affected the accuracy of the results ( Figure 2). Certainly, 0.15 was not an absolute critical value. The selection of critical value also needs to consider the sample size, variation range and so on (Vandesompele et al. 2002;Gu et al. 2011). In addition, the number of internal reference genes also needed to consider the test cost and the required test accuracy. As shown in Figure 2, except for all sample groups, the reliable results of gene expression analysis could be obtained only by introducing the first two internal reference genes in other words sample pools.

NormFinder analysis
The software is an alternative tool that is often used to screen stable internal reference genes. Similar to geNorm, NormFinder is also a VBA macro program written on the Microsoft Excel platform. By comparing the expression variation of each gene, the expression stability coefficient (expression values) of each gene under certain experimental conditions was calculated, and this coefficient was also negatively correlated with the stability of gene expression (Cruz et al. 2009;Hong et al. 2010;Gopalam et al. 2017). In addition, the software can group samples over a wide range of samples in order to obtain more accurate results. In this trial, when analysed all the samples, the 33 samples were divided into two groups: untreated and treated. According to the analysis of NormFinder software (Hu et al. 2014) (Table 3), TUB and EF1α exhibited the most stable expression pattern, whether grouped or not, although the order of other genes is slightly different. GAPDH and TUB expressed stable across different tissues, followed by EF1α and UBQ10. When the plant was treated with hormone, heavy metal, salt, heat, cold, flooding and drought, UBQ10/EF1α, UBQ10/ TIP41, PP2A/CYP, ACT/SamDC, EF1α/GAPDH, TUB/UBQ10 and TUA/PP2A exhibited the most important invention stable expression pattern.

BestKeeper analysis
BestKeeper is a program written for the analysis of internal reference genes and target gene expression. The judgment principle is that the smaller the standard deviation and the coefficient of variation, the better the stability of the internal reference gene, otherwise, the worse the stability; when SD > 1, the expression of the internal reference gene is unstable, and more stable genes have a lower value of CV ± SD. As shown in Table 4, when all 33 samples were taken together, TIP41 and EF1α showed the highest value for expression stability. In different tissues, the most stably expressed genes were TIP41 and SamDC, under the treatment of hormone, heavy metal, salt, heat, cold, flooding and drought, the stable reference genes were TIF5A/TIP41, ACT/ GAPDH, CYP/UBQ10, CYP/18S, TIF5A/TIP41, GAPDH/ACT and ACT/TIP41.
consistent. RefFinder can integrate the currently available major computational programs (geNorm, Normfinder, BestKeeper, and the comparative Delta-Ct method) to compare and rank the tested candidate reference genes Xie et al. 2012). Based on the rankings from each program, It assigns an appropriate weight to an individual gene and calculated the geometric mean of their weights for the overall final ranking. As shown in Table 5, when all 33 samples were taken together, TIP41 and UBQ10 were the reliable reference genes, while CYP and TUA exhibited the stable expression pattern across different tissue. In the treatment of hormone, heavy metal, salt, heat, cold, flooding and drought, the suitable reference gene was EF1α/TIF5A, UBQ10/18S, UBQ10/CYP, TUB/TIP41, TIF5A/EF1α, TUB/TIP41 and TIP41/PP2A. According to our results (Table 5), TIP41 and UBQ10 are the most stable genes for all samples, can be used as reference to normalise target gene expression.

Effects of stable or poor reference genes on gene expression analysis
In order to prove the influence of reference gene stability on the results of gene expression analysis, the two most stable and the most unstable genes were selected as reference genes, and the relative expression levels of the remaining eight genes were calculated. Compared the determination coefficient R 2 and correlation coefficient |r| of the two best reference genes and the most unstable genes were.  1  TIP41  TIP41  TIF5A  ACT  CYP  CYP  TIF5A  GAPDH  ACT  2  EF1α  SamDC  TIP41  GAPDH  UBQ10 18S  TIP41  ACT  TIP41  3  PP2A  EF1α  EF1α  PP2A  PP2A  PP2A  EF1α  UBQ10  EF1α  4  TUB  TUA  PP2A  CYP  EF1α  UBQ10 GAPDH TIP41  SamDC  5  GAPDH  GAPDH  UBQ10  TUA  TUB  SamDC UBQ10 EF1α  CYP  6  ACT  CYP  TUB  TIF5A  TIP41  TIP41  TUB  TUA  GAPDH  7  SamDC  TUB  SamDC  TIP41  TIF5A  TUB  SamDC TIF5A  PP2A  8  TIF5A  ACT  ACT  SamDC  TUA  TIF5A  ACT  CYP  TUA  9  TUA  UBQ10  CYP  UBQ10  ACT  EF1α  PP2A  SamDC  UBQ10  10  UBQ10  PP2A  GAPDH  TUB  SamDC TUA  TUA  PP2A  TIF5A  11  CYP  TIF5A  TUA  EF1α  GAPDH ACT  CYP  TUB  TUB  12 18S As shown in Table 6, the results of there was no significant difference in the reference genes between the two stable genes in each group, but the results of the two most unstable genes were quite different. This reflected two problems. First, the results of the stability evaluation of reference genes in this study were reliable. Second, the selection of unstable reference genes would have a great impact on the results of gene expression, and even the opposite result would be obtained.

Discussion
Because of its high specificity, sensitivity and high throughput, RT-qPCR has been widely used to analyse gene expression, animal and plant virus detection, etc. (Tang et al. 2007;Yee et al. 2018;Silva et al. 2019). However, if the unstable internal reference gene was selected in the experiment, it will seriously affect the accuracy of the experimental results, and even get the completely opposite results. The absolutely stable internal reference gene does not exist, even the same gene, in different tissues, at different developmental stages, and under different treatments, its expression stability is different. For example, ACT11 has been proved to be a stably expressed gene in Soybean (Hu et al. 2009). However, the expression of ACT11 in rice is unstable (Jain et al. 2006), while in ryegrass, its stability is the lowest (Dombrowski and Martin 2009), similar to the results of this experiment.
There was much software to evaluate the stability of gene expression (Pfaffl et al. 2004), and the softwares were based on different principles. Which software was authoritative has not been determined. Therefore, in this experiment, the three most commonly used softwares, geNorm, NormFinder and BestKeeper were selected evaluate the gene expression stability. It was shown that under most experimental conditions, the output of these three softwares was generally consistent. However, in different tissues and heavy metal stress, gene stability was quite different. The difference in results may be due to the difference in algorithms between the two tools, or to the existence of co-regulation between genes, because the results of geNorm would be interfered by gene co- regulation, while the results of NormFinder and BestKeeper will not (Willems et al. 2006). In order to test whether there was gene co-regulation in these genes, one of the 12 candidate genes was removed one by one, and then the data was introduced into the geNorm software for calculation, to see whether there was a significant change in the rank of other genes. The results showed that when removing any gene, there was no significant change in the remaining gene sequence. Therefore, the results of the three tools were different only because of their different algorithms (Tang et al. 2007;Expositorodriguez et al. 2008;Hong et al. 2010). As shown in Table 5, the order of stability of candidate reference gene expression was not consistent under different experimental conditions in Carnation. For example, TIP41 is the most stable internal control gene when looking at the expression data set as a whole; while in salt stress treatment, TIP41 was not at the top of the list. Before selecting the internal control genes, we should consider the development stages, tissues and treatment methods of plants, because the stability of gene expression was related to the experimental conditions.
To the best of our knowledge, this study is the first to analyse qPCR reference genes of carnation under different conditions. Based on the analysis results of geNorm, NormFinder, BestKeeper and NefFinder, when looking at the expression data set as a whole, the most important invention stable reference genes were TIP41 and UBQ10, while CYP and TUA ranked the top. When carnations were treated with hormone, heavy metal, salt, heat, cold, flooding and drought, the appropriate reference genes were EF1α/TIF5A, UBQ10/18S, UBQ10/CYP, TUB/TIP41, TIF5A/EF1α, TUB/TIP41 and TIP41/PP2A. This study is of great significance for accurate quantitative and expression analysis of carnation genes under different conditions. In the future research, it will play an important role in carnation breeding, such as the study of Tolerance Varieties and the analysis of chemical metabolism pathway.

Disclosure statement
No potential conflict of interest was reported by the author(s).