Selection and Validation of Appropriate Reference Genes for Real-Time Quantitative PCR Analysis in Needles of Larix olgensis under Abiotic Stresses

: Larix olgensis Henry is an important a ﬀ orestation species in northeastern China because of its fast juvenile growth, high-quality timber, and signiﬁcant economic and ecological values. The selection of appropriate reference genes is necessary for the normalization of gene expression determination during quantitative real-time polymerase chain reaction (qRT-PCR) experiments. In this study, qRT-PCR was used to study gene expression. Three software packages geNorm, NormFinder, BestKeeper were used, and a comprehensive ranking of candidate reference genes was produced based on their output to evaluate the expression stability of 16 candidate reference genes from L. olgensis under drought, salt, cold, and heat stress. PP2A-1 and GAPDH ranked as the most stable reference genes under drought and cold stress, PP2A-1 and UBQ10 were most stable under salt stress, and TIP41 and ACT2 were most stable under heat stress. The least stable gene was ADP , which ranked the last under all treatments. Expression proﬁle analysis of the antioxidant gene CAT using the two most stable and the single least stable reference genes under each stress further veriﬁed that the selected reference genes were suitable for gene expression normalization. This study provides an important foundation for the selection of suitable reference genes for the normalization and quantiﬁcation of L. olgensis gene expression under abiotic stress conditions.


Introduction
Gene expression analysis is an important tool for identifying key genes and understanding complex biological processes such as metabolic pathways, signal transduction, and plant development. Quantitative real-time polymerase chain reaction (qRT-PCR) is the most effective, simple, specific, inexpensive, and sensitive method for quantifying the expression of target genes [1,2]. Nevertheless, various factors such as RNA integrity, reverse transcription efficiency, cDNA quality, primer specificity, amplification efficiency, and the selection of reference genes (RGs) may significantly influence the reliability of qRT-PCR results [3,4]. Selection of inappropriate RGs will introduce inaccuracies into the experimental data, and screening for one or more suitable RGs is therefore important for the normalization of gene expression data.
Previous literature has described the selection of RGs for various species under different biotic and abiotic stresses and in different development stages and tissues. These species include Arabidopsis thaliana [5], Oryza sativa [6], Solanum lycopersicum [7], Malus domestica [8], and Populus euramericana cv [9]. RGs are often housekeeping genes that are associated with basic cellular processes and Forests 2020, 11,193 3 of 17

Plant Materials and Treatments
Mature seeds of L. olgensis Henry were collected from the Qingshan Forestry Bureau Seed Orchard in Heilongjiang Province and sown in plastic pots (11 × 11 cm) containing a vermiculite/soil mixture (1:1) and cultured in a growth chamber with 70% relative humidity under a 16-h/8-h light/dark photoperiod with a light intensity of 150 µmol m -2 s −1 [27]. Three-month-old seedlings with consistent growth were selected for abiotic treatments. For salt and drought experiments, pots were irrigated with 0.2 M NaCl or 20% PEG 6000 solutions for 0, 2, 6, 12, and 24 h. Plants were irrigated every 12 h during the experiment. For heat and cold stress treatments, three-month-old seedlings were transferred into light incubators and exposed to 40 or 4 • C for 0, 2, 6, 12, and 24 h [32]. All the treatments were performed with three biological replicates. Needles were carefully harvested from treated and untreated plants, frozen immediately in liquid nitrogen, and stored at −80 • C for RNA extraction.

