Next Article in Journal
Assessing Forest Structure and Composition along the Altitudinal Gradient in the State of Sikkim, Eastern Himalayas, India
Next Article in Special Issue
Non-Targeted Metabolomics Reveals Patterns of Metabolic Changes during Poplar Seed Germination
Previous Article in Journal
Space, Habitat and Isolation are the Key Determinants of Tree Colonization by the Carpenter Ant in Plantation Forests
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Selection of Suitable Reference Genes in Pinus massoniana Lamb. Under Different Abiotic Stresses for qPCR Normalization

Co-Innovation Center for Sustainable Forestry in Southern China, Nanjing Forestry University, Nanjing 210037, China
*
Author to whom correspondence should be addressed.
Forests 2019, 10(8), 632; https://doi.org/10.3390/f10080632
Submission received: 13 June 2019 / Revised: 17 July 2019 / Accepted: 25 July 2019 / Published: 27 July 2019
(This article belongs to the Special Issue Forest Genetics and Tree Improvement)

Abstract

:
The normalization of data by choosing suitable reference genes is fundamental for obtaining accurate and reliable results in quantitative real-time polymerase chain reaction (qPCR) analyses. In this study, the expression stability of 12 candidate reference genes of Pinus massoniana under different abiotic stresses was evaluated using four statistical algorithms: geNorm, NormFinder, BestKeeper, and RefFinder. The results indicate that the following genes could be used as reference genes under different treatments: Actin 2 (ACT2) and F-box family gene (F-box) for salinity treatment, cyclophilin (CYP) and alpha-tubulin (TUA) for ABA treatment, actin 7 (ACT7) and CYP for drought treatment, actin 1 (ACT1) and ACT7 for cold treatment, ACT1 and CYP for heat treatment, and TUA and ACT2 for the “Total” group. To validate the suitability of the selected reference genes in this study, the Short-Root protein (SHR), Alpha-pinene synthase (APS), and Pyrabactin resistance-like protein (PYL) gene expression patterns were analyzed. The expression patterns had significant biases when the most unstable reference genes were used for normalization, compared with when the optimum reference gene or gene combinations were used for normalization. These results will be beneficial for further studies on gene transcription in early-stage, unlignified seedlings of P. massoniana.

1. Introduction

Gene expression profiling maintains an important role in molecular biology research, such as aiding in the elucidation of complex regulatory networks of the genetic, signaling, and metabolic pathway mechanisms [1,2]. Quantitative real-time polymerase chain reaction (qPCR) has become one of the most prevalent technologies for monitoring gene expression patterns [3,4]. Absolute and relative quantification are two methods applied to quantify assays of gene expression [5]. Absolute quantification compares the transforming cycle value (Ct) and provides accurate copy numbers of genes with a standard curve [6]. For the relative quantification method, suitable internal control genes are required as references in qPCR [7], and with this method, researchers are more concerned about discrepancies in gene expression analyses rather than the exact copy numbers of genes. Therefore, relative quantification qPCR has been widely used and plays a significant role in gene expression analysis [8].
Despite being superior in reproducibility, sensitivity, and specificity, the relative quantification method requires stable internal controls as a prerequisite condition for gene expression analyses [1,9]. Typically, reference genes are selected as internal controls for normalizing gene expression data in experiments. The expression levels of ideal reference genes are not influenced by experimental conditions, such as tissues, developmental periods, biological stresses, or abiotic stresses. However, increasing evidence shows that ideal reference genes are often not available [1,10,11]. Therefore, it is critical to screen and validate suitable reference genes for the normalization of gene expression data under different experimental conditions [12,13].
Housekeeping genes associated with basic cellular processes are commonly chosen as candidate reference genes in plants. Actin and TUB were the most stable reference genes under abiotic stress and hormone stimuli in carrot [14]. GAPDH and EF-1A were the most stable reference genes under different abiotic sress in Caragana korshinskii [11]. F-box and CYP were suitable reference genes in soybean [15,16]. TUA showed the maximum stability under cold stress in Chinese prickly ash [17]. Besides traditional reference genes, AQP was considered as a new candidate gene in Stevia rebaudiana and olive [18,19]. One of the important evaluation criteria for selecting suitable reference genes under certain experimental conditions is stable expression of the candidate genes under the expected test conditions [1]. Considering the expression stability of the reference genes under different conditions or cross species are also not the same, thus, the reassessment of the expression stability of candidate reference genes in specific experimental conditions is required and essential to ensure the correctness of the relative quantification results in target genes studies [20,21].
Pinus massoniana (Masson pine) is an economically and ecologically important evergreen conifer in southern China and is one of the most significant forest tree species [22,23]. It is widely used in the production of fiber and solid wood [24,25]. Furthermore, more than 70% of resin and turpentine comes from Masson pines in China, which is the world’s largest exporter. Therefore, studying the molecular basis of the economic traits and physiological patterns of Masson pine is of great importance in promoting its breeding process. Illustration of the expression levels of key genes under different conditions is important. Wei et al. selected reference genes for gene expression analysis of P. massoniana under nematode-inoculation conditions [26]. Chen et al. selected reference genes in P. massoniana for different tissues and abiotic stresses [27]. However, to date, there has been no systematic validation of these selected reference genes for P. massoniana in such studies, and there have been no systematic analyses of reference gene screening in early-stage, unlignified seedlings of P. massoniana under different abiotic stresses.
In this study, 12 common reference genes (actin 1 gene (ACT1), actin 2 gene (ACT2), actin 4 gene (ACT4), actin 6 gene (ACT6), actin 7 gene (ACT7), F-box family gene (F-box), elongation factor-1α gene (EF-1A), alpha-tubulin gene (TUA), beta-tubulin gene (TUB), aquaporin protein gene (AQP), cyclophilin gene (CYP), and glyceraldehyde- 3-phosphate dehydrogenase gene (GAPDH)) were selected to assess the stability of gene expression under different abiotic stresses in P. massoniana. Four common statistical algorithms (geNorm [28], NormFinder [29], BestKeeper [30], and RefFinder [31]) were used to calculate the variability of the candidate reference genes and select the most suitable reference genes. The results will lay a foundation for future gene expression pattern research in unlignified seedlings of P. massoniana.

2. Methods

2.1. Plant Materials and Treatments

The seeds of P. massoniana were collected from the National P. massoniana Superior Variety Base of Xinyi, Guangdong, China. The seeds were germinated in pots containing a soil mixture (peat soil:perlite:vermiculite, 3:1:1 v/v) under the conditions of 24 °C and 16 h light/8 h dark. Forty-five-day-old unlignified seedlings in similar growth stages were selected for the following experiments.
For cold and heat treatments, seedlings were exposed to high temperature (42 °C) and low temperature (4 °C) for 0.5 h, 3 h, 6 h, 12 h, 24 h, and 48 h. For drought, salinity, and ABA treatments, washed seedlings were immersed in complete medium containing 10% (w/v) PEG6000, 200 mM NaCl, and 100 μM ABA, respectively, for 0.5 h, 3 h, 6 h, 12 h, 24 h, and 48 h. Drought, salinity, and ABA treatments were performed at room temperature. Untreated seedlings were used as the control. Each treatment was repeated three times independently as biological replicates. All samples were collected and immediately frozen in liquid nitrogen and then stored at −80 ℃ until RNA isolation.
In total, 93 samples were analyzed in the subsequent experiments, including 90 samples with five different stress treatments and three samples as controls.

2.2. RNA Isolation and cDNA Reverse Transcription

Total RNA was extracted using the RNAprep Pure Kit (DP441, Tiangen Biotech, Beijing, China) with a gDNA-removal step according to the manufacturer′s instructions. RNA concentration and purity were measured with a NanoDrop 2000 (Thermo Fisher Scientific, Waltham, Massachusetts, US), and RNA integrity was estimated by 1.2% agarose gel electrophoresis. cDNA (20 μL) was synthesized from 1000 ng of total RNA using the One-step gDNA Removal and cDNA Synthesis Kit (AT311, TransGen Biotech, Beijing, China) with Oligo (dT)18 Primer according to the manufacturer′s instructions.

2.3. Candidate Reference Gene Selection, Primer Design, and Gene Cloning

