Screening and Evaluation of Stable Reference Genes for Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) Analysis in Chinese Fir Roots under Water, Phosphorus, and Nitrogen Stresses

: Chinese ﬁr ( Cunninghamia lanceolata ) is an economical important timber species widely planted in southeastern Asia. Decline in yield and productivity during successive rotation is believed to be linked with abiotic stress, such as drought stress and nitrogen (N) and phosphorus (P) starvation. Molecular breeding could be an option to develop tolerant genotypes. For gene expression studies using quantitative real-time reverse transcription-polymerase chain reaction (qRT-PCR), stable reference genes are needed for normalization of gene expression under di ﬀ erent experimental conditions. However, there is no internal reference genes identiﬁed for Chinese ﬁr under abiotic stresses. Thus, nine internal reference genes based on transcriptome data were selected and analyzed in the root of Chinese ﬁr under drought stress and N and P starvation. Data were analyzed using geNorm, NormFinder, and BestKeeper, to screen and identify the best reference genes. The results showed that the UBQ and GAPDH genes were the two most stable genes under drought stress and the Actin1 and GAPDH were the two most stable genes under P starvation. Further, it was discovered that the Actin1 and UBC were the two most stable genes under N starvation among nine candidate reference genes. The gene expression of drought stress induced expression protein 14-3-3-4, the P transporter gene ClPht1;3 , and the nitrate transporter gene NRT1.1 were used to verify the stability of the selected reference genes under drought stress and P and N starvation, respectively, and the results revealed that the screened reference genes were su ﬃ cient to normalize expression of the target genes. In conclusion, the results demonstrate that the stability of reference genes was closely related to the external conditions and reference genes applied to the roots of Chinese ﬁr under di ﬀ erent abiotic stress treatments were di ﬀ erent. Our data will facilitate further studies on stress ecology and gene function analysis in Chinese ﬁr. and 1.0 mmol / L NH 4 NO 3 , representing N-starved, low N supply and normal N supply, and the roots were harvested after 0, 1, 3, 7, 15, and 25 days post-treatment. The gradient of drought stress was established by adding 0%, 10%, and 15% PEG6000, representing non-stress, moderately stress and severely stress condition, to one-third content of modiﬁed Hoagland solution, and the roots were harvest after 12, 24, and 48 h post-treatment. The ﬁne roots of all samples were harvested, washed and surface dried and then immediately frozen in liquid nitrogen and stored at − 80 ◦ C until analysis. All experimental treatments had three biological replicates. we recommend UBQ and the Actin1 gene for drought, P, and N stresses as superior internal controls for normalization of qRT-PCR, respectively. Additionally, our results showed that di ﬀ erent suitable reference genes or standardized reference gene combinations could be screened according to di ﬀ erent external conditions.

Gene expression studies on forest trees have become important to unravel the adaptation mechanisms to environmental conditions at molecular level. Drought stress and N and P starvation are key problems of forest tree growth in general and Chinese fir cultivation in particular, and the molecular research to address these problems depends on a prerequisite of screening and evaluation of the stability of reference gene for normalization of qRT-PCR gene expression analysis. Based on the transcriptome data of Chinese fir roots under drought stress and P and N starvation, we have chosen nine candidate reference genes including GAPDH, UBQ, 18S rRNA, 28S rRNA, EF1a, UBC, Actin1, elF-3, and CYP and examined the expression stability. Data were analyzed using geNorm, NormFinder, and BestKeeper for comprehensive and accurate determination of the most stability genes out of the nine reference genes under drought stress and N and P starvation in the root of Chinese fir seedlings [40]. To validate our results, the most stable and unstable genes were used to standardize the expression levels of 14-3-3-4, ClPht1;3 and NRT 1.1 genes that are related to drought stress and N and P starvation on different experimental conditions. Our results provide a basis for normalization of gene expression in other coniferous species, such as Douglas fir.