Selection of Candidate Reference Genes and Primer Design
This study took advantage of the unpublished L. olgensis genome information generated by our team: we had previously analyzed transcriptome for this species with and without watering (unpublished). We screened candidate reference genes according to the following criteria: the fragments per kilobase of exon model per million mapped reads (FPKM values) were appropriating in all samples, and the coefficient of variation (CV) of FPKM value was cut off less than 0.5. Based on previous studies the criteria, 16 commonly used reference genes including glyceraldehyde 3-phosphate dehydrogenase (GAPDH), 18S ribosomal RNA (18S), actin 2 (ACT2), tubulin beta-6 (TUB), eukaryotic translation initiation factor 4α-1 (eIF-4α), Ef 1alpha (EF-1α), tubulin alpha-2 (TUA), ubiquitin-conjugating enzyme 9 (UBC9), TIP41-like protein (TIP41), protein phosphatase 2A-1 (PP2A-1), polyubiquitin 10 (UBQ10), polypyrimidine tract binding protein (PTBP1), ADP-ribosylation factor (ADP, or ARF), histone (HIS), ubiquitin-like protein RUB2 (UBQ7), and actin protein coding 12 (ACT12) genes were used as candidate genes to identify the most stable RGs under different treatments. Detailed information on each gene is presented in Table 1. The coding sequences (CDs) of the 16 Arabidopsis genes were downloaded from the TAIR database (http://www.arabidopsis.org) to identify their homologs in the L. olgensis transcriptome. The Bioedit Sequence Alignment Editor was used to perform a local BLAST to conduct the blastn search of the L. olgensis transcriptome using the Arabidopsis query sequences. We identified 16 candidate RGs and one target gene, CAT ( Table 2). The CDs sequences from the L. olgensis transcriptome were submitted to GenBank. Primer5 was used to design specific primers for each gene using the following criteria: GC content 44%-60%, optimal Tm 58-60 • C, primer length 20-22 bp, and amplicon length 80-220 bp ( Table 2). The specificity of all primer pairs was checked by standard PCR using cDNA as a template with KOD FX (Toyobo, Osaka, Japan), and amplified products were verified with 2% agarose gel and sequenced to confirm their identity.

Reverse Transcription Quantitative PCR (qRT-PCR) Analysis
Quantitative RT-PCR reactions were conducted in 96-well plates with a qTOWER 3G Cycler and qPCR software (Analytik Jena, Jena Germany) using the TransStart Top Green qPCR SuperMix (TransGen Biotech, Beijing, China). The 20 µL reaction mix contained 10 µL TransStart Top Green qPCR SuperMix, 7 µL nuclease-free water, 1 µL diluted cDNA, and 1 µL of each specific primer (final concentration 10 µM). PCR conditions were 94 • C for 30 s, followed by 45 cycles of 95 • C for 5 s, 59 • C for 15 s, and 72 • C for 10 s. Three technical replicates were performed for each sample to ensure the accuracy of the results. All primers used in the study are listed in Table 2.

Gene Expression Stability Analysis
Three different Microsoft Excel-based software programs, geNorm [3], NormFinder [15], and BestKeeper [13,34], were used to analyze the expression stability of the candidate RGs under different experimental conditions. The Cq (PCR cycle threshold) data were used directly in the BestKeeper program, but for geNorm and NormFinder, Cq values were converted into relative values and imported into geNorm and NormFinder program analyzed as described previously [13]. The output of each software program permitted us to rank the expression stability of RGs in different treatment groups. A comprehensive ranking of RGs was also generated as described previously [35,36].

Validation of RGs by qRT-PCR
QRT-PCR was also performed to analyze the expression levels of the target gene CAT under different experimental conditions. The top two best RGs and worst ranked RG for each experimental condition were selected normalization of CAT expression in L. olgensis. The quantitative variation between replicates was calculated with the relative quantification method (2 −∆∆CT ) [37]. Graphs were generated using Excel and GraphPad Prism7. The primer for CAT is listed in Table 2.

Selection of Candidate Reference Genes
Based on sequence homology with Arabidopsis genes, 16 candidate RGs were identified from the L. olgensis transcriptome (unpublished) ( Table 1). Subsequently, primer specificity was confirmed based on agarose gel electrophoresis. Amplicons of a single band with the expected size indicate good primer specificity ( Figure S1). In qRT-PCR reactions, melting curve analysis of every candidate RG showed a single peak. Amplification efficiency (E) and correlation coefficient values (R 2 ) of the RG standard curves varied from 96.85 to 108.72% and from 0.9918 to 0.9997, respectively (Table 2 and Figure S2).

Expression Profiles of Candidate RGs
The distribution of raw Cq values for the 16 candidate RGs in 20 L. olgensis samples (five samples from each stress condition) is shown in Figure 1 (Table 2).
Forests 2020, 11,193 7 of 18 heat stress. The ranking of gene stability by average CV values across all treatments was PP2A (Table 2).     (Table 2).

Analysis of Gene Expression Stability
The software programs geNorm, NormFinder, and BestKeeper were used to analyze the expression stability of the 16 candidate RGs in L. olgensis plants that had been subjected to four experimental treatments. Data from each treatment were analyzed separately and in combination. geNorm analysis: for geNorm analysis, the M-value of candidate RGs was calculated and used to rank RG stability. A gene with an M-value below 1.5 is considered to be stably expressed, and the smaller the M-value, the more stable the gene [3]. In our study, the M-values of all candidate RGs were lower than 1.5 in individual treatments and when the treatments were combined, with the exception of ADP under heat stress. The most suitable RG differed among treatments. As shown in Table 3 and Figure S3, GAPDH and 18S under drought stress, UBQ10 and PP2A-1 under salt stress, PP2A-1 and GAPDH under cold stress, and TIP41 and ACT2 under heat stress were ranked as the most stable genes with the lowest M-values. PP2A-1 and ACT2 had the most stable expression when data from all samples were combined, whereas, ADP was the least stable gene in all treatment conditions.  For qRT-PCR, the selection of a greater number of RGs can permit more accurate quantification of target gene expression. geNorm was used to calculate the pairwise variation (Vn/n+1), which permits the determination of an optimal number of RGs for each treatment group using a threshold value of 0.15. There is no need for additional RGs if the variation value is below 0.15 [3]. As shown in Figure 3, the pairwise variation values V 2 /V 3 for drought, salt, and cold stress samples were all less than 0.15, indicating that two suitable RGs were adequate for the normalization of data from these treatments. For the heat stress samples, V4/5 was 0.1498, indicating that four RGs (TIP41, ACT2, PTBP1, and eIF-4α) were required. When all samples were combined, V3/4 was 0.1315, indicating that three RGs (ACT2, PP2A-1, and TIP41) were required.
NormFinder analysis: the NormFinder approach provides a stability value for each gene based on inter-and intra-group variation in expression [15]. As shown in Table 4, the stability ranks calculated in NormFinder were consistent with those calculated in geNorm under cold stress (GAPDH). The top three most stable genes were GAPDH, 18S, and PP2A-1 under drought stress; ACT2, 18S, UBC9, PPA2-1, and UBQ10 under salt stress; and 18S, PTBP1, and ACT2 under heat stress. These rankings were broadly similar to those of geNorm, although there were some differences. For instance, in the stability rankings for the combined samples, NormFinder suggested that 18S and UBC9 were the two most suitable two RGs, a result that differed from that of geNorm. Similar to geNorm, NormFinder indicated that ADP was the least stable gene in all treatments. Although there were some differences in Forests 2020, 11,193 9 of 17 rankings between geNorm and NormFinder, the top five most stable genes were relatively consistent between them. NormFinder analysis: the NormFinder approach provides a stability value for each gene based on inter-and intra-group variation in expression [15]. As shown in Table 4, the stability ranks calculated in NormFinder were consistent with those calculated in geNorm under cold stress (GAPDH). The top three most stable genes were GAPDH, 18S, and PP2A-1 under drought stress; ACT2, 18S, UBC9, PPA2-1, and UBQ10 under salt stress; and 18S, PTBP1, and ACT2 under heat stress. These rankings were broadly similar to those of geNorm, although there were some differences. For instance, in the stability rankings for the combined samples, NormFinder suggested that 18S and UBC9 were the two most suitable two RGs, a result that differed from that of geNorm. Similar to geNorm, NormFinder indicated that ADP was the least stable gene in all treatments. Although there were some differences in rankings between geNorm and NormFinder, the top five most stable genes were relatively consistent between them. BestKeeper analysis: another method for detection of suitable RGs is BestKeeper. Unlike the other two programs, BestKeeper ranks RGs based on the standard deviation (SD) and coefficient of  BestKeeper analysis: another method for detection of suitable RGs is BestKeeper. Unlike the other two programs, BestKeeper ranks RGs based on the standard deviation (SD) and coefficient of variation (CV) of the Cq values from the qRT-PCR assay. The smaller the CV, the better the stability of the gene was [13,34]. The results of BestKeeper analysis are also listed in Table 5. PPA2-1 was ranked first under salt stress and cold stress, with CV ± SD values of 0.17 ± 0.65 and 0.40 ± 1.51, respectively. GAPDH was ranked first under drought stress, with a CV ± SD of 0.17 ± 0.79 and UBQ7 was ranked first under heat stress, with a CV ± SD of 0.29 ± 1.22. BestKeeper suggested that ACT2 was the most suitable RG when data from all samples were combined. Few genes had an SD greater than 1.0, indicating that most of the candidate RGs were relatively stable. Similar to the results of geNorm and NormFinder, ADP was the least stable RG for all sample sets, with the exception of the combined samples, in which EF-1α was the least stable.

Comprehensive Analysis
To reduce the effect of any limitations and biases associated with individual algorithms, a comprehensive stability analysis was performed by taking the geometric mean of the geNorm, NormFinder, and BestKeeper rankings in order to identify the best RGs [38]. For comprehensive analysis, two RGs were selected for further normalization. As shown in Table 6 and Figure 4, PP2A-1 was ranked among of the top two most stable RGs under drought, salt, and cold stress, and GAPDH under cold stress, UBQ10 under salt stress were also the most stable RGs, respectively. TIP41 and ACT2 were the most stable RGs under heat stress, and ACT2 and PP2A-1 were the most stable RGs when all samples were combined. The comprehensive RG rankings for single treatments (drought, salt, and cold) and for combined samples were consistent with the results obtained by geNorm and BestKeeper. ADP was the least stable RG under all experimental conditions (Table S1).

Validation of CAT Reference Gene
To verify the reliability of the selected reference genes, CAT was selected as a target gene. The expression of the antioxidant CAT gene is induced by many abiotic stresses, including chilling, drought, osmotic stress, and salt stress [39][40][41]. The CAT sequence was obtained from L. olgensis transcriptome data and showed 48.75% nucleotide identity with CAT from Arabidopsis. We used the comprehensive ranking results to select the top two stable RGs for normalization of CAT expression in L. olgensis needles. The selected RGs were PP2A-1 and GAPDH for drought and cold stress, PP2A-1 and UBQ10 for salt stress, and TIP41 and ACT2 for heat stress. As shown in Figure 5, the expression patterns of CAT showed few differences when different stable RGs were used. Similar expression patterns were obtained under drought and cold stress when GAPDH and PP2A-1 were used for normalization. CAT expression reached a peak at 12 h for both drought stress (<29-fold) and cold stress (<5-fold). When using the most stable RGs for salt stress PP2A-1 and UBQ10, CAT expression patterns were consistent, and the highest expression was observed at 24 h (no more than 10-fold). Under heat stress, CAT expression levels reached a similar peak at 6 h when normalized with TIP41 and ACT2. However, large differences in expression patterns were detected when the least stable RG, ADP, was used for normalization under the same stresses. Specifically, the expression levels of CAT were overestimated.

Validation of CAT Reference Gene
To verify the reliability of the selected reference genes, CAT was selected as a target gene. The expression of the antioxidant CAT gene is induced by many abiotic stresses, including chilling, drought, osmotic stress, and salt stress [39][40][41]. The CAT sequence was obtained from L. olgensis transcriptome data and showed 48.75% nucleotide identity with CAT from Arabidopsis. We used the comprehensive ranking results to select the top two stable RGs for normalization of CAT expression in L. olgensis needles. The selected RGs were PP2A-1 and GAPDH for drought and cold stress, PP2A-1 and UBQ10 for salt stress, and TIP41 and ACT2 for heat stress. As shown in Figure 5, the expression patterns of CAT showed few differences when different stable RGs were used. Similar expression patterns were obtained under drought and cold stress when GAPDH and PP2A-1 were used for normalization. CAT expression reached a peak at 12 h for both drought stress (<29-fold) and cold stress (<5-fold). When using the most stable RGs for salt stress PP2A-1 and UBQ10, CAT expression patterns were consistent, and the highest expression was observed at 24 h (no more than 10-fold). Under heat stress, CAT expression levels reached a similar peak at 6 h when normalized with TIP41 and ACT2. However, large differences in expression patterns were detected when the least stable RG, ADP, was used for normalization under the same stresses. Specifically, the expression levels of CAT were overestimated.

Discussion
Quantitative RT-PCR is an important technology that permits rapid and reliable quantification of gene expression [41]. Due to its high sensitivity, many experimental variables can easily affect the reliability of qRT-PCR gene expression results [4]. Therefore, to obtain more reliable results, an appropriate normalization strategy is very important. Among several methods [1,2], selection of one or more RGs as normalization factor(s) is the most common approach for different experimental conditions such as abiotic stress [38,42,43]. However, no research on RGs for L. olgensis under abiotic stress has been reported. In this study, we identified appropriate L. olgensis RGs under different abiotic stress conditions using qRT-PCR and parallel calculations in geNorm, Normfinder, Bestkeeper, and comprehensive analysis.
Sixteen potential RGs were identified from L. olgensis transcriptome data generated in our laboratory. Three statistical algorithms, geNorm, NormFinder, and BestKeeper, were used to assess the expression stability of candidate RGs for accurate normalization in gene expression studies [20]. In the geNorm analysis, M-values of the RGs under different experimental conditions were all less than 1.5 (except ADP in heat stress), indicating their potential expression stability [44]. The candidate RGs exhibited differential stability and relative rankings in response to different stresses. In two treatments, the most stable genes calculated by geNorm and Bestkeeper were consistent, namely GAPDH in drought and PP2A-1 in cold (Tables 3 and 5). In salt and heat stress, the most stable genes calculated by the three algorithms differed. For heat stress, TIP41 and ACT2 were ranked as the most stable RGs according to the M-values calculated in geNorm. NormFinder regarded 18S as the most stable RG, and BestKeeper identified UBQ7 as the most stable RG; however, UBQ7 was ranked lower by geNorm and Normfinder (Tables 3 and 4). The differences in ranking were mainly caused by variations in the algorithms of the three programs. Previous work has also shown that the three programs generate different results under abiotic stress in Vitis vinifera [45]. A comprehensive analysis is therefore necessary to provide ultimate stability rankings for RGs under different treatments, as reported previously [46,47].

Discussion
Quantitative RT-PCR is an important technology that permits rapid and reliable quantification of gene expression [41]. Due to its high sensitivity, many experimental variables can easily affect the reliability of qRT-PCR gene expression results [4]. Therefore, to obtain more reliable results, an appropriate normalization strategy is very important. Among several methods [1,2], selection of one or more RGs as normalization factor(s) is the most common approach for different experimental conditions such as abiotic stress [38,42,43]. However, no research on RGs for L. olgensis under abiotic stress has been reported. In this study, we identified appropriate L. olgensis RGs under different abiotic stress conditions using qRT-PCR and parallel calculations in geNorm, Normfinder, Bestkeeper, and comprehensive analysis.
Sixteen potential RGs were identified from L. olgensis transcriptome data generated in our laboratory. Three statistical algorithms, geNorm, NormFinder, and BestKeeper, were used to assess the expression stability of candidate RGs for accurate normalization in gene expression studies [20]. In the geNorm analysis, M-values of the RGs under different experimental conditions were all less than 1.5 (except ADP in heat stress), indicating their potential expression stability [44]. The candidate RGs exhibited differential stability and relative rankings in response to different stresses. In two treatments, the most stable genes calculated by geNorm and Bestkeeper were consistent, namely GAPDH in drought and PP2A-1 in cold (Tables 3 and 5). In salt and heat stress, the most stable genes calculated by the three algorithms differed. For heat stress, TIP41 and ACT2 were ranked as the most stable RGs according to the M-values calculated in geNorm. NormFinder regarded 18S as the most stable RG, and BestKeeper identified UBQ7 as the most stable RG; however, UBQ7 was ranked lower by geNorm and Normfinder (Tables 3 and 4). The differences in ranking were mainly caused by variations in the algorithms of the three programs. Previous work has also shown that the three programs generate different results under abiotic stress in Vitis vinifera [45]. A comprehensive analysis is therefore necessary to provide ultimate stability rankings for RGs under different treatments, as reported previously [46,47].
Our comprehensive ranking results showed that PP2A-1 and GAPDH were the top two most stable RGs under drought and cold stress. PP2A-1 also ranked first under salt stress in the comprehensive analysis. Our findings confirm several previous studies in which PP2A-1 was also shown to be the most stable RG. These include studies of Pennisetum glaucum under different abiotic stresses and hormonal stimuli [48], Agrostis stolonifera under different abiotic stresses [49], P. euphratica under drought stress [50], and M. charantia under salt stress [51]. GAPDH also ranked as a highly suitable RG under drought and cold stress in this study. This result is in agreement with previous research on Caragana korshinskii under heat stress [46] and Coffea arabica leaves under drought stress and GA, SA, and MeJA treatments [52]. However, GAPDH was ranked as the least stable RG in various tissues and under abiotic stresses in Sorghum bicolor and in various tissues and under PEG6000 and MeJA treatments in Peucedanum praeruptorum [53,54]. These results demonstrate that there is no universal RG that is stably expressed in all treatments and tissues.
UBQ10 was also the most stable RG pair with PP2A-1 under salt stress. UBQ10 was also the most constitutively expressed polyubiquitin gene in A. thaliana across all development processes [44], in Prunus persica fruit at multiple developmental stages [55] and in Platycladus orientalis under NaCl and ABA treatments [56]. However, in O. sativa, UBQ10 was the least stable RG in different tissues, cell types, and developmental stages, a result similar to that observed in Glycine max [57]. In the present study, TIP41 and ACT2 showed highly stable expression under heat stress. In a previous work, TIP41 was a stable RG in M. charantia under UV and CuSO 4 treatments [51] and in P. praeruptorum under NaCl stress [54]. No previous research has shown that TIP41 is the most stable RG under heat stress. ACT2 is a traditional housekeeping gene widely used in many species, such as Betula platyphylla under salt stress and Panicum virgatum L. leaves and roots [58,59]. In this study, ACT2 was an appropriate RG under heat stress, consistent with results from P. praeruptorum under heat stress [54]. Nonetheless, geNorm, NormFinder, and Bestkeeper all showed that ADP was the least stable RG under different experimental conditions. Previous studies have shown that ADP is the most stably RG in many species, such as Triticeae [60], Triticum aestivum L [61], and Swingle citrumelo [62]. However, ADP was recognized as the least stable gene by three algorithms for all samples in this study. Several papers have reported similar results. For example, GAPDH ranked the worst in T. aestivum but showed the most stability in V. vinifera cv. Cabernet Sauvignon [63,64]; UBI and ACT showed instability in S. lycopersicum [7], but performed suitably in T. aestivum [63]. These results indicate that the expression stabilities of RGs vary in different tissues and under environmental stresses.
Plant growth and development are generally affected by abiotic stresses (cold, heat, drought, and high salinity), and previous research has shown that CAT is the predominant enzyme controlling H 2 O 2 levels. CAT has often been used as an abiotic-stress-inducible gene, upregulated by drought [40], cold [65], and salt [66] treatments. CAT1 exhibited increased expression in response to exogenous ABA and osmotic stress in maize [39], and its expression was activated by drought, ABA and salt stress in A. thaliana [67]. In this study, CAT expression was upregulated gradually and reached a maximum value at 24 h (≈10-fold upregulated) under salt stress when using the two most stable RGs, PP2A-1 and UBQ10, for normalization. The CAT expression pattern was consistent with previously reported responses of PgCAT1 to salt stress, which showed the highest PgCAT1 expression at 24 h [66] and the highest transcript levels as 10-fold [14]. However, the expression level of CAT was overestimated when the least stable RG, ADP, was used for normalization ( Figure 5). Like salt stress, large differences were detected in CAT expression pattern when the least RG rather than the two most stable RGs was used for normalization under drought, cold, and heat stress ( Figure 5). These results demonstrate that it is necessary to select a suitable stable RG for normalization of target gene expression under different conditions.
In this study, needles were used as the single material to test RGs in L. olgensis. The identified RGs may not apply to a broad range of developmental tissue samples of L. olgensis. For further studies, we will focus on exploring appropriate reference genes among different tissues and different developmental stages for accurate determination of gene expression in L. olgensis in further related research.

Conclusions
We selected 16 candidate RGs to validate suitable RGs using three statistical algorithms (geNorm, NormFinder, BestKeeper, and comprehensive analysis) for gene expression normalization in L. olgensis under drought, salt, cold, and heat stress. The results were compared and ranked using a comprehensive analysis. Based on the comprehensive analysis of gene stability, we identified PP2A-1 and GAPDH as the most stable RGs under drought and cold stresses and PP2A-1 and UBQ10 as the two most stable RGs under salt stress. TIP41 and ACT2 were the most stable RGs under heat stress, and ADP was the least stable RG under all stresses. Furthermore, the expression profiles of CAT confirmed the importance of using two suitable RGs rather than one unstable RG under different stresses. Selection of suitable RGs provides a foundation for functional genomic studies in L. olgensis and other woody plants under abiotic stress conditions. Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/2/193/s1, Figure S1: Agarose gel electrophoresis of PCR products for each of the 16 RGs and the one target genes in Larix olgensis., Figure S2: Melting curve for 16 RGs and one target gene CAT in L. olgensis, Figure S3: Average expression stability values (M-values) of the 16 RGs in L. olgensis calculated by geNorm. Table S1: Expression stability ranking of the 16 candidate reference genes under different abiotic stresses in L. olgensis.