1. Introduction
Chinese fir (Cunninghamia lanceolata (Lamb.) Hook.) is an evergreen conifer, which has been widely planted for timber in southeastern Asia for more than 1000 years [1]. It is the most commercially important afforestation tree species in southern China with fast growth [2], high yield, and good material quality and pest resistance [3,4]. Because of its fast-growing characteristics, Chinese fir has a high nutrient requirement for its optimal growth, especially nitrogen and phosphorus, which are scarce in forest soils in southern China [5,6]. Furthermore, successive rotation of Chinese fir plantations amplifies the lack of available N and P in forest soil due to clear-cutting and subsequent burning practice [7]. It had been found that the yield of Chinese fir plantation decreased by 5–15% with successive planting in recent years and the N and P deficiency was the most important factor limiting growth of the species. Thus, application of N and P fertilizers has been commonly used in Chinese fir plantations, which is costly for tree growers and causing serious issues of soil and water pollution. In addition, the extremely harsh weather conditions typical of southern China coupled with seasonal droughts due to global climate change seriously affect the survival of seedlings at the beginning of reforestation [8]. Most previous studies on Chinese fir focused on morphological and physiological adaptation mechanisms to nutrient stress including root morphological plasticity [9], root secretion of organic acids [10], formation of cortical aerenchyma [11], and physiological plasticity and biochemical changes [12]. However, studies on gene expression and molecular biology in Chinese fir are scarce. Thus, it is essential to understand the mechanisms underlying the response of Chinese fir to abiotic stresses using gene expression analysis in order to develop stress tolerant genotypes.
Transcriptome analysis is mainly used to detect the gene expression profiles and quantitative real-time reverse transcription-polymerase chain reaction (qRT-PCR) is a stable method with the advantages of high precision, accuracy, and specificity for detecting target gene expression [13]. It has become one of the most powerful tools used in molecular biology research [14]. However, the results of qRT-PCR are easily prone to interference with experimental errors such as RNA quality, sample loading, and reverse transcription efficiency. Therefore, it requires internal reference genes to correct these experimental errors between samples [15]. The ideal reference genes are thought to be equally expressed in different samples, developmental stages, and hence they can be used to measure the expressions of the target genes [16]. The accuracy and reliability of qRT-PCR results is greatly dependent on a suitable reference gene to normalize gene expression and avoid errors caused by different experimental procedure [17,18]. The more stable the expression of the internal reference gene is, the better the accuracy of expression level of the target gene will be [18]. Usually, housekeeping genes like Actin, EF1a (elongation factor 1-a), GAPDH (glyceraldehyde-3-phosphate dehydrogenase), 18S rRNA (18S ribosomal RNA) etc. are used as reference genes [19,20] but previous studies have reported that the stability of expression of these reference genes is relative, and the stability of expression of housekeeping genes varies considerably in different plant species, different tissues, and different experimental conditions [15,21]. Particularly in some non-model plant species, like Chinese fir, for which there is a lack of genome sequence and poor genetic background information, the reference genes are usually used by the orthologous sequence of common housekeeping genes reported in model plant species [22], which is the probable case of deviation of the accuracy in qRT-PCR experiment [23,24]. Therefore, it is important to screen and evaluate the stability of reference gene expression according to experimental conditions to improve the accuracy of qRT-PCR analyses [25,26,27].
Due to the importance of reference gene selection, there are an increasing number of attempts to identify reference genes in different experimental conditions. Nicot et al. [28] analyzed the expression stability of seven housekeeping genes of Actin, APRT (Adeninephosphorribosyltransferase), 18S rRNA, EF1a, TUB (Tubulin), CYP (cytochrome), and Ribosomal protein L2 under biological stress (late blight) and abiotic stress (salt, drought) in potato, and found variation in the expression level of target gene hsp20.2 according to the standard of reference genes with different stability. Hong et al. [29] analyzed the expression stability of nine reference genes in 21 tissues and organs at different developmental stages under abiotic stresses (drought, salt, cold, heat) and hormone conditions in Brachpodium distachyon, and reported that UBC18 (ubiquitin-conjugating enzyme) was the most stable expression in all samples, and the expression stability of UBQ4 (ubiquitin) and UBQ10 were the best in different tissues, organs and developmental stages, and SamDC was the most stable gene in environmental stress samples. Additionally, reference gene validation has been done in many plant species, such as Gossypium barbadense [30], Saccharum officinarum [31], Lycopersicon esculentum [32], Vitis viniferag [33], Glycine max [34], Musa acuminata [35], Oryza sativa [36], Linum usitatissimum [37], Kosteletzkya virginica [20], Sesamum indicum [38], and Populus euphratica [39]. However, few studies have been conducted on the selection and expression stability of reference genes in Chinese fir, such as selection of reference genes under cold, high temperature, abscisic acid (ABA), polyethylene glycol (PEG), and sodium chloride (NaCl), and the selected reference was different for different experimental conditions [40,41].
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.
2. Materials and Methods
2.1. 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 KH2PO4, 5 mmol/L KNO3, 2 mmol/L MgSO4-7H2O, 5 mmol/L Ca(NO3)2-H2O, 1 mmol/L Na2-EDTA, 0.001 mmol/L FeSO4-7H2O, and micro nutrients including 0.463 mmol/L H3BO3, 0.003 mmol/L CuSO4-5H2O, 0.007 mmol/L ZnSO4-7H2O, 0.091 mmol/L MnCl2-4H2O, 0.005 mmol/L H8MoN2O4, 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 KH2PO4, 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 (NH4)2SO4, 2 mmol/L KNO3, and 1 mmol/L NH4NO3 were used to represent NH4+, NO3−, 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 NH4NO3, 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.
2.2. 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 A260/A280 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 GoScriptTM 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.
2.3. 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 (R2).
2.4. Quantitative Real-Time PCR (qRT-PCR)
Real-time fluorescence quantitative PCR was performed on the Thermo Fisher’s StepOnePlus real-time fluorescence quantitative PCR instrument. Promega’sGoTagqPCR Master Mix kit was used for PCR. The reaction system was 20 µL, including: 10 µL 2 × qPCR Master Mix, 0.4 µL upstream primer (10 µM), 0.4 µL downstream primer (10 µM), 1 µL DNA, 0.2 µL ROX reference dye (containing SYBR), and 8 µL ddH2O. 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).
2.5. 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 NO3− transporter, and also was the substrates of polypeptides, amino acids, sugar ligands, and plant hormones [43].
2.6. 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].
3. Results
3.1. 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.
3.2. Stability of Candidate Genes
3.2.1. 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.
3.2.2. P Starvation Experiment
In the P starvation experiment, UBQ and UBC genes had lowest M values in P-starved root samples while 28S rRNA, 18S rRNA, and EF1a had highest M value. GAPDH and 18S rRNA had the lowest M values in low P-stressed roots while EF1a and elF-3 had the highest M values. Actin1 and elF-3 had the lowest M value in P-replete root samples while EF1a and UBQ had the highest M values (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.
3.2.3. 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 NO3− stress samples (Figure 2J) and the UBC and Actin1 were the most stable genes in NH4+ stress samples (Figure 2K). The least stable genes were CYP and 18S rRNA in NO3− and NH4+ 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 NH4+ stressed roots, whereas the least stable gene in N-starved root was 28S rRNA and 18S rRNA in NH4+ stressed roots. The expression of CYP was the least stable in normal N and NO3− supplied root, whereas the best stable gene in NO3− 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 NO3−-N-treated root, and CYP as highly stable genes for N-starved and NH4+-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.
3.3. 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).
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 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.
4. 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.
5. 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.
Author Contributions
Conceptualization, X.M. and M.L.; Formal analysis, R.C. and W.C.; Funding acquisition, M.L.; Investigation, R.C., W.C., W.Z., and Y.L.; Methodology, R.C., W.C., W.Z., Y.L., X.M., and M.L.; Supervision, M.L.; Validation, M.T.; Visualization, R.C.; Writing—original draft, R.C.; Writing—review and editing, M.T.
Funding
This work was financially supported by the National Natural Science Foundation of China (31971674).
Acknowledgments
We are grateful for Weimin Zhong, Wei Liu, Shaoliang Wang and Peng Shuai for their assistance during the lab work.
Conflicts of Interest
The authors declare no conflict of interest.
Figures and Tables
Enlarge this image.Figure 1. Expression levels (Ct values) of the nine candidate reference genes. The lines across the box indicate median values. The upper and lower edges bordered in each box are interquartile range, which indicate the 25th and 75th percentiles.
Enlarge this image.Figure 2. Average expression stability (M values) and ranking of candidate reference genes in Chinese fir root with different stress treatment calculated by geNorm. A lower M value indicates a more stable expression. (A): 0%-Root: 0% PEG6000 treatment; (B): 10%-Root: 10% PEG6000 treatment; (C): 15%-Root: 15% PEG6000 treatment; (D): NO P-Root: P-starved supply treatment; (E): Low P-Root: low P supply treatment; (F): Normal P-Root: normal P supply treatment; (G): NO N-Root: N-starved supply treatment; (H): Low N-Root: low N supply treatment; (I): Normal N-Root: normal N supply treatment; (J): NO3−-N-Root: NO3− stress treatment; (K): NH4+-N-Root: NH4+ stress treatment.
Enlarge this image.Figure 3. Relative expression of 14-3-3-4 under drought stress for different time period. 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.
Enlarge this image.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.
Enlarge this image.Figure 5. Relative expression of NRT1.1 under N stress treatments at different time period. Actin1 and UBC were used as two most stable reference genes, 28S rRNA was used as the least stable reference gene. NO N: N-starved supply treatment; LOW N: low N supply treatment; NORMAL N: normal N supply treatment; NO3− N: NO3− stress treatment; NH4+-N: NH4+ stress treatment.
Primer sequences of the candidate reference genes and target genes (ClPht1;3, NRT1.1, 14-3-3-4) together with polymerase chain reaction (PCR) amplification efficiency and the regression coefficient.
| Gene Name | Primer Sequence (5′ to 3′) | TM (°C) | GC Content (%) | Amplification Efficiency (%) | Reg. Coeff. (R2) |
|---|---|---|---|---|---|
| GAPDH | F: GCACCTATGTTTGTGGTTGGAGTA | 60.3 | 45.8 | 102.018 | 0.976 |
| R: ACCGTCTTCTGTGTAGCTGTTGTT | 60.3 | 45.8 | |||
| UBQ | F: AATAAATGCTTCAAATGTCAGGCTA | 55.4 | 32.0 | 95.653 | 0.987 |
| R: TGAGATGGTCTGGTGATGTCGTGG | 63.7 | 54.2 | |||
| 18S rRNA | F: GGGAACATTATCACGGACAGCATCAAC | 63.5 | 48.1 | 101.547 | 0.994 |
| R: TGCGACTAATGGTCCAGATAGACTCCT | 63.5 | 48.1 | |||
| 28S rRNA | F: AGTTTGTAGAAGCGTCCTCAG | 58.0 | 47.6 | 103.030 | 0.988 |
| R: GGCGCCCCCTTCCA | 58.8 | 78.6 | |||
| EF1a | F: TGGCAAGGAGCTTGAGAAAGAACCCA | 63.6 | 50.0 | 98.927 | 0.994 |
| R: ACCCCAACAGCAACAGTCTGACGCAT | 65.1 | 53.8 | |||
| UBC | F: CTCCTGATGATAGAGGTCTAACG | 60.2 | 47.8 | 96.603 | 0.991 |
| R: ATGACATGGAGGTCGGTACATAG | 60.2 | 47.8 | |||
| Actin1 | F: CTCTCTCAGCACCTTCGAGCAG | 63.8 | 59.1 | 100.458 | 0.965 |
| R: TCCACATACAACCGCTCCACTG | 61.9 | 54.5 | |||
| elF-3 | F: TCTGCATCCCAATGATTTGT | 53.7 | 53.7 | 94.829 | 0.999 |
| R: TCTTCAGATGTTGCTTGCCT | 55.8 | 45.0 | |||
| CYP | F: TGTGTCAGGGTGGCGACTTCA | 61.9 | 57.1 | 94.931 | 0.953 |
| R: TTCGCCATGGAGAGGATGCC | 61.9 | 60.0 | |||
| ClPht1;3 | F: CCAATACAACAACTTTCGTCGTG | 58.4 | 43.5 | 93.960 | 0.951 |
| R: GTCCTCGTGCTTATTCTGATTTG | 58.4 | 43.5 | |||
| NRT1.1 | F: GCCAAGAACACCGTGAACAT | 57.8 | 50 | 98.658 | 0.950 |
| R: TCAGCCCGAAGAAGTGATGT | 57.8 | 50 | |||
| 14-3-3-4 | F: CATGCCATGGATGAGTACCGAGAAG | 63.6 | 52 | 94.995 | 0.944 |
| R: CGGGATCCTCACATCTCTTCATCT | 62.0 | 50 |
Expression stability values and ranking of reference genes for drought experiment calculated by the NormFinder.
| Intra Group | ||||||||
|---|---|---|---|---|---|---|---|---|
| Rank | 0%-Root | 10%-Root | 15%-Root | Overall | ||||
| Gene | Stab | Gene | Stab | Gene | Stab | Gene | Stab | |
| 1 | UBQ | 0.185 | UBQ | 0.367 | 18S | 0.221 | UBQ | 0.907 |
| 2 | GAPDH | 0.394 | 28S | 0.900 | GAPDH | 0.343 | GAPDH | 1.638 |
| 3 | EF1a | 0.563 | GAPDH | 0.901 | UBQ | 0.355 | Actin1 | 2.226 |
| 4 | UBC | 0.625 | UBC | 1.134 | 28S | 0.370 | 18S | 2.258 |
| 5 | Actin1 | 0.670 | Actin1 | 1.177 | EF1a | 0.371 | EF1a | 2.397 |
| 6 | 18S | 0.786 | CYP | 1.218 | Actin1 | 0.379 | UBC | 2.416 |
| 7 | elF-3 | 1.046 | 18S | 1.251 | UBC | 0.657 | 28S | 2.703 |
| 8 | 28S | 1.433 | elF-3 | 1.391 | elF-3 | 0.827 | elF-3 | 3.264 |
| 9 | CYP | 1.479 | CYP | 1.463 | CYP | 1.531 | CYP | 4.228 |
| Best gene | UBQ | UB Q | 18S | UBQ | ||||
| 0.185 | 0.367 | 0.221 | 0.907 | |||||
0%-Root: 0% PEG6000 treatment; 10%-Root: 10% PEG6000 treatment; 15%-Root: 15% PEG6000 treatment; PEG: Polyethylene glycol treatment of molecular weight 6000.
Table 3Expression stability values and ranking of reference genes for drought experiment calculated by BestKeeper.
| Rank | 0%-Root | 10%-Root | 15%-Root | Average | ||||
|---|---|---|---|---|---|---|---|---|
| Gene | CV | Gene | CV | Gene | CV | Gene | CV | |
| 1 | UBQ | 5.85 | UBQ | 5.27 | UBQ | 4.00 | UBQ | 5.04 |
| 2 | GAPDH | 7.32 | 18S | 5.29 | GAPDH | 4.35 | GAPDH | 5.92 |
| 3 | EF1a | 8.03 | GAPDH | 6.10 | 18S | 5.44 | EF1a | 6.92 |
| 4 | UBC | 8.36 | 28S | 6.57 | EF1a | 5.79 | 18S | 6.96 |
| 5 | 28S | 9.58 | EF1a | 6.95 | 28S | 8.12 | UBC | 8.01 |
| 6 | Actin1 | 9.98 | UBC | 7.35 | UBC | 8.31 | 28S | 8.09 |
| 7 | 18S | 10.15 | Actin1 | 8.45 | Actin1 | 8.84 | Actin1 | 9.09 |
| 8 | elF-3 | 12.19 | elF-3 | 10.37 | elF-3 | 11.84 | elF-3 | 11.47 |
| 9 | CYP | 14.31 | CYP | 10.99 | CYP | 16.96 | CYP | 14.09 |
Expression stability values and ranking of reference genes for phosphorus experiment calculated by the NormFinder.
| Intra Group | ||||||||
|---|---|---|---|---|---|---|---|---|
| Rank | No P-Root | Low P-Root | Normal P-Root | Overall | ||||
| Gene | Stab | Gene | Gene | Stab | Gene | Gene | Stab | |
| 1 | Actin1 | 0.773 | Actin1 | 0.079 | CYP | 0.700 | Actin1 | 2.319 |
| 2 | CYP | 0.776 | GAPDH | 0.079 | Actin1 | 1.467 | GAPDH | 2.422 |
| 3 | GAPDH | 0.800 | 28S | 0.534 | GAPDH | 1.543 | CYP | 2.512 |
| 4 | elF-3 | 0.915 | CYP | 1.036 | elF-3 | 1.553 | elF-3 | 3.879 |
| 5 | UBC | 1.069 | UBC | 1.077 | 18S | 1.891 | UBC | 4.133 |
| 6 | UBQ | 1.083 | 18S | 1.159 | UBC | 1.987 | UBQ | 4.686 |
| 7 | EF1a | 1.713 | UBQ | 1.313 | 28S | 2.123 | 18S | 4.861 |
| 8 | 18S | 1.811 | elF-3 | 1.411 | UBQ | 2.290 | 28S | 2.512 |
| 9 | 28S | 2.737 | EF1a | 3.048 | EF1a | 2.726 | EF1a | 7.487 |
| Best gene | Actin1 | Actin1 | CYP | Actin1 | ||||
| 0.773 | 0.079 | 0.700 | 2.319 | |||||
NO P-Root: P-starved supply treatment; Low P-Root: low P supply treatment; Normal P-Root: normal P supply treatment.
Table 5Expression stability values and ranking of reference genes for phosphorus experiment calculated by BestKeeper.
| Rank | No P-Root | Low P-Root | Normal P-Root | Average | ||||
|---|---|---|---|---|---|---|---|---|
| Gene | CV | Gene | CV | Gene | CV | Gene | CV | |
| 1 | GAPDH | 2.28 | GAPDH | 0.59 | Actin1 | 0.97 | Actin1 | 1.32 |
| 2 | Actin1 | 2.39 | Actin1 | 0.61 | GAPDH | 2.36 | GAPDH | 1.74 |
| 3 | EF1a | 3.36 | 18S | 0.63 | 28S | 2.85 | CYP | 5.40 |
| 4 | CYP | 5.80 | UBC | 2.11 | elF-3 | 3.07 | 18S | 6.78 |
| 5 | elF-3 | 8.71 | UBQ | 5.50 | CYP | 4.28 | elF-3 | 7.16 |
| 6 | UBC | 9.77 | CYP | 6.11 | 18S | 7.89 | UBC | 7.56 |
| 7 | UBQ | 10.29 | 28S | 8.41 | UBQ | 8.17 | 28S | 7.82 |
| 8 | 18S | 11.81 | elF-3 | 9.70 | UBC | 10.80 | UBQ | 7.99 |
| 9 | 28S | 12.20 | EF1a | 17.47 | EF1a | 17.58 | EF1a | 12.80 |
Expression stability values and ranking of reference genes for nitrogen experiment calculated by the NormFinder.
| Rank | No N-Root | Low N-Root | Normal N-Root | NO3−-N-Root | NH4+-N-Root | Overall | ||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Gene | Stab | Gene | Stab | Gene | Stab | Gene | Stab | Gene | Stab | Gene | Stab | |
| 1 | Actin1 | 0.048 | UBC | 0.052 | Actin1 | 0.230 | UBC | 0.204 | Actin1 | 0.077 | UBC | 0.638 |
| 2 | EF1a | 0.060 | Actin1 | 0.063 | UBC | 0.230 | elF-3 | 0.224 | UBC | 0.104 | Actin1 | 0.774 |
| 3 | GAPDH | 0.094 | elF-3 | 0.082 | 18S | 0.243 | Actin1 | 0.356 | GAPDH | 0.117 | elF-3 | 1.594 |
| 4 | UBC | 0.121 | EF1a | 0.107 | elF-3 | 0.248 | UBQ | 0.609 | UBQ | 0.170 | UBQ | 1.594 |
| 5 | UBQ | 0.162 | 18S | 0.117 | UBQ | 0.308 | EF1a | 0.686 | 28S | 0.322 | GAPDH | 1.964 |
| 6 | elF-3 | 0.234 | 28S | 0.232 | GAPDH | 0.329 | 18S | 0.687 | CYP | 0.527 | EF1a | 2.354 |
| 7 | CYP | 0.341 | UBQ | 0.345 | 28S | 0.570 | 28S | 0.715 | elF-3 | 0.806 | 18S | 2.554 |
| 8 | 18S | 0.446 | CYP | 0.385 | EF1a | 0.582 | GAPDH | 0.815 | EF1a | 0.919 | CYP | 3.416 |
| 9 | 28S | 1.664 | GAPDH | 0.609 | CYP | 0.620 | CYP | 1.543 | 18S | 1.061 | 28S | 3.503 |
| Best gene | Actin1 | UBC | Actin1 | UBC | Actin1 | UBC | ||||||
| 0.048 | 0.052 | 0.230 | 0.204 | 0.077 | 0.638 | |||||||
Expression stability values and ranking of reference genes for nitrogen experiment calculated by BestKeeper.
| Rank | No N-Root | Low N-Root | Normal N-Root | NO3−-N-Root | NH4+-N-Root | Average | ||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Gene | CV | Gene | CV | Gene | CV | Gene | CV | Gene | CV | Gene | CV | |
| 1 | CYP | 0.33 | Actin1 | 0.36 | Actin1 | 1.02 | Actin1 | 0.73 | CYP | 0.28 | Actin1 | 0.55 |
| 2 | Actin1 | 0.36 | CYP | 1.48 | 18S | 1.80 | UBQ | 0.93 | UBC | 0.35 | UBC | 1.29 |
| 3 | UBC | 0.84 | UBC | 1.49 | UBC | 2.26 | UBC | 1.52 | Actin1 | 0.65 | 18S | 1.46 |
| 4 | 18S | 0.97 | 18S | 1.49 | UBQ | 2.30 | 28S3 | 1.78 | GAPDH | 1.06 | UBQ | 1.58 |
| 5 | UBQ | 1.05 | EF1a | 1.51 | elF-3 | 2.76 | 18S | 1.81 | 18S | 1.25 | CYP | 2.29 |
| 6 | EF1a | 1.19 | elF-3 | 1.99 | GAPDH | 2.94 | elF-3 | 2.04 | UBQ | 1.36 | GAPDH | 2.31 |
| 7 | GAPDH | 1.54 | UBQ | 2.28 | 28S | 3.63 | GAPDH | 3.18 | elF-3 | 3.41 | elF-3 | 2.418 |
| 8 | elF-3 | 1.89 | GAPDH | 2.81 | CYP | 3.99 | EF1a | 3.83 | 28S | 3.54 | EF1a | 2.91 |
| 9 | 28S | 2.87 | 28S | 3.23 | EF1a | 4.46 | CYP | 5.37 | EF1a | 3.55 | 28S | 2.96 |
© 2019 by the authors.
Details
; Zhong, Weimin 1 ; Li, Yushan 1 ; Ma, Xiangqing 3 ; Li, Ming 3 1 Forestry College, Fujian Agriculture and Forestry University, Fuzhou 350002, China;
2 Forestry College, Fujian Agriculture and Forestry University, Fuzhou 350002, China;
3 Forestry College, Fujian Agriculture and Forestry University, Fuzhou 350002, China;