Plant Material and Experimental Setup
One year old seedlings of Chinese fir clone Yang 061 were bought from Yangkou State-owned forest farm in Fujian Province, China and cultivated in the greenhouse in Fujian Agriculture and Forestry University. The seedlings were cultivated in polyethylene pots with a volume of 1 L filled with cleaned sand and supplied with a modified (one-third content) Hoagland solution for two weeks. The solution contained 1 mmol/L KH 2 PO 4 , 5 mmol/L KNO 3 O 4 , and the pH of the nutrient solution was adjusted to 5.6. After two weeks of growth, the seedlings were transferred to another pot filled with thoroughly cleaned sand to make sure that traces of N and P were negligible. The concentration gradient of nutrient solution in P starvation treatment was as follows: 0, 0.5, and 1.0 mmol/L KH 2 PO 4 , representing P-starved, low P supply and normal P supply, respectively, and the roots were harvested after 0, 1, 3, 7, and 15 days post-treatment. The N treatment included different nitrogen sources and concentrations of nitrogen. As N sources, 1 mmol/L (NH 4 ) 2 SO 4 , 2 mmol/L KNO 3 , and 1 mmol/L NH 4 NO 3 were used to represent NH 4 + , NO 3 − , and mixed nitrogen sources, respectively. The concentration gradient of nutrient solution in N starvation test was as follows: 0, 0.5, and 1.0 mmol/L NH 4 NO 3 , representing N-starved, low N supply and normal N supply, and the roots were harvested after 0, 1, 3, 7, 15, and 25 days post-treatment. The gradient of drought stress was established by adding 0%, 10%, and 15% PEG6000, representing non-stress, moderately stress and severely stress condition, to one-third content of modified Hoagland solution, and the roots were harvest after 12, 24, and 48 h post-treatment. The fine roots of all samples were harvested, washed and surface dried and then immediately frozen in liquid nitrogen and stored at −80 • C until analysis. All experimental treatments had three biological replicates.

RNA Isolation and cDNA Synthesis
Total RNAs were extracted from Chinese fir roots using the RNA purification reagent (TIANGEN, Beijing, China). The frozen specimens were ground in liquid nitrogen to a fine powder with a pestle and mortar, and the powder was completely dissolved in the plant lysate, following the manufacturer's protocol. The integrity of obtained RNA samples was examined by 1% agarose gel electrophoresis. The concentration and purity of total RNA was examined using an ultramicro spectrophotometer (Eppendorf, Hamburg, Germany). The A 260 /A 280 ratio of RNA between 1.8 and 2.0 was considered to be the required quality for further experiments. Complementary DNA (cDNA) was synthesized using the GoScript TM Reverse Transcription System Kit (Promega, Madison, WI, USA), following the manufacturer's protocol. The reverse transcription system was based on 2 µg of total RNA in a 20 µL reaction volume. The resulting cDNAs were diluted to 100 ng/µL with nuclease-free water and stored at −20 • C.

Selection of Candidate Reference Genes
Based on the transcriptome data of Chinese fir, the commonly used internal reference genes were selected as the candidate internal reference genes for this experiment. They were GAPDH, UBQ, 18S rRNA, 28S rRNA, EF1a, UBC, Actin1, elF-3, and CYP (Table 1). According to the sequence of each gene, fluorescent quantitative PCR primers were designed using MEGA6.0. The primer sequence was synthesized by Qingke Bioengineering (Fuzhou, China) Co., Ltd. A five-fold cDNA dilution series with three replicates per concentration was used to make standard curves for estimation of amplification efficiency (E) and the correlation coefficient (R 2 ).  primer (10 µM), 0.4 µL downstream primer (10 µM), 1 µL DNA, 0.2 µL ROX reference dye (containing SYBR), and 8 µL ddH 2 O. The procedure of PCR reaction was as follows: 95 • C for 5 min, 40 cycles of 95 • C for 15 s, and 60 • C for 1 min, 95 • C for 15 s, 60 • C for 15 s, temperature gradually increased to 95 • C. Each reaction had three biological replicates and three technical replicates (i.e., each biological replicate was analyzed three times).

