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

: 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 di ﬀ erent 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 di ﬀ erent 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 signiﬁcant 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 beneﬁcial for further studies on gene transcription in early-stage, unligniﬁed seedlings of P. massoniana. and 48 h. For drought, salinity, and ABA treatments, washed seedlings were immersed in complete medium containing 10% (w / v) PEG 6000 , 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.


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.

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-fiveday-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) PEG 6000 , 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°C 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.

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.

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 (R 2 ) [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 ddH 2 O, 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.

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 ddH 2 O. 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 ddH 2 O and total RNA to replace the cDNA.

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].

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].

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 R 2 varying from 0.989 to 0.999 (Table 1 and Appendix A Figure A3).

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. 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.

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.
Ct Value C C 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.

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.

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. 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).

Figure 2.
Pairwise variation (V) of the candidate reference genes analyzed by geNorm. Pairwise variation (V n /V n+1 ) was analyzed between the normalization factors (NF n and NF n+1 ) by geNorm to determine the optimal number of reference genes. V n /V n+1 values below 0.15 suggest that there was no need to introduce another gene.

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.

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).

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).

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).

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  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.

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.

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.  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.

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.          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.