The sequences of 12 candidate reference genes (ACT1, ACT2, ACT4, ACT6, ACT7, TUA, TUB, EF1A, F-box, AQP, GAPDH, and CYP) were acquired from the GenBank database (https://www.ncbi.nlm.nih.gov/genbank/). Primers were designed by Primer Premier 5.0 software. The primer sequences of candidate reference genes used in this study are summarized in Table 1. The primer specificities and amplicon sizes were verified by 2% agarose gel electrophoresis. A two-fold cDNA dilution series with three replicates per concentration was used to make standard curves for estimation of amplification efficiency (E= (10[−1/slope] −1) × 100%) and the correlation coefficient (R2) [11]. The sequences of 12 candidate reference genes from P. massoniana were cloned using 2 × Rapid Taq Master Mix (P222, Vazyme Biotech, Nanjing, Jiangsu, China) as the polymerase. The PCR (20 μL) mixture contained 10 μL of 2 × Rapid Taq Master Mix, 7 μL of ddH2O, 1 μL of the template cDNA, and 1 μL of each primer (10 μM). The amplification conditions were as follows: 3 min at 95 °C for denaturation; 35 cycles of 20 s at 95 °C (denaturation), 20 s at 58 °C (annealing), and 60 s at 72 °C (extension); and a final step of 5 min at 72 °C for extension.

2.4. Quantitative Real-Time PCR Assay

qRT-PCRs were performed in a StepOne Plus System (Applied Biosystems, Foster City, USA) using SYBR Green Real-time PCR Master Mix (QPK-201, Toyobo Bio-Technology, Shanghai, China). Each PCR mixture (10 μL) contained 1 μL of diluted cDNA (20 × dilution), 5 μL of SYBR Green Real-time PCR Master Mix, 0.4 μL of each primer (10 μM), and 3.2 μL of ddH2O. The amplification conditions were as follows: 95 °C for 60 s for predenaturation, 40 cycles at 95 °C for 10 s for denaturation, and 58 °C for 30 s for annealing and extension. Melting curves were analyzed from 60 °C to 95 °C to confirm primer specificity and lack of primer dimers. Each reaction was repeated three times. The negative controls were included on each plate and for each sample, with ddH2O and total RNA to replace the cDNA.

2.5. Data Analysis

The raw Ct data are shown in Appendix A Table A1. The analysis of the candidate reference gene expression stability was performed by four common software programs: geNorm [28], NormFinder [29], BestKeeper [30], and RefFinder [31]. For analysis with the geNorm and NormFinder algorithms, the raw Ct values were converted into relative quantities. However, for BestKeeper and RefFinder analysis, the Ct data were not transformed.
GeNorm computes the expression stability measure (M-value), analyzes the pairwise variation (V value) in each candidate reference gene, and then excludes the most unstable genes with the highest M-value. Furthermore, pairwise variation (Vn/Vn+1; 0.15 is recommended threshold) determines the optimal number of reference genes for normalization [20,32].
NormFinder computes the expression stability value (SV) on the basis of intra- and intergroups for each reference gene [29]. A lower SV indicates a higher expression stability of this gene.
BestKeeper computes the stability of candidate reference genes based on standard deviation (SD), Pearson correlation coefficient (r), and coefficient of variation (CV) with the Ct data of all candidate genes. The lowest SD and CV values represent the most stable gene. The range of variation in SD should not surpass 1 [20,32].
RefFinder assigns an appropriate weight to an individual candidate gene and calculates the geometric mean of their weights for the overall final ranking, based on the rankings from the currently available major programs (geNorm, Normfinder, BestKeeper, and the comparative ΔCt method) [31].

2.6. Validation of Reference Genes

To validate the reliability of selected reference genes, the two most stable and two unstable reference genes and the optimum internal reference gene combinations were used to normalize the relative expression patterns of SHR, APS, and PYL in each experimental condition, according to the rankings from RefFinder and geNorm, respectively. The relative expression levels were calculated by the 2−ΔΔCt method [8].

3. Results

3.1. Selection of Reference Genes, Amplification Specificity, and Efficiency

Based on the GenBank database, 12 candidate reference genes (ACT1, ACT2, ACT4, ACT6, ACT7, TUA, TUB, EF1A, F-box, AQP, GAPDH, and CYP) were cloned from P. massoniana. The primer pairs of all candidate reference genes and target genes were designed for qPCR, and the amplicon lengths were between 95 and 114 bp. Agarose gel electrophoresis (Appendix A Figure A1) and melting curve analysis (Appendix A Figure A2) were used to determine primer specificity. The amplification efficiency of qPCR across all 12 reference genes and three validated genes varied from 91.0% to 104.8%, with R2 varying from 0.989 to 0.999 (Table 1 and Appendix A Figure A3).

3.2. Expression Profiles of the Reference Genes

The expression levels of candidate reference genes were quantified by Ct values; the high expression stability of the genes was reflected by the low SV. The raw Ct values for all samples in this study are listed in Appendix A Table A1 (there were no Ct values in the negative controls), and a box and whiskers plot was used to describe the raw Ct value distribution (Figure 1 and Appendix A Table A1). The 12 candidate reference genes had a wide range of expression across all samples in this study (20.89 ≤ Ct ≤ 33.96). The results indicate that there were five genes (ACT1, CYP, TUA, GAPDH, and AQP) with average Ct values between 23 and 25 cycles that displayed high expression levels; the other seven genes (EF1A, ACT4, TUB, F-box, ACT7, ACT6, and ACT2), which presented Ct values between 25 and 31, displayed intermediate expression levels. ACT1 was the most highly expressed gene among the 12 candidate genes (mean Ct of 23.32). On the other hand, ACT2 was the least highly expressed gene (mean Ct of 30.86). In addition, every reference gene had different coefficients of variation (lower values represent less variability) in this study, as shown in Figure 1; TUB had the most variation, whereas F-box had the least variation, which indicated that TUB was the most unstable gene and F-box was the most stable gene of all the reference genes.

3.3. Expression Stability Analysis of the Reference Genes

In this study, five experimental conditions were performed, including cold, heat, drought, ABA, and salinity treatments. Furthermore, these experimental conditions were sorted into the “Total” group (all experimental conditions). To select appropriate reference genes for these experimental conditions, four software programs (geNorm [28], NormFinder [29], BestKeeper [30], and RefFinder [31]) were used to evaluate the stability of the 12 candidate reference genes.

3.4. geNorm Analysis

The geNorm analysis results are presented in Table 2. Under different experimental conditions, the reference genes with the highest stability were different. The twelve reference genes in the five experimental conditions and “Total” group showed high expression stabilities with threshold values below 1.5. For cold treatment, ACT1 and TUA (M = 0.157) were the most stable genes, while EF-1A (M = 0.634) was the most variable gene. For heat treatment, CYP and GAPDH (M = 0.309) were identified as the most stable genes, whereas the ACT6 (M = 0.47) gene was the least stable. For salinity treatment, ACT2 and F-box (M = 0.275) were the most stable genes, while AQP (M = 0.787) was the most variable gene. For drought treatment, ACT4 and TUB (M = 0.280) were identified as the most stable genes, whereas ACT6 (M = 0.785) was the least stable. For ABA treatment, CYP and TUA (M = 0.246) were the most stable genes, while TUB (M = 0.663) was the most variable gene. A combination of individual samples in the “Total” group showed that ACT1 and TUA (M = 0.437) featured the most stable expression, whereas the most variable genes included AQP (M = 0.882) and ACT6 (M = 0.932).
The optimal number of reference genes for normalization depends on pairwise variation (Vn/n+1). When Vn/n+1 <0.15, it suggests that an additional reference gene is not necessary for normalization. For four experimental conditions (cold, salinity, drought, and ABA), two reference genes were sufficient for accurate normalization (Figure 2); the most stable gene pairs for cold, salinity, drought, and ABA conditions were ACT1 and TUA, ACT2 and F-box, ACT4 and TUB, and CYP and TUA, respectively (Appendix A Figure A4). For heat treatment and the “Total” group, V2/3 > 0.15, and V3/4 < 0.15 (Figure 2). Therefore, the genes CYP, GAPDH, and ACT2, and ACT1, TUA, and ACT7 were chosen, respectively (Appendix A Figure A4).

3.5. NormFinder Analysis

The SV of each candidate reference gene was also analyzed by NormFinder, wherein a lower SV indicates a higher expression stability. In this study, the results of NormFinder analysis were similar to those of geNorm analysis (Table 2). For cold treatment, GAPDH was the most stable reference gene, and EF1A was the most variable gene. For heat treatment, ACT1 was the most stable reference gene, and ACT6 was the most variable gene. For salinity treatment, F-box was the most stable reference gene, and AQP was the most variable gene. For drought treatment, ACT7 was the most stable reference gene, and ACT6 was the most variable gene. For ABA treatment, CYP was the most stable reference gene, and TUB was the most variable gene. For the “Total” group, ACT2 was the most stable reference gene, and ACT6 was the most variable gene.

3.6. BestKeeper Analysis

The results of the BestKeeper analysis are also shown in Table 2. The results indicated that all candidate reference genes were markedly stable when expressed under three experimental conditions (cold, drought, and ABA) in this study. The rankings by BestKeeper analysis showed that the most stable reference genes were ACT1 (CV ± SD = 0.99 ± 0.24) and TUA (CV ± SD = 1.20 ± 0.29) for cold treatment. F-box (CV ± SD = 0.84 ± 0.24) and ACT7 (CV ± SD = 1.06 ± 0.31) were the most stable genes for drought treatment. GAPDH (CV ± SD = 1.78 ± 0.42) and F-box (CV ± SD = 1.51 ± 0.42) were the most stable genes for ABA treatment. For salinity treatment and the “Total” group, 10 reference genes were expressed stably; the most stable genes were ACT2 (CV ± SD = 1.13 ± 0.34) and EF1A (CV ± SD = 1.24 ± 0.30) for salinity treatment. F-box (CV ± SD = 1.57 ± 0.44) and CYP (CV ± SD = 2.53 ± 0.59) showed the most stable expression under cold treatment. For heat treatment, only seven reference genes were expressed stably; the most stable genes were ACT4 (CV ± SD = 1.92 ± 0.55) and TUA (CV ± SD = 2.06 ± 0.50) (Table 2).

3.7. RefFinder Analysis

The geometric mean rankings were estimated from geNorm, NormFinder, and BestKeeper programs by RefFinder algorithms [31]. This allowed us to generate a recommended comprehensive ranking of reference genes for accurate transcript normalization in each experimental set. The rankings by RefFinder analysis showed that ACT2 and F-box were the most stable genes for salinity treatment, that CYP and TUA were the most stable genes for ABA treatment, that ACT7 and CYP were the most stable genes for drought treatment, that ACT1 and ACT7 were the most stable genes for cold treatment, that ACT1 and CYP were the most stable genes for heat treatment, and that TUA and ACT2 were the most stable genes for the “Total” group (Table 2).

3.8. Reference Gene Validation

To validate the accuracy of selected reference genes, the relative expression levels of SHR, APS, and PYL were analyzed in all the experimental conditions involved in this study. For each experimental condition, the two most stable and two unstable reference genes according to the combined analytic results obtained by geNorm, NormFinder, and BestKeeper, and the reference gene combination according to geNorm, were selected for normalization.
Under cold treatment, SHR was upregulated, and the highest expression levels were observed at 3 h (3.91-fold change) (Figure 3A). APS was downregulated at 3 h and 48 h and upregulated at 0.5 h, 6 h, 12 h, and 24 h, especially at 6 h (2.28-fold change) (Figure 3B). There was no significant difference in the expression of PYL, and PYL had the highest expression at 3 h (1.56-fold change) (Figure 3C).
Under drought treatment, SHR was downregulated at 6 h and upregulated at 0.5 h, 3 h, 12 h, 24 h, and 48 h, especially at 24 h (3.34-fold change) (Figure 3D). APS was upregulated, and the highest expression levels were observed at 24 h (21.02-fold change) (Figure 3E). PYL was downregulated at 0.5 h, 3 h, 6 h, and 48 h and upregulated at 12 h and 24 h, and the highest expression levels were observed at 24 h (3.47-fold change) (Figure 3F).
Under ABA treatment, SHR was downregulated at 6 h and 12 h and upregulated at 0.5 h, 3 h, 24 h, and 48 h, especially at 48 h (2.68-fold change) (Figure 3G). APS was upregulated, and the highest expression levels were observed at 3 h (7.56-fold change) (Figure 3H). PYL was downregulated at 0.5 h, 3 h, 6 h, and 12 h and upregulated at 24 h and 48 h, and the highest expression levels were observed at 48 h (2.45-fold change) (Figure 3I).
Under salinity treatment, SHR was upregulated, and the highest expression levels were observed at 3 h (4.26-fold change) (Figure 3J). APS was upregulated, and the highest expression levels were observed at 3 h (15.58-fold change) (Figure 3K). PYL was downregulated at 6 h and upregulated at 0.5 h, 3 h, 12 h, 24 h, and 48h, and the highest expression levels were observed at 24 h (3.21-fold change) (Figure 3L).
Under heat treatment, SHR was upregulated, and the highest expression levels were observed at 48 h (6.03-fold change) (Figure 3M). APS was downregulated at 6 h and upregulated at 0.5 h, 3h, 12 h, 24 h, and 48 h, especially at 0.5 h (4.06-fold change) (Figure 3N). PYL was upregulated, and the highest expression levels were observed at 24 h (3.11-fold change) (Figure 3O).
Our results confirm that using different reference genes for normalization causes great differences among the expression patterns of SHR, APS, and PYL. When the stable reference genes and optimum reference gene combinations were used for normalization, the expression patterns of SHR, APS, and PYL were similar. However, the expression patterns of SHR, APS, and PYL had significant biases when the most unstable reference genes were used for normalization compared with when the optimum reference gene combinations were used for normalization. The results illustrate that a stably-expressed reference gene was essential for the accurate normalization of target gene expression.

4. Discussion

Level and pattern analyses of gene expression is essential for the functional analysis of genes in different organisms and conditions [33]. Although gene expression levels can be obtained using different technologies, qPCR is currently the prevailing tool for the analysis of gene expression patterns in molecular biology research because of its speed, reproducibility, sensitivity, and specificity [34,35,36]. The use of reliable reference genes with steady expression levels and proper expression abundance is necessary for target gene expression pattern evaluation in different species or under different conditions by correcting the errors of RNA quantity and reverse transcription efficiency [2,37]. Ideal reference genes should be expressed at a steady and continuous level across different experimental conditions. An increasing number of studies have indicated that it is important to select an appropriate reference gene or gene combination using a systematic experiment due to the lack of available ideal reference genes [32]. Many studies on reference gene selection have been reported in plant organisms under different cultivars, developmental stages, biotic stresses, and abiotic stresses; for instance, rice [38], wheat [39], peach [40], poplar [41], Prunus mume [42], Coffea Arabica [43], Metasequoia [44], P. pinaster [45], and Picea abies [45], among others. However, to the best of our knowledge, the selection of reference genes has been performed in only nematode-inoculated and lignified materials of P. massoniana.
In this study, 12 commonly used reference genes (ACT1, ACT2, ACT4, ACT6, ACT7, TUA, TUB, EF1A, F-box, AQP, GAPDH, and CYP) were selected as candidate reference genes for analysis under five abiotic stress conditions for the lignified seedlings of P. massoniana (45-days-old). The twelve candidate reference genes used in this study achieved an appropriate expression abundance (20.89 ≤ Ct ≤ 33.96), which allowed further assessment of their expression stability. To date, this study is the first report of a systematic analysis of reference genes in different abiotic stress conditions in lignified plant materials of P. massoniana. Many researchers have chosen several statistical methods simultaneously to assess the optimal reference genes under different experimental conditions and avoided a one-sided analysis of the expression stability of reference genes that may occur when using only one algorithm [46,47]. In the present study, four of the most popular statistical programs (geNorm, NormFinder, BestKeeper, and RefFinder) were used to screen and determine appropriate reference genes. The four statistical programs generated inconsistent stability rankings in each experimental condition, since these algorithms had different sensitivities to coregulation genes [28,29,30,31]. Many published studies indicated that the stability ranking results from different statistical algorithms were approximately the same, and the most divergent ranking of gene stability was achieved with BestKeeper [32]. In this study, for the ABA treatment, CYP and TUA were identified as the most stable genes by geNorm, NormFinder, and RefFinder, while BestKeeper indicated that F-box and GAPDH were the best reference genes, although F-box and GAPDH were ranked as the 11th and seventh genes by both geNorm and NormFinder.
There is a possibility to obtain declinational gene expression analysis for normalization by a single reference gene [28,48]. Therefore, the experimental error could be reduced using two or more internal control genes as references [49]. In this study, the optimal number of reference genes for gene expression analysis under different abiotic stress conditions was determined using the geNorm program. Our results showed that under cold, drought, salinity, and ABA treatments, the pairwise variation was V2/3 < 0.15, which indicated that two reference genes were sufficient for optimal normalization. However, in the heat treatment and “Total” groups, V2/3 > 0.15, which indicates that more genes were needed. However, although the use of two or more reference genes can make gene expression analyses more exact, this practice is not an essential standard [28].
The suitability of the selected reference genes was assessed by analyzing the expression levels in three target genes related to a GRAS family transcription factor (SHR), alpha-pinene biosynthesis (APS), and an ABA receptor (PYL). The SHR protein, as a member of the GRAS transcription factor family, is related to the process of Masson pine root development; APS is a key enzyme-coding gene related to alpha-pinene biosynthesis; and PYL is an ABA receptor in the ABA signaling pathway. Their expression levels may be directly related to the process of growth and development, resin biosynthesis, and response to drought stress. In our study, the expression levels of SHR, APS, and PYL were different under different abiotic stresses.
Furthermore, we compared the most stable, unstable, and optimal reference gene combinations in transcriptional expression analysis of SHR, APS, and PYL, and the results were quite different. The results indicate that when the most stable reference genes and gene combinations were used for normalization, the transcriptional expression patterns of SHR, APS, and PYL were similar. However, the expression patterns of SHR, APS, and PYL had significant biases when the least stable reference genes were used for normalization compared with when the most stable reference genes and gene combinations were used for normalization. These results indicate that the reference genes screened in this study are reliable.
Significantly, the expression stability of these candidate reference genes had not been evaluated for various tissues from other developmental stages of P. massoniana or other plant species under different abiotic stresses in our study. Some reported findings suggest that the stability of reference genes exhibited organ-specificity as well as species-specificity in plants, generally [50]. Tu et al. reported that the suitable reference genes for flower organs and vegetative organs were not consistent [51]. A similar result had been found in Petunia hybrid [52]. Yang et al. reported GAPDH was the most stable reference gene under heat stress in Caragana korshinskii [11]. By contrast, GAPDH had a low expression stability in pisum sativum [53]. Similarly, GAPDH was an unsuitable reference gene under different abiotic stresses in P. massonina in our study.
The selected stable reference genes in this study will be beneficial for more accurate quantification of gene expression levels in early-stage, unlignified seedlings of P. massoniana under different abiotic stresses.

5. Conclusions

In this study, we selected specific reference genes for qPCR analysis of unlignified seedlings in P. massoniana different abiotic stresses. The most suitable combination for qPCR analysis across all samples was TUA and ACT1. The following genes could be used as reference genes under different abiotic stresses: ACT2 and F-box for salinity treatment, CYP and TUA for ABA treatment, ACT7 and CYP for drought treatment, ACT1 and ACT7 for cold treatment, and ACT1 and CYP for heat treatment. In conclusion, this study provides an important framework for quantifying reference gene expression levels in P. massoniana.

Author Contributions

Experimental design: K.J. and P.Z.; material collection and performing the experiments: P.Z. and Y.M.; data analysis: P.Z.; writing—original draft: P.Z. and L.Z.; writing—review and editing: K.J., L.Z., Y.C., and R.L.; all authors approved the final draft.

Funding

This study was supported by the National Key R&D Program of China (2017YFD0600304) and a project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).