Validation of the Candidate Reference Genes
In order to verify the results of our experiments, the most stable reference genes were selected to validate the expression of the 14-3-3-4, ClPht1;3, and NRT1.1 genes in drought stress, P, and N starvation, respectively. The time and degree of experimental treatment were consistent with the three stress treatments mentioned above. The gene 14-3-3-4 belongs to 14-3-3 family, which binds to signal molecules to regulate plant response to drought stress [42]. The gene ClPht1;3 belongs to PHT1 family, which is the key transporter for uptake and transport of inorganic phosphorus from soil solution to the cell membranes of root. Moreover, the gene NRT1.1 belongs to NPF family and functions as NO 3 − transporter, and also was the substrates of polypeptides, amino acids, sugar ligands, and plant hormones [43].

Data Analysis
Quantitative real-time PCR was performed and analyzed using GoTagqPCR Master Mix (Promega) and a StepOnePlus™ real-time PCR instrument (Corbett Research, Australia). The expression stability of the nine candidate reference genes under different conditions were analyzed by geNorm, NormFinder, and BestKeeper, the most commonly used statistical procedures for screening reference genes. Expression levels of each reference gene were shown by Ct values. Before applying the geNorm and NormFinder analyses, the raw Ct values, weighted by amplification efficiency (Table 1), were used to calculate relative quantities by the equation: where E is the amplification efficiency and ∆Ct is the difference between the sample with the lowest Ct (highest expression) from the data set and the Ct value of the sample in question.
The geNorm analysis compares the stability of the candidate gene by calculating the average stability index, M, of the candidate gene. The stability boundary is usually set at M value equal to 1.5 and if the M value is lower than 1.5, the stability of the candidate gene will be higher. According to this principle, the geNorm analysis will step-wisely exclude the least stable gene and be repeated until only two genes remain. The lower the M value is the more stable gene expression will be [44][45][46]. NormFinder was used to evaluate the expression stability of the candidate reference genes on the experimental samples. Raw Ct values were first index-transformed and used as input for the NormFinder. This algorithm takes into account the intra-and inter-group variations for normalization factors (NFs), which requires input data from a minimum of three candidate reference genes and a minimum of two tested samples per group. The calculated results of this analysis were not influenced by random co-regulated genes and the best candidate reference gene displayed lower stability values that were close to zero. The more stable reference gene will have lower stability value and inter-and intra-group variation. In our case, we calculated the intra group variations as there are no distinct groups. BestKeeper analysis was applied to rank the stability of candidate genes by calculating coefficient of variance (CV). The candidate reference genes are considered to be the most stably expressed genes when they present the lowest CV. The relative expression level of the target genes, ClPht1;3, NRT1.1, 14-3-3-4, in Chinese fir roots were calculated by 2 −∆∆Ct method, which is applicable when the amplification efficiencies of the target and reference genes are approximately equal [47].

Expression Levels of Candidate Reference Genes
The expression levels (Ct values) of all nine reference genes are shown in Figure 1. The Ct values varied from 16.882 (28S rRNA) to 34.947 (CYP) across all samples, and mean Ct ranged from 25.757 (UBQ) to 29.517 (Actin1). Moreover, the expression levels of Actin1 were the most variable with 3.354 Ct, while 28S rRNA showed the least variable levels with 0.291 Ct. Considering the gene expression levels are negatively correlated with Ct values, the 28S rRNA had high expression level and CYP with low expression level. applicable when the amplification efficiencies of the target and reference genes are approximately equal [47].

Expression Levels of Candidate Reference Genes
The expression levels (Ct values) of all nine reference genes are shown in Figure 1. The Ct values varied from 16.882 (28S rRNA) to 34.947 (CYP) across all samples, and mean Ct ranged from 25.757 (UBQ) to 29.517 (Actin1). Moreover, the expression levels of Actin1 were the most variable with 3.354 Ct, while 28S rRNA showed the least variable levels with 0.291 Ct. Considering the gene expression levels are negatively correlated with Ct values, the 28S rRNA had high expression level and CYP with low expression level.

Drought Stress Experiment
The geNorm analysis revealed that GAPDH and EF1a genes had lowest M values in non-stressed (0% PEG6000) root samples, whereas the CYP genes had highest M values in non-stressed roots ( Figure 2A). In moderately stressed (10% PEG6000) root samples, elF-3 and CYP had the lowest M values, whereas the EF1a had the highest M value ( Figure 2B). In severely stressed roots (15% PEG6000), GAPDH and 18S rRNA had the lowest M values, whereas the CYP gene had the highest M value ( Figure 2C). Comprehensive analysis of M value of each gene under drought stress showed that GAPDH, 18S rRNA, and UBQ were the most stable reference genes, whereas the least stable gene was CYP.
The NormFinder also revealed differences in stability of expression of reference genes among drought stress treatments ( Table 2). The best reference gene in non-stressed (0% PEG6000) and moderately stressed (10% PEG6000) condition was UBQ, whereas the least stable gene in both treatments was CYP. For severely water stressed roots (15% PEG6000), the best reference gene was 18S rRNA, whereas the least stable gene was CYP. Overall, UBQ and GAPDH were the most stable genes, and CYP and elF-3 were the least stable genes for the drought stress sample sets.
The BestKeeper analysis also identified differences in stability of reference gene expression among different drought treatments (Table 3). BestKeeper analysis based on CV values identified UBQ and GAPDH as the two most stable reference genes in non-water stressed (0% PEG6000) and severely stressed (15% PEG6000) roots, and UBQ and 18S rRNA as highly stable genes in moderately (10% PEG6000) water stressed roots, respectively ( Table 3). The average expression level for drought stress sample sets revealed that UBQ and GAPDH were the best reference genes, and CYP and elF-3

Drought Stress Experiment
The geNorm analysis revealed that GAPDH and EF1a genes had lowest M values in non-stressed (0% PEG6000) root samples, whereas the CYP genes had highest M values in non-stressed roots ( Figure 2A). In moderately stressed (10% PEG6000) root samples, elF-3 and CYP had the lowest M values, whereas the EF1a had the highest M value ( Figure 2B). In severely stressed roots (15% PEG6000), GAPDH and 18S rRNA had the lowest M values, whereas the CYP gene had the highest M value ( Figure 2C). Comprehensive analysis of M value of each gene under drought stress showed that GAPDH, 18S rRNA, and UBQ were the most stable reference genes, whereas the least stable gene was CYP.
The NormFinder also revealed differences in stability of expression of reference genes among drought stress treatments ( Table 2). The best reference gene in non-stressed (0% PEG6000) and moderately stressed (10% PEG6000) condition was UBQ, whereas the least stable gene in both treatments was CYP. For severely water stressed roots (15% PEG6000), the best reference gene was 18S rRNA, whereas the least stable gene was CYP. Overall, UBQ and GAPDH were the most stable genes, and CYP and elF-3 were the least stable genes for the drought stress sample sets.
The BestKeeper analysis also identified differences in stability of reference gene expression among different drought treatments (Table 3). BestKeeper analysis based on CV values identified UBQ and GAPDH as the two most stable reference genes in non-water stressed (0% PEG6000) and severely stressed (15% PEG6000) roots, and UBQ and 18S rRNA as highly stable genes in moderately (10% PEG6000) water stressed roots, respectively ( Table 3). The average expression level for drought stress sample sets revealed that UBQ and GAPDH were the best reference genes, and CYP and elF-3 were the least stable genes, which is similar to the results obtained using geNorm and NormFinder analysis. were the least stable genes, which is similar to the results obtained using geNorm and NormFinder analysis.      Figure 2D-F). According to the M value of each gene in all samples under phosphorus stress, the most stable internal reference genes were Actin1, elF-3, and GAPDH. EF1a was the least stable gene.
According to the NormFinder analysis, the best reference gene was Actin1 in P-starved and low P supply roots, whereas 28S rRNA and EF1a genes were the least stable genes, respectively. In normal P supply root, the best reference gene was CYP, whereas the least stable gene was EF1a. Overall, Actin1 and GAPDH were the most stable genes, and EF1a and 28S rRNA were the least stable genes for the P starvation sample sets (Table 4). BestKeeper analysis also identified GAPDH as the best reference gene for P-starved and low P root, and Actin1 as highly stable genes in normal P supplied root ( Table 5). As a whole, Actin1 and GAPDH were the most stable genes and EF1a and UBQ were the least stable genes for the P starvation sample sets.