Acknowledgments

We thank Jing Ye for the suggestions about data analysis.

Conflicts of Interest

The authors declare that they have no competing interests.

Appendix A

Figure A1. Amplification of q-PCR specific primers. Amplified fragments of 12 candidate reference genes and three validate genes in Pinus massonina were separated by 2.0% agarose gel electrophoresis. Lane 1 and 17 (GsDL2502, Shanghai Generay Biotech, Shanghai, China); Lane 2: ACT1(100bp); Lane 3: ACT2(97bp); Lane 4: ACT4(106bp); Lane 5: ACT6(99bp); Lane 6: ACT7(114bp); Lane 7: AQP(103bp); Lane 8: CYP(97bp); Lane 9: EF1A(97bp); Lane 10: F-box(109bp); Lane 11: GAPDH(96bp); Lane 12: TUA(95bp); Lane 13: TUB(108bp); Lane 14: SHR(97bp); Lane 15: APS(105bp); Lane 16: PYL(110bp).
Figure A1. Amplification of q-PCR specific primers. Amplified fragments of 12 candidate reference genes and three validate genes in Pinus massonina were separated by 2.0% agarose gel electrophoresis. Lane 1 and 17 (GsDL2502, Shanghai Generay Biotech, Shanghai, China); Lane 2: ACT1(100bp); Lane 3: ACT2(97bp); Lane 4: ACT4(106bp); Lane 5: ACT6(99bp); Lane 6: ACT7(114bp); Lane 7: AQP(103bp); Lane 8: CYP(97bp); Lane 9: EF1A(97bp); Lane 10: F-box(109bp); Lane 11: GAPDH(96bp); Lane 12: TUA(95bp); Lane 13: TUB(108bp); Lane 14: SHR(97bp); Lane 15: APS(105bp); Lane 16: PYL(110bp).
Forests 10 00632 g0a1
Figure A2. Melting curve for 12 candidate reference genes and three validate genes in P. massonina. The melting curve for ACT1, ACT2, ACT4, ACT6, ACT7, AQP, CYP, EF1A, F-box, GAPDH, TUA, TUB, SHR, APS, and PYL, with single peak obtained from three technical replicates of different cDNA with no template control.
Figure A2. Melting curve for 12 candidate reference genes and three validate genes in P. massonina. The melting curve for ACT1, ACT2, ACT4, ACT6, ACT7, AQP, CYP, EF1A, F-box, GAPDH, TUA, TUB, SHR, APS, and PYL, with single peak obtained from three technical replicates of different cDNA with no template control.
Forests 10 00632 g0a2
Figure A3. Standard curves of 12 candidate reference genes and three validate genes in P. massonina. Standard curves of ACT1, ACT2, ACT4, ACT6, ACT7, AQP, CYP, EF1A, F-box, GAPDH, TUA, TUB, SHR, APS, and PYL. The linear correlation (R2) and PCR efficiencies (E = (10[-1/slope] -1) × 100%) were calculated from the standard curve.
Figure A3. Standard curves of 12 candidate reference genes and three validate genes in P. massonina. Standard curves of ACT1, ACT2, ACT4, ACT6, ACT7, AQP, CYP, EF1A, F-box, GAPDH, TUA, TUB, SHR, APS, and PYL. The linear correlation (R2) and PCR efficiencies (E = (10[-1/slope] -1) × 100%) were calculated from the standard curve.
Forests 10 00632 g0a3
Figure A4. The average expression stability values (M value) of the candidate reference genes analyzed by geNorm. Expression stability was evaluated in samples from P. massonina submitted to cold stress, ABA stress, heat stress, drought stress, salinity stress, and total group.
Figure A4. The average expression stability values (M value) of the candidate reference genes analyzed by geNorm. Expression stability was evaluated in samples from P. massonina submitted to cold stress, ABA stress, heat stress, drought stress, salinity stress, and total group.
Forests 10 00632 g0a4
Table A1. Raw Ct values of 12 reference genes in three validate genes in P. massonina.
Table A1. Raw Ct values of 12 reference genes in three validate genes in P. massonina.
ACT-1ACT-2ACT-4ACT-6ACT-7AQPCYPEF1AF-boxGAPDHTUATUBSHRAPSPYL
Control0h23.93 30.57 28.69 31.83 29.96 24.71 23.93 22.77 28.81 23.84 23.87 27.82 30.40 27.88 27.86
23.83 30.49 28.79 31.63 29.76 24.92 23.81 22.90 28.40 24.08 23.99 27.93 30.27 27.87 27.62
23.78 30.86 28.54 31.72 29.80 22.84 23.79 24.90 28.38 23.87 23.80 28.52 30.17 27.72 27.54
ABA0.5h22.55 29.77 27.96 29.97 27.74 23.50 22.84 23.86 28.85 23.93 23.41 25.76 29.62 24.55 27.83
22.88 30.61 28.16 30.29 28.19 23.80 23.31 24.40 28.79 24.34 23.54 26.84 28.84 25.86 27.48
22.94 30.39 27.49 30.46 28.38 23.93 23.44 24.41 28.31 23.63 22.96 26.53 28.79 24.74 27.02
3h23.22 31.19 26.75 30.49 28.27 24.27 23.23 24.57 28.07 23.29 22.89 26.56 28.90 24.36 27.90
22.90 28.40 25.85 30.81 28.63 23.01 22.76 24.29 27.94 22.91 22.86 26.42 28.36 23.70 27.58
23.18 28.87 26.10 30.68 28.86 23.29 22.98 24.51 28.01 23.07 22.88 26.52 28.55 23.90 27.72
6h22.76 28.47 26.27 30.46 28.41 23.25 22.66 23.67 28.16 22.40 21.94 26.37 28.91 23.65 26.95
22.84 28.89 26.83 30.34 28.83 23.12 22.58 24.14 28.66 23.69 22.64 26.83 29.30 24.07 27.61
21.84 29.26 26.66 29.39 27.08 22.30 21.80 22.94 28.24 22.92 22.29 25.78 29.02 23.88 27.26
12h21.54 29.02 26.57 28.84 26.93 21.89 21.64 22.72 27.86 22.86 21.91 25.85 28.68 23.81 26.92
21.89 29.04 26.01 29.04 27.09 22.07 21.91 23.07 27.52 22.43 21.34 25.83 27.90 22.80 26.40
21.70 28.92 26.86 28.95 27.01 21.96 21.82 22.92 27.14 22.49 21.54 25.73 28.25 23.06 26.72
24h22.85 29.80 27.66 30.75 27.96 22.68 22.92 24.05 27.85 23.19 23.44 26.83 29.01 24.89 26.87
22.65 29.62 27.85 30.61 27.88 22.36 22.41 23.97 27.76 23.52 23.62 26.57 28.97 25.42 27.00
23.60 30.47 27.18 31.50 28.91 23.48 23.55 25.05 27.22 22.78 22.96 27.38 28.33 24.75 26.39
48h23.00 31.88 27.99 31.76 29.23 23.94 23.76 24.95 27.85 23.66 24.51 29.42 29.12 28.27 26.80
23.50 31.65 27.44 31.98 29.79 24.19 23.90 25.65 27.04 23.41 23.84 29.82 28.83 26.09 26.20
23.05 31.40 27.78 31.53 29.18 23.70 23.54 25.10 27.53 23.55 24.21 29.39 29.04 26.90 26.56
Cold0.5h22.96 29.99 28.37 31.67 30.30 22.90 23.45 24.76 27.88 22.89 24.56 27.95 29.73 28.80 28.15
23.59 30.80 28.68 30.38 29.01 23.55 22.34 25.20 28.52 23.80 23.17 27.93 28.55 26.77 26.79
23.74 30.93 28.95 30.73 29.49 23.65 22.64 25.41 28.61 23.89 23.66 28.23 28.96 27.42 27.26
3h24.89 32.25 29.87 31.72 30.25 23.70 22.86 25.71 28.71 24.88 24.95 29.87 29.01 29.99 27.58
23.86 30.77 28.68 31.85 30.40 22.86 23.02 24.92 27.62 23.65 24.82 28.89 28.89 29.35 27.86
24.46 31.40 29.25 31.68 30.34 23.45 22.88 25.44 28.16 24.33 24.73 29.48 28.94 29.40 27.68
6h23.89 30.88 28.30 31.79 29.53 23.30 23.82 25.83 28.91 23.88 23.92 26.63 29.86 25.91 28.37
23.36 30.64 27.69 31.75 29.94 22.77 23.46 25.40 28.39 23.55 23.81 27.83 29.47 27.88 27.85
23.55 30.48 27.87 31.70 29.45 22.92 23.34 25.49 28.48 23.62 23.71 27.02 29.36 26.21 27.66
12h23.90 31.01 28.85 31.91 30.57 23.88 23.91 26.30 28.08 24.38 24.44 27.93 29.91 27.87 27.77
23.69 30.51 28.40 31.32 30.00 23.82 23.25 25.90 27.61 23.76 23.84 27.40 28.88 26.40 27.51
23.68 31.63 28.65 31.40 29.84 23.66 23.15 26.01 27.74 23.90 23.77 27.34 28.97 26.67 27.34
24h24.85 32.94 29.88 32.85 30.79 25.30 24.13 27.89 29.89 25.47 24.94 29.83 29.79 28.87 28.94
23.65 31.64 28.55 31.67 30.47 23.88 23.17 27.07 28.75 24.87 23.76 28.86 28.84 26.22 27.65
23.75 31.36 28.66 32.10 30.51 23.90 23.53 26.99 28.87 24.87 24.13 28.85 29.31 27.06 28.04
48h23.94 31.64 29.35 30.52 29.73 23.89 22.75 26.88 28.82 23.92 23.89 28.99 28.86 29.82 27.80
23.59 31.08 28.66 30.64 29.17 23.82 22.63 26.39 28.42 23.45 23.67 28.72 28.39 29.77 27.72
23.65 31.71 28.93 30.62 29.44 23.76 22.54 26.65 28.53 23.65 23.61 28.77 28.61 29.62 27.69
Salinity0.5h22.62 30.31 26.61 30.82 27.58 22.89 22.66 23.73 27.86 23.66 22.84 25.37 28.65 24.67 26.60
22.82 30.54 27.84 30.90 28.61 23.23 23.01 24.01 28.21 24.53 23.15 27.10 28.82 25.39 26.47
22.60 30.08 26.92 30.58 27.94 22.70 22.54 23.50 27.75 23.92 22.88 25.94 28.65 24.90 26.33
3h21.73 30.04 26.68 29.81 28.40 22.89 22.89 24.36 27.81 22.02 21.62 26.54 27.76 23.76 26.35
20.91 29.22 26.33 29.09 27.87 22.32 22.46 23.88 27.28 21.76 20.89 25.91 26.93 22.40 25.78
21.08 29.93 26.81 29.39 28.52 22.87 22.94 24.46 27.67 21.71 21.12 26.36 27.14 22.82 25.87
6h22.51 29.79 27.02 30.63 27.85 22.77 22.92 24.01 27.93 23.51 22.59 26.81 28.73 23.94 27.10
22.89 29.78 26.90 30.75 27.53 22.85 22.73 23.94 27.71 23.88 22.90 26.90 29.81 24.54 27.55
22.57 30.50 27.49 30.95 27.99 23.26 23.43 24.56 28.29 23.50 22.62 27.15 28.97 24.08 27.21
12h22.85 30.57 27.90 29.87 28.45 23.92 23.39 24.80 27.89 23.28 22.91 27.38 28.98 23.95 26.91
22.57 30.04 27.29 29.78 27.97 23.44 22.90 23.95 27.28 22.36 22.75 26.99 28.31 23.76 26.96
22.77 30.24 27.57 29.69 28.30 23.56 23.09 24.30 27.48 22.47 22.75 27.17 28.55 23.85 27.00
24h23.42 31.12 28.91 30.89 28.92 25.98 23.85 24.88 29.17 23.75 23.80 29.42 29.67 25.79 26.80
22.55 30.80 28.49 29.88 28.70 26.04 24.19 24.46 28.83 22.47 22.92 28.93 28.59 23.85 26.11
22.92 30.95 28.81 30.26 28.90 25.88 24.05 24.63 28.96 22.84 23.22 29.00 28.98 24.45 26.30
48h23.45 30.70 28.69 31.45 28.77 25.87 22.95 23.95 28.79 22.84 24.55 28.85 29.81 26.96 27.71
23.83 30.67 28.57 31.69 28.87 25.80 23.19 23.83 28.69 23.20 24.72 29.07 29.46 27.89 27.89
23.85 30.90 28.87 31.35 29.06 26.16 23.50 24.15 28.84 23.06 24.80 29.12 29.63 27.62 27.80
Drought0.5h22.94 31.58 27.86 30.91 29.28 25.66 23.79 25.46 29.42 23.91 22.82 27.55 29.91 24.95 27.30
22.19 30.87 27.34 27.79 28.53 24.91 23.59 24.95 28.90 23.59 22.63 26.82 28.84 23.98 27.70
22.18 31.66 27.92 27.78 28.94 25.19 23.92 25.40 29.39 23.47 22.59 27.28 28.95 24.14 27.37
3h23.03 30.84 27.78 31.05 28.90 25.21 22.93 24.86 28.88 23.85 23.14 26.87 28.93 24.88 27.84
21.85 28.97 26.08 30.48 28.02 24.40 22.40 23.89 28.17 22.51 21.87 26.00 28.83 23.38 26.88
22.23 29.87 26.81 30.71 28.71 24.89 22.92 24.63 28.79 22.90 22.29 26.59 28.87 23.83 27.19
6h23.37 30.69 27.82 30.91 29.12 25.54 23.90 25.47 29.04 24.09 22.96 26.70 29.80 24.17 28.25
23.75 30.58 28.00 31.54 29.22 25.38 23.86 24.88 28.95 24.27 23.77 27.90 30.19 24.84 27.97
23.33 30.36 27.46 30.92 28.92 24.93 23.46 24.68 28.54 23.99 23.19 26.81 29.86 24.37 27.97
12h22.49 30.46 27.89 29.82 28.73 24.78 22.93 24.74 28.65 22.88 22.79 26.88 28.87 22.56 26.65
22.88 31.28 28.47 30.16 29.35 25.28 23.56 24.91 28.76 23.10 23.46 28.32 28.68 24.94 26.82
22.46 31.22 28.50 29.70 29.25 25.41 23.61 24.99 28.88 22.64 22.88 27.91 28.59 23.29 26.45
24h22.30 31.22 28.81 30.24 29.07 25.66 23.92 25.56 28.66 23.05 22.91 27.89 28.85 23.64 26.58
21.89 31.81 28.46 30.31 29.31 25.56 23.97 25.80 28.85 22.70 23.33 28.49 28.57 23.70 25.79
21.34 31.52 28.78 29.69 29.36 25.90 24.20 25.95 29.00 22.09 22.86 28.59 28.17 22.96 25.38
48h24.54 31.72 28.89 31.93 29.71 24.96 23.67 25.49 28.34 24.75 25.10 29.41 29.92 27.72 27.88
23.86 30.47 27.88 31.21 29.24 24.31 22.96 24.86 27.93 24.03 23.76 27.88 29.61 25.58 27.75
23.92 31.43 28.34 31.39 29.68 24.81 23.46 25.22 28.15 23.97 23.97 28.38 29.41 25.94 27.38
Heat0.5h23.95 31.52 27.82 29.53 29.46 25.92 23.92 25.90 29.24 22.94 22.74 28.40 28.47 25.47 26.67
23.43 30.66 27.67 29.49 28.96 24.74 23.23 25.01 28.52 22.92 22.87 27.84 29.00 25.36 26.52
23.58 30.75 27.98 29.94 28.84 24.95 23.39 25.25 28.72 23.08 23.04 27.93 28.93 25.53 26.86
3h24.55 30.89 28.19 31.92 29.94 25.51 23.47 25.87 28.72 23.87 24.69 29.72 29.86 29.25 26.53
24.61 31.15 27.52 31.15 29.40 26.20 23.54 25.63 28.64 22.71 23.90 29.76 28.78 26.65 25.98
24.90 31.48 27.63 31.11 29.57 26.34 23.88 25.96 28.95 23.03 24.05 30.39 29.08 27.35 26.15
6h23.78 30.59 28.78 32.62 29.62 25.70 23.94 25.78 27.66 24.38 24.84 28.63 29.88 28.78 27.60
23.88 31.50 28.43 32.14 29.79 25.87 23.78 25.89 28.13 23.88 23.85 29.08 29.40 29.26 27.00
23.90 31.06 28.65 32.58 29.87 25.67 23.83 25.90 27.96 23.98 24.20 29.08 29.57 28.95 27.12
12h24.32 31.95 27.65 28.33 31.81 26.56 24.02 27.92 28.38 23.87 23.88 29.76 28.59 26.92 26.72
23.90 31.63 27.87 27.74 30.89 25.95 23.64 27.36 27.86 23.80 23.82 29.08 28.95 26.63 27.17
24.30 31.96 27.61 27.86 31.48 26.52 24.02 27.89 28.08 23.74 23.74 29.60 28.73 26.85 26.92
24h25.91 33.44 29.89 32.81 32.93 28.37 26.43 28.82 29.31 26.75 25.61 31.56 30.22 28.16 28.44
25.34 33.59 28.73 31.96 31.83 27.90 25.84 27.99 28.81 25.79 24.78 30.61 29.98 30.09 27.89
25.20 32.89 28.96 32.50 31.97 27.58 25.85 27.95 28.67 25.86 24.87 30.90 30.00 28.86 27.94
48h24.67 33.64 28.79 31.05 32.03 26.87 25.94 27.86 28.18 25.86 24.31 30.01 29.93 29.32 27.83
25.60 33.96 29.40 31.20 32.62 27.86 26.66 28.23 28.49 25.73 24.88 30.61 29.52 30.24 28.46
24.90 33.72 29.01 31.12 32.01 27.07 26.04 27.63 28.01 25.72 24.57 30.11 29.80 29.84 28.05