N Starvation Experiment
In the N starvation experiment, the EF1a and Actin1 were the most stable genes in N-starved and low N supply root samples, whereas the 18S rRNA and UBC were the most stably expressed genes in roots supplied with normal N level based on their low M values ( Figure 2G-I). The least stable genes were 28S rRNA, GAPDH, and CYP in N-starved roots, roots supplied with low and normal N levels, respectively. Comparison between N sources revealed that 18S rRNA and 28S rRNA genes were stable in NO 3 − stress samples ( Figure 2J) and the UBC and Actin1 were the most stable genes in NH 4 + stress samples ( Figure 2K). The least stable genes were CYP and 18S rRNA in NO 3 − and NH 4 + stress samples based on their high M values, respectively. According to the M value of each gene in all samples under nitrogen stress, the most stable internal reference genes were Actin1 and UBC, and CYP and EF1a were the most unstable genes. NormFinder analysis revealed that Actin1 was the best reference gene in N-starved, normal N supplied, and NH 4 + stressed roots, whereas the least stable gene in N-starved root was 28S rRNA and 18S rRNA in NH 4 + stressed roots. The expression of CYP was the least stable in normal N and NO 3 − supplied root, whereas the best stable gene in NO 3 − supplied root was UBC. In low N supplied root, the best reference gene was UBC, whereas the least stable gene was GAPDH. Overall, UBC and Actin1 were the most stable genes and 28S rRNA and CYP were the least stable genes for the N starvation sample sets (Table 6). BestKeeper analysis also identified Actin1 as the most stable reference gene for low N-, normal N-, and NO 3 − -N-treated root, and CYP as highly stable genes for N-starved and NH 4 + -N supplied roots ( Table 7). As a whole, Actin1 and UBC were the most stable genes and 28S rRNA and EF1a were the least stable genes for N starvation sample sets. In most tested samples, geNorm analysis results were roughly the same as NormFinder and BestKeeper analyses, with slight variations in the ranking sequence of genes. In addition, regardless of the variations in ranking, these analyses identified the same most stably expressed gene in all experimental samples.

Validation of Stability of Reference Genes
To verify the reliability of selected reference genes, expression profiles of 14-3-3-4, ClPht1;3, and NRT1.1 genes were determined in drought, P, and N stress experiments, respectively. Relative expression levels were normalized using two most stable reference genes (UBQ and GAPDH, Actin1 and GAPDH, Actin1 and UBC) and the least stable gene (CYP, EF1a, 28S rRNA) in each abiotic stress experiment. The expression of 14-3-3-4 gene showed similar level when a single reference gene was used to normalize the expression. When UBQ and GAPDH were used as reference genes, the expression of 14-3-3-4 gene was up regulated with increasing drought stress level without a marked difference among different drought stress times ( Figure 3). However, when CYP (unstable gene) was used for normalization of the target gene, relative expression profile of 14-3-3-4 was different between stress periods compared to normalization of expression done using the two most stable reference genes identified in our study (UBQ and GAPDH). Because the expression of 14-3-3-4 gene was up regulated by drought stress [48], UBQ and GAPDH could be used as reference genes to obtain relatively stable results.
When Actin1 was used as internal reference gene for P starvation experiment, the expression profile of ClPht1;3 was the highest under P deficiency after 1, 3, and 7 days (Figure 4). After 15 days of P treatment, the expression of the target gene (ClPht1;3) slightly decreased in P deficient roots. However, when EF1a (unstable gene) was used for normalization of the target gene, relative expression profile of ClPht1;3 was different compared to normalization of expression done using the two most stable reference genes identified in our study (Actin1 and GAPDH).
difference among different drought stress times (Figure 3). However, when CYP (unstable gene) was used for normalization of the target gene, relative expression profile of 14-3-3-4 was different between stress periods compared to normalization of expression done using the two most stable reference genes identified in our study (UBQ and GAPDH). Because the expression of 14-3-3-4 gene was up regulated by drought stress [48], UBQ and GAPDH could be used as reference genes to obtain relatively stable results. UBQ and GAPDH were used as two most stable reference genes, CYP was used as the least stable reference gene. PEG0%: 0% PEG6000 treatment; PEG10%: 10% PEG6000 treatment; PEG15%: 15% PEG6000 treatment; PEG: Polyethylene Glycol treatment of molecular weight 6000.
When Actin1 was used as internal reference gene for P starvation experiment, the expression profile of ClPht1;3 was the highest under P deficiency after 1, 3, and 7 days (Figure 4). After 15 days of P treatment, the expression of the target gene (ClPht1;3) slightly decreased in P deficient roots. However, when EF1a (unstable gene) was used for normalization of the target gene, relative expression profile of ClPht1;3 was different compared to normalization of expression done using the two most stable reference genes identified in our study (Actin1 and GAPDH). UBQ and GAPDH were used as two most stable reference genes, CYP was used as the least stable reference gene. PEG0%: 0% PEG6000 treatment; PEG10%: 10% PEG6000 treatment; PEG15%: 15% PEG6000 treatment; PEG: Polyethylene Glycol treatment of molecular weight 6000.  Actin1 and GAPDH were used as two most stable reference genes. EF1a was used as the least stable reference gene. NO P: P-starved supply treatment; LOW P: low P supply treatment; NORMAL P: normal P supply treatment.
NRT1.1 participated in the absorption of nitrate by plant roots with low affinity, involving transport of different nitrogen concentrations [20]. When Actin1 and UBC were used as reference genes, the relative expression of NRT1.1 was basically the same with the elongation of time at the same concentration ( Figure 5). At the same time, the expression level of NRT1.1 was up regulated with the increase in concentration of N and the highest expression level was observed in the treatment of NO3 − and the lowest was in the treatment of NH4 + . The relative expression of the target gene (NRT1.1) was significantly different from the former two when least stable gene (28S rRNA) was used for normalization of expression of target gene. This indicated that the expression stability of 28S in Figure 4. Relative expression of ClPht1;3 under P stress treatments at different time period. Actin1 and GAPDH were used as two most stable reference genes. EF1a was used as the least stable reference gene. NO P: P-starved supply treatment; LOW P: low P supply treatment; NORMAL P: normal P supply treatment.
NRT1.1 participated in the absorption of nitrate by plant roots with low affinity, involving transport of different nitrogen concentrations [20]. When Actin1 and UBC were used as reference genes, the relative expression of NRT1.1 was basically the same with the elongation of time at the same concentration ( Figure 5). At the same time, the expression level of NRT1.1 was up regulated with the increase in concentration of N and the highest expression level was observed in the treatment of NO 3 − and the lowest was in the treatment of NH 4 + . The relative expression of the target gene (NRT1.1) was significantly different from the former two when least stable gene (28S rRNA) was used for normalization of expression of target gene. This indicated that the expression stability of 28S in different nitrogen treated roots was poor and it was not suitable for accurate correction of the relative expression of NRT1.1. As a whole, the results showed that different reference genes were used to correct the expression of target genes in fluorescence quantitative test and different results would be obtained. If the reference genes were not selected properly, the relative expression of target genes could be misestimated.

Discussion
The qRT-PCR is widely used in gene expression analysis because of its high specificity and sensitivity. qRT-PCR analysis results are closely related to the selection of reference genes and improper selection of reference genes may lead to deviation of analysis results and even wrong conclusions [48][49][50][51]. Therefore, to ensure the validity of the gene expression analysis results, it is necessary to screen the reference genes before using qRT-PCR to analyze the expression of the target genes. A large number of studies have shown that the screening of reference genes should be accurately evaluated and validated according to plant type [39], tissue location [32], and experimental conditions [21], so as to establish an evaluation system of multiple reference genes, which is an important step towards improving the accuracy of qRT-PCR test [52,53].
In this study, nine commonly used candidate reference genes (GAPDH, UBQ, 18S, 28S, EF1a, UBC, Actin1, elF-3, and CYP) were screened from fine roots of Chinese fir under three abiotic stresses (drought, P, and N stress treatments) at different time and concentrations. In a large number of papers, the stability of these commonly used reference genes is different in different plants, different