References

  1. Zhou, X.; Liu, J.; Zhuang, Y. Selection of appropriate reference genes in eggplant for quantitative gene expression studies under different experimental conditions. Sci. Hortic. 2014, 176, 200–207. [Google Scholar] [CrossRef]
  2. Wan, H.; Zhao, Z.; Qian, C.; Sui, Y.; Malik, A.A.; Chen, J. Selection of appropriate reference genes for gene expression studies by quantitative real-time polymerase chain reaction in cucumber. Anal. Biochem. 2010, 399, 257–261. [Google Scholar] [CrossRef]
  3. Derveaux, S.; Vandesompele, J.; Hellemans, J. How to do successful gene expression analysis using real-time PCR. Methods 2010, 50, 227–230. [Google Scholar] [CrossRef] [PubMed]
  4. Ginzinger, D.G. Gene quantification using real-time quantitative PCR. Exp. Hematol. 2002, 30, 503–512. [Google Scholar] [CrossRef]
  5. Ye, J.; Jin, C.F.; Li, N.; Liu, M.H.; Fei, Z.X.; Dong, L.Z.; Li, L.; Li, Z.Q. Selection of suitable reference genes for qRT-PCR normalisation under different experimental conditions in Eucommia ulmoides Oliv. Sci. Rep. 2018, 8, 15043. [Google Scholar] [CrossRef] [PubMed]
  6. Leong, D.T.; Gupta, A.; Bai, H.F.; Wan, G.; Yoong, L.F.; Too, H.P.; Chew, F.T.; Hutmacher, D.W. Absolute quantification of gene expression in biomaterials research using real-time PCR. Biomaterials 2007, 28, 203–210. [Google Scholar] [CrossRef] [PubMed]
  7. Pfaffl, M.W. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001, 29, e45. [Google Scholar] [CrossRef] [PubMed]
  8. Livak, K.J.; Schmittgen, T.D. Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef] [PubMed]
  9. Bustin, S.A.; Benes, V.; Garson, J.A.; Hellemans, J.; Huggett, J.; Kubista, M.; Mueller, R.; Nolan, T.; Pfaffl, M.W.; Shipley, G.L.; et al. The MIQE guidelines: Minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 2009, 55, 611–622. [Google Scholar] [CrossRef] [PubMed]
  10. Chen, J.; Huang, Z.; Huang, H.; Wei, S.; Liu, Y.; Jiang, C.; Zhang, J.; Zhang, C. Selection of relatively exact reference genes for gene expression studies in goosegrass (Eleusine indica) under herbicide stress. Sci. Rep. 2017, 7, 46494. [Google Scholar] [CrossRef]
  11. Yang, Q.; Yin, J.; Qi, L.; Yang, F.; Wang, R.; Li, G. Reference gene selection for qRT-PCR in Caragana korshinskii Kom. under different stress conditions. Mol. Biol. Rep. 2014, 41, 2325–2334. [Google Scholar] [CrossRef] [PubMed]
  12. Dheda, K.; Huggett, J.F.; Bustin, S.A.; Johnson, M.A.; Rook, G.; Zumla, A. Validation of housekeeping genes for normalizing RNA expression in real-time PCR. Biotechniques 2004, 37, 112–119. [Google Scholar] [CrossRef] [PubMed]
  13. Huggett, J.; Dheda, K.; Bustin, S.; Zumla, A. Real-time RT-PCR normalisation; strategies and considerations. Genes Immun. 2005, 6, 279–284. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Tian, C.; Jiang, Q.; Wang, F.; Wang, G.L.; Xu, Z.S.; Xiong, A.S. Selection of suitable reference genes for qPCR normalization under abiotic stresses and hormone stimuli in carrot leaves. PLoS ONE 2015, 10, e0117569. [Google Scholar] [CrossRef] [PubMed]
  15. Jian, B.; Liu, B.; Bi, Y.; Hou, W.; Wu, C.; Han, T. Validation of internal control for gene expression study in soybean by quantitative real-time PCR. BMC Mol. Biol. 2008, 9, 59. [Google Scholar] [CrossRef] [PubMed]
  16. Libault, M.; Thibivilliers, S.; Bilgin, D.D.; Radwan, O.; Benitez, M.; Clough, S.J.; Stacey, G. Identification of Four Soybean Reference Genes for Gene Expression Normalization. Plant Genome J. 2008, 1, 44–54. [Google Scholar] [CrossRef]
  17. Fei, X.; Shi, Q.; Yang, T.; Fei, Z.; Wei, A. Expression stabilities of ten candidate reference genes for RT-qPCR in zanthoxylum bungeanum Maxim. Molecules 2018, 23, 802. [Google Scholar] [CrossRef] [PubMed]
  18. Ray, D.L.; Johnson, J.C. Validation of reference genes for gene expression analysis in olive (Olea europaea) mesocarp tissue by quantitative real-time RT-PCR. BMC Res. Notes 2014, 7, 802. [Google Scholar] [CrossRef]
  19. Lucho, S.R.; do Amaral, M.N.; Benitez, L.C.; Milech, C.; Kleinowski, A.M.; Bianchi, V.J.; Braga, E.J.B. Validation of reference genes for RT-qPCR studies in Stevia rebaudiana in response to elicitor agents. Physiol. Mol. Biol. Plants 2018, 24, 767–779. [Google Scholar] [CrossRef]
  20. Li, M.-Y.; Wang, F.; Jiang, Q.; Wang, G.-L.; Tian, C.; Xiong, A.-S. Validation and Comparison of Reference Genes for qPCR Normalization of Celery (Apium graveolens) at Different Development Stages. Front. Plant Sci. 2016, 7, 313. [Google Scholar] [CrossRef]
  21. Wu, Z.J.; Tian, C.; Jiang, Q.; Li, X.H.; Zhuang, J. Selection of suitable reference genes for qRT-PCR normalization during leaf development and hormonal stimuli in tea plant (Camellia sinensis). Sci. Rep. 2016, 6, 19748. [Google Scholar] [CrossRef] [PubMed]
  22. Ni, Z.; Zhou, P.; Xu, M.; Xu, L.A. Development and characterization of chloroplast microsatellite markers for Pinus massoniana and their application in Pinus (Pinaceae) species. J. Genet. 2018, 97, 53–59. [Google Scholar] [CrossRef]
  23. Maleki, S.S.; Mohammadi, K.; Ji, K.S. Study on factors influencing transformation efficiency in Pinus massoniana using Agrobacterium tumefaciens. Plant Cell Tissue Organ Cult. 2018, 133, 437–445. [Google Scholar] [CrossRef]
  24. Ni, Z.X.; Ye, Y.J.; Bai, T.; Xu, M.; Xu, L.A. Complete chloroplast genome of pinus massoniana (pinaceae): Gene rearrangements, loss of ndh genes, and short inverted repeats contraction, expansion. Molecules 2017, 22, 1528. [Google Scholar] [CrossRef] [PubMed]
  25. Xu, C.; Wu, X. Physiological and Proteomic Analysis of Mycorrhizal Pinus massoniana Inoculated with Lactarius insulsus Under Drought Stress. Физиология Растений 2016, 63, 754–762. [Google Scholar] [CrossRef]
  26. Wei, Y.; Liu, Q.; Dong, H.; Zhou, Z.; Hao, Y.; Chen, X.; Xu, L. Selection of Reference Genes for Real-Time Quantitative PCR in Pinus massoniana Post Nematode Inoculation. PLoS ONE 2016, 11, e0147224. [Google Scholar] [CrossRef] [PubMed]
  27. Chen, H.; Yang, Z.; Hu, Y.; Tan, J.; Jia, J.; Xu, H.; Chen, X. Reference genes selection for quantitative gene expression studies in Pinus massoniana L. Trees-Struct. Funct. 2016, 30, 685–696. [Google Scholar] [CrossRef]
  28. Vandesompele, J.; De Preter, K.; Pattyn, F.; Poppe, B.; Van Roy, N.; De Paepe, A.; Speleman, F.; Åman, P.; Semb, H.; Powers, D.; et al. The multifunctional FUS, EWS and TAF15 proto-oncoproteins show cell type-specific expression patterns and involvement in cell spreading and stress response. Genome Biol. 2002, 3, 37. [Google Scholar] [CrossRef]
  29. Andersen, C.L.; Jensen, J.L.; Ørntoft, T.F. Normalization of real-time quantitative reverse transcription-PCR data: A model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004, 64, 5245–5250. [Google Scholar] [CrossRef]
  30. Pfaffl, M.W.; Tichopad, A.; Prgomet, C.; Neuvians, T.P. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper—Excel-based tool using pair-wise correlations. Biotechnol. Lett. 2004, 26, 509–515. [Google Scholar] [CrossRef]
  31. Xie, F.; Xiao, P.; Chen, D.; Xu, L.; Zhang, B. miRDeepFinder: A miRNA analysis tool for deep sequencing of plant small RNAs. Plant Mol. Biol. 2012, 80, 75–84. [Google Scholar] [CrossRef] [PubMed]
  32. Petriccione, M.; Mastrobuoni, F.; Zampella, L.; Scortichini, M. Reference gene selection for normalization of RT-qPCR gene expression data from Actinidia deliciosa leaves infected with Pseudomonas syringae pv. Actinidiae. Sci. Rep. 2015, 5, 16961. [Google Scholar] [CrossRef] [PubMed]
  33. Klie, M.; Debener, T. Identification of superior reference genes for data normalisation of expression studies via quantitative PCR in hybrid roses (Rosa hybrida). BMC Res. Notes 2011, 4, 518. [Google Scholar] [CrossRef] [PubMed]
  34. Kubista, M.; Andrade, J.M.; Bengtsson, M.; Forootan, A.; Jonák, J.; Lind, K.; Sindelka, R.; Sjöback, R.; Sjögreen, B.; Strömbom, L.; et al. The real-time polymerase chain reaction. Mol. Asp. Med. 2006, 27, 95–125. [Google Scholar] [CrossRef] [PubMed]
  35. Bustin, S. Quantification of mRNA using real-time reverse transcription PCR (RT-PCR): Trends and problems. J. Mol. Endocrinol. 2002, 29, 23–39. [Google Scholar] [CrossRef] [PubMed]
  36. Nolan, T.; Hands, R.E.; Bustin, S.A. Quantification of mRNA using real-time RT-PCR. Nat. Protoc. 2006, 1, 1559–1582. [Google Scholar] [CrossRef]
  37. Udvardi, M.K.; Czechowski, T.; Scheible, W.-R. Eleven Golden Rules of Quantitative RT-PCR. Plant Cell Online 2008, 20, 1736–1737. [Google Scholar] [CrossRef] [Green Version]
  38. Jain, M.; Nijhawan, A.; Tyagi, A.K.; Khurana, J.P. Validation of housekeeping genes as internal control for studying gene expression in rice by quantitative real-time PCR. Biochem. Biophys. Res. Commun. 2006, 345, 646–651. [Google Scholar] [CrossRef]
  39. Scholtz, J.J.; Visser, B. Reference gene selection for qPCR gene expression analysis of rust-infected wheat. Physiol. Mol. Plant Pathol. 2013, 81, 22–25. [Google Scholar] [CrossRef]
  40. Tong, Z.; Gao, Z.; Wang, F.; Zhou, J.; Zhang, Z. Selection of reliable reference genes for gene expression studies in peach using real-time PCR. BMC Mol. Biol. 2009, 10, 71. [Google Scholar] [CrossRef]
  41. Xu, M.; Zhang, B.; Su, X.; Zhang, S.; Huang, M. Reference gene selection for quantitative real-time polymerase chain reaction in Populus. Anal. Biochem. 2011, 408, 337–339. [Google Scholar] [CrossRef] [PubMed]
  42. Wang, T.; Lu, J.; Xu, Z.; Yang, W.; Wang, J.; Cheng, T.; Zhang, Q. Selection of suitable reference genes for miRNA expression normalization by qRT-PCR during flower development and different genotypes of Prunus mume. Sci. Hortic. 2014, 169, 130–137. [Google Scholar] [CrossRef]
  43. Barsalobres-Cavallari, C.F.; Severino, F.E.; Maluf, M.P.; Maia, I.G. Identification of suitable internal control genes for expression studies in Coffea arabica under different experimental conditions. BMC Mol. Biol. 2008, 10, 1. [Google Scholar] [CrossRef] [PubMed]
  44. Wang, J.J.; Han, S.; Yin, W.; Xia, X.; Liu, C. Comparison of reliable reference genes following different hormone treatments by various algorithms for qRT-PCR analysis of Metasequoia. Int. J. Mol. Sci. 2019, 20, 34. [Google Scholar] [CrossRef] [PubMed]
  45. De Vega-Bartol, J.J.; Santos, R.R.; Simões, M.; Miguel, C.M. Normalizing gene expression by quantitative PCR during somatic embryogenesis in two representative conifer species: Pinus pinaster and Picea abies. Plant Cell Rep. 2013, 32, 715–729. [Google Scholar] [CrossRef] [PubMed]
  46. Monteiro, F.; Sebastiana, M.; Pais, M.S.; Figueiredo, A. Reference Gene Selection and Validation for the Early Responses to Downy Mildew Infection in Susceptible and Resistant Vitis vinifera Cultivars. PLoS ONE 2013, 8, e72998. [Google Scholar] [CrossRef] [PubMed]
  47. Yang, H.; Liu, J.; Huang, S.; Guo, T.; Deng, L.; Hua, W. Selection and evaluation of novel reference genes for quantitative reverse transcription PCR (qRT-PCR) based on genome and transcriptome data in Brassica napus L. Gene 2014, 538, 113–122. [Google Scholar] [CrossRef]
  48. Zhu, J.; He, F.; Song, S.; Wang, J.; Yu, J. How many human genes can be defined as housekeeping with current expression data? BMC Genomics 2008, 9, 172. [Google Scholar] [CrossRef]
  49. Schmid, H.; Cohen, C.D.; Henger, A.; Irrgang, S.; Schlöndorff, D.; Kretzler, M. Validation of endogenous controls for gene expression analysis in microdissected human renal biopsies. Kidney Int. 2003, 64, 356–360. [Google Scholar] [CrossRef] [Green Version]
  50. Bao, W.; Qu, Y.; Shan, X.; Wan, Y. Screening and validation of housekeeping genes of the root and cotyledon of cunninghamia lanceolata under abiotic stresses by using quantitative real-time PCR. Int. J. Mol. Sci. 2016, 17, 1198. [Google Scholar] [CrossRef]
  51. Tu, Z.; Hao, Z.; Zhong, W.; Li, H. Identification of suitable reference genes for RT-qPCR assays in Liriodendron chinense (Hemsl.) Sarg. Forests 2019, 10, 441. [Google Scholar] [CrossRef]
  52. Mallona, I.; Lischewski, S.; Weiss, J.; Hause, B.; Egea-Cortines, M. Validation of reference genes for quantitative real-time PCR during leaf and flower development in Petunia hybrida. BMC Plant Biol. 2010, 10, 4. [Google Scholar] [CrossRef] [PubMed]
  53. Die, J.V.; Román, B.; Nadal, S.; González-Verdejo, C.I. Evaluation of candidate reference genes for expression studies in Pisum sativum under different experimental conditions. Planta 2010, 232, 145–153. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Distribution of transforming cycle (Ct) values of the twelve candidate reference genes across all experimental samples of P. massoniana. The outside box is determined from 25th to 75th percentiles, and the inside ‘+’ represents the mean value. The line across the box is the median. The whiskers represent the 5th to 95th percentiles, and the asterisks represent the outliers.