Discussion
The qRT-PCR is widely used in gene expression analysis because of its high specificity and sensitivity. qRT-PCR analysis results are closely related to the selection of reference genes and improper selection of reference genes may lead to deviation of analysis results and even wrong conclusions [48][49][50][51]. Therefore, to ensure the validity of the gene expression analysis results, it is necessary to screen the reference genes before using qRT-PCR to analyze the expression of the target genes. A large number of studies have shown that the screening of reference genes should be accurately evaluated and validated according to plant type [39], tissue location [32], and experimental conditions [21], so as to establish an evaluation system of multiple reference genes, which is an important step towards improving the accuracy of qRT-PCR test [52,53].
In this study, nine commonly used candidate reference genes (GAPDH, UBQ, 18S, 28S, EF1a, UBC, Actin1, elF-3, and CYP) were screened from fine roots of Chinese fir under three abiotic stresses (drought, P, and N stress treatments) at different time and concentrations. In a large number of papers, the stability of these commonly used reference genes is different in different plants, different tissues, different developmental stages and different experimental treatments [42,54]. Our study also showed that the reference genes of Chinese fir roots under different stress treatments were not the same, and the stability of gene expression is highly dependent on external conditions. We identified different stable reference genes for drought stress and nutrient starvation, which effectively normalize the expression of target genes. A previous study screened and validated the stability of five housekeeping genes (Actin, 18S, UBQ, EF1a, and GAPDH) in root and leaf tissues of Chinese fir under different abiotic stresses and found Actin and GAPDH as stable reference genes for root and leaf tissues of Chinese fir, respectively [41]. We also found GAPDH and UBQ as the two most stable genes under drought stress, Actin1 and GAPDH as the two most stable genes under P starvation, and Actin1 and UBC as the two most stable genes under N starvation among the nine candidate reference genes. Thus, not all housekeeping genes are equally valuable to serve as reference genes to normalize the expression level of target genes, as reported previously [55].
To ensure accurate normalization, some authors have recommended that multiple reference genes be used in gene expression analysis [18,[56][57][58]. In this experiment, we selected the best two reference genes and the worst one to verify the expression levels of genes related to drought and nutrient stresses. The expression levels of genes related to three abiotic stresses were verified by examining the expression level of 14-3-3-4, ClPht1;3 and NRT1.1 in drought stressed, P, and N stressed roots, respectively. According to the verification results of target gene expression, UBQ and GAPDH were selected as single standardized reference genes under drought stress; Actin1 and GAPDH were selected as single standardized reference genes under P starvation; and Actin1 and UBC were selected as single standardized reference genes under N starvation. The expression patterns of 14-3-3-4, ClPht1;3 and NRT1.1 were similar under the same abiotic stress. The result also showed that the expression pattern of 14-3-3-4, ClPht1;3, and NRT1.1 were almost the same as that of one or two reference genes after normalization. Thus, our results suggest that UBQ gene can be used to standardize Chinese fir root under drought stress, while the Actin gene is the best internal control gene of Chinese fir root under N and P starvation. The stable reference gene selected in this study will be very useful for revealing the gene expression profiles in Chinese fir and related conifer species under abiotic stress; promoting the realization of it at cellular and gene levels. As abiotic stresses, such as drought and low nutrient availability, are major factors limiting tree growth, understanding the genetic basis for coping with such stresses has paramount importance in identifying stress tolerant genotypes. Genes that are up or down regulated during stress have not been fully understood in trees, which necessities further gene expression studies. For gene expression analysis to be accurate, it needs to be normalized using stable reference genes. In this regard, our study will lay the foundation for further investigation of stress-tolerance mechanisms at the molecular level in conifers by providing stable reference genes for normalization of gene expression levels in qRT-PCR analysis.

Conclusions
In summary, this study used qRT-PCR technology to determine the optimal reference genes for the relative quantification of transcript abundance in Chinese fir under different abiotic stresses. The expression stability of nine candidate reference genes in Chinese fir root samples was tested under three abiotic treatments with different time intervals and concentrations. As a consequence, we recommend UBQ and the Actin1 gene for drought, P, and N stresses as superior internal controls for normalization of qRT-PCR, respectively. Additionally, our results showed that different suitable reference genes or standardized reference gene combinations could be screened according to different external conditions. The findings reported here will facilitate further studies on stress ecology and gene function analysis in Chinese fir. Our results also provide a basis for normalization of gene expression in other coniferous species, such as Douglas fir.