Figure 1. Distribution of transforming cycle (Ct) values of the twelve candidate reference genes across all experimental samples of P. massoniana. The outside box is determined from 25th to 75th percentiles, and the inside ‘+’ represents the mean value. The line across the box is the median. The whiskers represent the 5th to 95th percentiles, and the asterisks represent the outliers.
Forests 10 00632 g001
Figure 2. Pairwise variation (V) of the candidate reference genes analyzed by geNorm. Pairwise variation (Vn/Vn+1) was analyzed between the normalization factors (NFn and NFn+1) by geNorm to determine the optimal number of reference genes. Vn/Vn+1 values below 0.15 suggest that there was no need to introduce another gene.
Figure 2. Pairwise variation (V) of the candidate reference genes analyzed by geNorm. Pairwise variation (Vn/Vn+1) was analyzed between the normalization factors (NFn and NFn+1) by geNorm to determine the optimal number of reference genes. Vn/Vn+1 values below 0.15 suggest that there was no need to introduce another gene.
Forests 10 00632 g002
Figure 3. Effect of different reference genes for normalization on the relative expression of the SHR, APS, and PYL genes. (A–C) cold treatment; (D–F) drought treatment; (G–I) ABA treatment; (J–L) salinity treatment; and (M–O) heat treatment.
Figure 3. Effect of different reference genes for normalization on the relative expression of the SHR, APS, and PYL genes. (A–C) cold treatment; (D–F) drought treatment; (G–I) ABA treatment; (J–L) salinity treatment; and (M–O) heat treatment.
Forests 10 00632 g003
Table 1. Candidate reference genes and target genes descriptions, primer sequences, and amplicon characteristics in this study.
Table 1. Candidate reference genes and target genes descriptions, primer sequences, and amplicon characteristics in this study.
Gene SymbolGene DescriptionAccession NumberPrimer Sequence (5′-3′)Amplicon Size (bp)PCR Efficiency (%)Regression Coefficient (R2)Tm (°C)
Reference genes
ACT1Actin 1 geneKM496527.1CCGTATGAGCAAGGAAATCAC10095.6790.99985.0
AGAACCTCCAATCCAGACACT
ACT2Actin 2 geneKM496525.1CACGGAATAGGCAGAAGTTGG9798.1450.98980.8
TGGGCATAAAGTGTTAGAATAGC
ACT4Actin 4 geneKM496528.1ATTTATGAGGGATACGCTTTG10697.5050.99584.6
AGGTGTACCCACGTTCTGTAA
ACT6Actin 6 geneKM496530.1AACTCCTGCCATCCTCATCTT9996.5490.99683.1
CTGTTCCAGCCTTGCTTTCA
ACT7Actin 7 geneKM496529.1TGGGATGCTATGGAAGATTTG114104.8010.99282.9
TACGCCCTTTGGAGTAAGAAG
AQPAquaporin protein geneKF582038.1CACCTTGCCACAATTCCTATCA10395.4810.99886.3
TCCAATGGTCATCCCAAACAC
CYPCyclophilin geneKM496534.1CGAGAAGTTTGCCGATGAGAA9791.0420.99787.5
GAATTGCGAGCCGTTAGTGTT
EF1AElongation factor 1-alpha geneKM496532.1GGATTTGAAACGTGGGTATGT9799.3830.99883.5
CAGGGTGGTTCATTATGATTACT
F-boxF-box family protein geneKM496542.1TATTATTGTTGCAGGTGGGTT109100.1010.99681.6
AGAATGTTGAAGTTCGGCTAT
GAPDHGlyceraldehyde 3-Phosphatase geneKM496531.1GGATTTGGTCGTATTGGGAGG9692.1320.99882.9
TTTGGCATCAATGAAAGGGTC
TUATubulin alpha geneKM496535.1CAAACTTGGTCCCGTATCCTC9592.0110.99983.7
CACAGAAAGCTGCTCATGGTAA
TUBTubulin beta geneKM496536.1CTGCGACTATGAGTGGAGTGA10898.9230.99385.6
AGAAATGAAGACGAGGGAATG
Target genes
SHRShort-Root geneMK153765GCCTGTGAGGATTCTGAAGTT9797.6360.99685.3
CACTGAAAGCAGCATGTATGA
APSAlpha-pinene synthase geneKF547035TGGATCGCCAGTGGTGAGGTG10594.6940.99186.2
GTCGGTCGTCAGAATGGGTTG
PYLPyrabactin resistance-like geneMK953936GAGTCCGAGTATGTGTGGAGGC110100.0510.99786.2
ACTAATGACCAAACCAGATGAA
Table 2. The stability ranking of candidate reference genes by geNorm, NormFinder, BestKeeper, and RefFinder. SD [±Ct], standard deviation of the Ct; CV [%Ct], coefficient of variance expressed as a percentage of the Ct level.
Table 2. The stability ranking of candidate reference genes by geNorm, NormFinder, BestKeeper, and RefFinder. SD [±Ct], standard deviation of the Ct; CV [%Ct], coefficient of variance expressed as a percentage of the Ct level.
TreatmentRankgeNormNormFinderBestKeeperRefFinder
GeneStabilityGeneStabilityGeneSD(±Ct)CV(%Ct)GeneStability
Salinity1ACT20.275 F-box0.145 EF1A0.30 1.24 ACT21.41
2F-box0.275 ACT20.164 ACT20.34 1.13 F-box1.68
3CYP0.315 ACT40.177 CYP0.39 1.69 CYP3.46
4ACT40.410 CYP0.264 F-box0.49 1.74 ACT44.36
5ACT70.453 ACT10.266 ACT70.54 1.91 ACT15.48
6ACT10.502 ACT70.347 ACT10.58 2.53 ACT75.48
7ACT60.549 TUA0.398 GAPDH0.60 2.58 EF1A6.04
8TUA0.583 ACT60.405 ACT60.66 2.15 ACT67.48
9GAPDH0.633 TUB0.501 TUA0.76 3.28 TUA7.97
10TUB0.687 GAPDH0.572 ACT40.78 2.81 GAPDH8.91
11EF1A0.734 EF1A0.588 TUB1.03 3.75 TUB9.72
12AQP0.787 AQP0.656 AQP1.12 4.67 AQP12.00
ABA1CYP0.246 CYP0.085 GAPDH0.42 1.78 CYP1.41
2TUA0.246 TUA0.142 F-box0.42 1.51 TUA2.38
3ACT60.319 AQP0.193 ACT10.47 2.07 ACT13.94
4ACT70.345 ACT10.202 CYP0.53 2.29 AQP4.24
5ACT10.361 ACT70.257 EF1A0.61 2.52 GAPDH4.30
6AQP0.375 ACT60.267 AQP0.62 2.67 ACT75.14
7GAPDH0.423 GAPDH0.282 ACT70.69 2.43 ACT65.73
8ACT20.466 ACT20.321 TUA0.70 3.03 F-box7.18
9ACT40.507 ACT40.400 ACT20.75 2.49 ACT28.24
10EF1A0.556 EF1A0.453 ACT60.75 2.44 EF1A8.41
11F-box0.611 F-box0.578 ACT40.75 2.74 ACT49.46
12TUB0.663 TUB0.580 TUB0.99 3.67 TUB12.00
Drought1ACT40.280 ACT70.093 F-box0.24 0.84 ACT71.57
2TUB0.280 CYP0.223 ACT70.31 1.06 CYP2.63
3ACT70.360 TUA0.264 CYP0.34 1.43 ACT43.36
4CYP0.427 ACT40.296 AQP0.39 1.56 TUA4.24
5ACT20.458 ACT20.328 ACT20.45 1.45 TUB4.36
6TUA0.490 TUB0.344 TUA0.50 2.16 F-box4.45
7F-box0.557 GAPDH0.398 EF1A0.51 2.04 ACT25.00
8AQP0.600 F-box0.403 ACT40.52 1.84 AQP7.14
9EF1A0.640 AQP0.466 GAPDH0.56 2.36 GAPDH8.43
10GAPDH0.679 ACT10.512 TUB0.60 2.16 EF1A9.34
11ACT10.722 EF1A0.555 ACT10.74 3.21 ACT110.49
12ACT60.785 ACT60.692 ACT60.83 2.71 ACT612.00
Cold1ACT10.157 GAPDH0.111 ACT10.24 0.99 ACT11.57
2TUA0.157 ACT70.133 TUA0.29 1.20 ACT72.34
3ACT70.232 ACT10.145 ACT40.30 1.04 TUA2.51
4GAPDH0.281 ACT20.150 F-box0.30 1.06 GAPDH2.91
5ACT40.325 TUA0.215 ACT70.35 1.18 ACT44.82
6ACT20.350 ACT40.224 GAPDH0.37 1.54 ACT25.89
7ACT60.394 AQP0.283 AQP0.39 1.64 F-box6.93
8AQP0.427 F-box0.284 CYP0.40 1.72 AQP7.24
9F-box0.453 ACT60.330 ACT60.42 1.34 ACT68.45
10CYP0.476 CYP0.395 ACT20.45 1.43 CYP9.46
11TUB0.523 TUB0.455 TUB0.70 2.47 TUB11.00
12EF1A0.634 EF1A0.782 EF1A0.87 3.41 EF1A12.00
Heat1CYP0.309 ACT10.187 TUA0.50 2.06 ACT12.30
2GAPDH0.309 TUA0.275 ACT40.55 1.92 CYP2.63
3ACT20.416 CYP0.284 ACT10.58 2.39 TUA3.13
4ACT70.496 TUB0.325 TUB0.89 3.02 GAPDH3.64
5AQP0.589 GAPDH0.329 AQP0.94 3.58 ACT25.05
6TUB0.635 ACT20.356 GAPDH0.96 3.96 TUB5.18
7ACT10.658 AQP0.405 ACT20.99 3.10 F-box5.62
8TUA0.686 ACT40.419 F-box1.13 4.07 AQP6.19
9ACT40.724 ACT70.482 ACT61.23 3.96 ACT46.45
10F-box0.798 F-box0.682 ACT71.32 4.35 ACT77.54
11EF1A0.874 EF1A0.847 EF1A1.38 5.24 EF1A11.24
12ACT61.019 ACT61.145 CYP1.46 5.82 ACT611.74
Total1ACT10.437 ACT20.267 F-box0.44 1.57 TUA2.11
2TUA0.437 TUA0.306 CYP0.59 2.53 ACT22.30
3ACT70.625 ACT10.316 GAPDH0.65 2.76 ACT12.71
4ACT20.650 CYP0.335 ACT40.72 2.57 CYP3.72
5GAPDH0.671 ACT70.358 TUA0.75 3.19 ACT75.10
6CYP0.686 ACT40.372 ACT10.77 3.31 GAPDH5.21
7ACT40.703 GAPDH0.394 ACT20.77 2.50 F-box5.62
8TUB0.747 TUB0.508 ACT60.85 2.76 ACT45.63
9EF1A0.778 EF1A0.541 ACT70.89 3.03 TUB8.66
10F-box0.827 F-box0.580 EF1A1.00 3.96 EF1A9.24
11AQP0.882 AQP0.698 TUB1.16 4.16 ACT610.84
12ACT60.932 ACT60.711 AQP1.22 4.99 AQP11.24

Share and Cite

MDPI and ACS Style

Zhu, P.; Ma, Y.; Zhu, L.; Chen, Y.; Li, R.; Ji, K. Selection of Suitable Reference Genes in Pinus massoniana Lamb. Under Different Abiotic Stresses for qPCR Normalization. Forests 2019, 10, 632. https://doi.org/10.3390/f10080632

AMA Style

Zhu P, Ma Y, Zhu L, Chen Y, Li R, Ji K. Selection of Suitable Reference Genes in Pinus massoniana Lamb. Under Different Abiotic Stresses for qPCR Normalization. Forests. 2019; 10(8):632. https://doi.org/10.3390/f10080632

Chicago/Turabian Style

Zhu, Peihuang, Yinyan Ma, Lingzhi Zhu, Yu Chen, Rong Li, and Kongshu Ji. 2019. "Selection of Suitable Reference Genes in Pinus massoniana Lamb. Under Different Abiotic Stresses for qPCR Normalization" Forests 10, no. 8: 632. https://doi.org/10.3390/f10080632

APA Style

Zhu, P., Ma, Y., Zhu, L., Chen, Y., Li, R., & Ji, K. (2019). Selection of Suitable Reference Genes in Pinus massoniana Lamb. Under Different Abiotic Stresses for qPCR Normalization. Forests, 10(8), 632. https://doi.org/10.3390/f10080632

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop