Selection of Suitable Reference Genes for Gene Expression Normalization Studies in Dendrobium huoshanense

Dendrobium huoshanense is a kind of precious herb with important medicinal and edible value in China, which is widely used in traditional Chinese medicine for various diseases. Recent studies have paid close attention to the genetic expression of the biosynthetic pathway of the main active components (polysaccharides, alkaloids, and flavonoids), and real-time polymerase chain reaction (qPCR) is one of the most widely used methods for doing so. However, so far, no reference gene selections have been reported in D. huoshanense. In this study, 15 reference gene candidates (GAPDH, eIF, EF-1α, PP2A, UBCE, RPL5, TBP, APT1, MDH, PTBP3, PEPC, CYP71, NCBP2, TIP41, and F-box) were selected and evaluated for their expression stability in D. huoshanense under various experimental conditions, including in different tissues (root, stem, and leaf), abiotic stresses (oxidative, drought, cold, and UV), and hormone treatment (methyl jasmonate) using three statistical programs (geNorm, NormFinder, and BestKeeper). Then, the RefFinder program was employed to comprehensively validate the stability of the selected reference genes. Finally, the expression profiles of the CESA and GMPP genes were further analyzed, and these results indicated that TBP, NCBP2, and CYP71 were the top three most stable reference genes after comprehensive comparison, which could be used as stable reference genes for normalizing the genes expression in D. huoshanense. This study described here provides the first data regarding on reference gene selection in D. huoshanense, which will be extremely beneficial for future research on the gene expression normalization in D. huoshanense.


Introduction
D. huoshanense, as named by C.Z. Tang and S.J. Cheng, is a perennial epiphytic herb listed in Chinese Pharmacopoeia as having high edible and medicinal value; it belongs to the Dendrobium genus of Orchidaceae [1]. The wild D. huoshanense has been in an extremely endangered state, caused by the low natural reproduction rate and long-term destructive picking [2]. However, it has been greatly protected and industrialized through tissue culture and artificial cultivation technology [3]. As one of the most important traditional Chinese medicines, it is also known as "soft gold", and has been used as tonic for centuries [4,5]. Its stem is used as the medicinal part and it is generally considered that the plants with more juice and sticky teeth are the best [1]. For thousands of years, it was widely used as an ingredient in herbal medicines for nourishing yin and clearing heat, benefiting the stomach, and generating fluid [6]. The major active components of D. huoshanense are polysaccharides, flavonoids, sesquiterpenes, phenols, etc. [3,[7][8][9]. Many pharmacological experiments showed that it has the effect of enhancing immunity, anti-oxidation, antiaging, anti-tumor, and that it plays a role in reducing blood sugar, protecting the liver,

Plant Materials
The plants used in the experiment were potted and planted in the greenhouse of the medicinal botanical garden of West Anhui University, which were identified as D. huoshanense by Professor Bangxing Han from West Anhui University. One-year-old D. huoshanense plants were used as experimental subjects. Specifically, the roots, stems and leaves of fresh D. huoshanense were used as samples of different tissues. For oxidative stress, the plant leaves were treated with 50 mM hydrogen peroxide (H 2 O 2 ) once for 24 h. For drought treatment, 200 mL 25% PEG 6000 was used to treat plants every day for seven consecutive days. For temperature treatment, the plants were placed in a light incubator at 4 • C for two days. For UV treatment, the plants were exposed 15 cm from the UV light source (Philips TUV 30 W, 92 µW/cm 2 , 253 nm) for 15 min of radiation [47], and then a dark culture for 2 days. For hormone treatment, plant leaves were sprayed once with 25 mM MeJA and sampled after 6 h. All collected samples were quickly washed with distilled water and immediately frozen in liquid nitrogen and stored at −80 • C until RNA was extracted. Each group of samples had three biological replicates, and untreated D. huoshanense was used as control.

Total RNA Extraction and cDNA Synthesis
The total RNA of the cryopreservation samples from the refrigerator at −80 • C was extracted using an EASYspin Plus Complex Plant RNA kit (Aidlab, China), referring to the manufacturer's instructions, and DNA was eliminated with DNase I (Aidlab, China). The purity and concentration of RNA samples were detected using 1% agarose gel electrophoresis and a Synergy H1 multifunction microplate detector (BioTek, American), respectively. Additionally, 0.2 µg of total RNA with a 260/280 ratio between 1.9 and 2.1 was used to synthesize first-strand cDNA using HiScript ® II Q RT SuperMix for qPCR (Vazyme, Nanjing, China) in accordance with the manufacture's manual, and then the sample was diluted five times for qPCR analyses.

Selection of Candidate Reference Genes and Primers Design
Here, a total of 15 genes were selected as candidate reference genes by comparison with the TAIR database (http://www.arabidopsis.org (accessed on 7 February 2021)) using protein sequences of common housekeeping genes (GAPDH, eIF, EF-1α, PP2A, UBCE, RPL5, TBP, APT1, MDH, PTBP3, PEPC, CYP71, NCBP2, TIP41, and F-box) from Arabidopsis thaliana L. as templates. Then, the local BLAST program of BioEdit was employed to obtain the nucleotide sequences of putative D. huoshanense homologs based on the transcriptome data (access No. is PRJNA577972), and the information on 15 gene sequences of D. huoshanense was shown in Table S1. Next, Primer Premier 5 software was applied to design all the primers using the following principle: (1) the primer sequence length is 18-24 bp; (2) the amplification length is 80-300 bp; (3) the melting temperature (Tm) is 58-62 • C; (4) the GC content is 40-60%. Finally, all primers were synthesized by General Biol Company (Anhui, China) and checked by regular PCR products with 2% agarose gel electrophoresis before qPCR analysis. The melting curve was obtained to further determine primer specificity under the reaction condition at 95 • C for 15 s, 60 • C for 60 s, and 95 • C for 15 s using an AceQ qPCR SYBR Green Master Mix (Vazyme, Nanjing, China). The PCR efficiency (E) and correlation coefficient (R 2 ) of each gene was obtained directly from the StepOne TM Realtime PCR system (Applied Biosystems, Waltham, MA, USA) according to a standard curve generated from the 5-fold dilution cDNA series. All the information of the gene-specific primer pairs used in this study is listed in Table 1.  Figure S1.

qPCR Analysis
According to the introduction of the AceQ qPCR SYBR Green Master Mix, the qPCR reaction system was conducted in a 20 µL mixture with 10 µL AceQ SYBR Green Master Mix (High ROX Premixed), 0.2 µM of forward and reverse primers, 2 µL diluted cDNA template Genes 2022, 13, 1486 5 of 16 mentioned in item 2.2, and 7.2 µL of RNase-free ddH 2 O. The reaction condition followed the illustration from the manufacturer, which is 95 • C for 5 min for pre-denaturation, 40 cycles of 95 • C for 10 s for denaturation, and 60 • C for 30 s for annealing/extension. All experiments were performed with three independent biological and technical repetitions.

Gene Expression Stability Analysis
In order to evaluate the stability of each candidate gene under various experimental conditions, including different tissues (root, stem, and leaf), abiotic stresses (oxidation, drought, cold, and UV) and hormone treatment (MeJA), three different algorithms, geNorm, NormFinder, and BestKeeper, were employed to carry out statistical analysis. Specifically, geNorm calculated the expression stability values (M) and pairwise variation comparison (Vn/Vn + 1). NormFinder assessed the reliability of reference genes based on their variations in both the intra-group and inter-group. Different to geNorm and NormFinder, which needed the 2 −∆Ct value [43,44], BestKeeper directly used the Ct values to rank the gene expression stability by calculating CV ± SD (coefficient of variation ± standard deviation) to choose the most stable genes.

Comprehensive Analysis and Validation of Selected Reference Genes
Firstly, the website tool RefFinder (https://www.heartcure.com.au/reffinder/ (accessed on 1 August 2021)) was used to comprehensively reassess the results of the candidate genes stability analysis from the geNorm, NormFinder, and BestKeeper softwares. Then, in order to examine the reliability of the selected reference genes, the two most stable and one least stable internal reference genes were used for normalizing the expression levels of two target genes (GMPP and CESA) [20,46] related to polysaccharide biosynthesis in different tissues and hormone treatment. Sample collections and experiments were conducted as mentioned above. The expression levels of two target genes were normalized by the 2 −∆∆Ct method [48].

Primer Specificity Verification and PCR Efficiency
To investigate the reference gene of D. huoshanense, 15 candidate genes (GAPDH, eIF, EF-1α, PP2A, UBCE, RPL5, TBP, APT1, MDH, PTBP3, PEPC, CYP71, NCBP2, TIP41, and F-box) were chosen based on previous studies [36,37] and the TAIR database. The specific amplification of primer pairs of all candidate reference genes was firstly checked by regular PCR, which exhibited a single PCR band meet the expected size in 2.0% agarose gel electrophoresis analysis ( Figure 1). Moreover, the melting curve analysis showed a single amplification peak which further confirmed the specificity of the genes ( Figure 2). According to the standard curve of each internal reference gene obtained by diluting cDNA in a 5-fold gradient, the amplification efficiency ranged from 90.502% for PEPC to 104.77% for TBP, and all the correlation coefficients (R 2 ) were greater than 0.990 (Table 1).

Expression Profile of the Reference Genes
The raw Ct values of the reference genes were collected and are exhibited in Figure  3. It can be observed that the number of cycles correspond to the threshold of fluorescence

Expression Profile of the Reference Genes
The raw Ct values of the reference genes were collected and are exhibited in Figure 3. It can be observed that the number of cycles correspond to the threshold of fluorescence level which can be detected. The average Ct values of the 15 reference genes ranged from 22.00 to 28.58. According to the average Ct values, it can be judged that the expression abundance of 15 internal reference genes from high to low was Low Ct values correspond to high gene expression abundance, so the expression abundance of MDH was the highest and TIP41 was the lowest. Here, TBP, APT1, NCBP2, and MDH, have relatively narrow Ct variation ranges, indicating that their expression levels may be more stable. However, considering the complexity and diversity of the plant growth environment, the stability of the selected internal reference genes still needs to be further investigated under different environmental conditions.
Low Ct values correspond to high gene expression abundance, so the expression abundance of MDH was the highest and TIP41 was the lowest. Here, TBP, APT1, NCBP2, and MDH, have relatively narrow Ct variation ranges, indicating that their expression levels may be more stable. However, considering the complexity and diversity of the plant growth environment, the stability of the selected internal reference genes still needs to be further investigated under different environmental conditions.

The Analysis of Expression Stability of Candidate Reference Genes
In order to improve the accuracy of analysis, three frequently used algorithms, namely geNorm, NormFinder, and BestKeeper, were separately employed to further assess the stability of candidate reference genes under different treatments and tissues.

geNorm Analysis
In geNorm analysis, measurement (M) values of expression stability of each reference gene is generated for each pair of genes. This ranks the expression stability by M values, and the lower M value represents the higher stability of gene expression. Generally, if M > 1.5, it could be regarded as an unsuitable reference gene. As shown in Figure 4, the stability trend of the 15 candidate internal reference genes in the tissue group, MeJA group, and different abiotic stress groups were arranged from high to low by the geNorm analysis, and the most stable internal reference genes given by all groups had some differences. Concretely, the five most stable internal reference genes in each group were as follows: for the total group. More uniformly, for all samples, except for PTBP3 and F-box in the cold group, the two most unstable genes in other groups were TIP41 and GAPDH. In addition, geNorm can also recommend the required number of optimal internal reference genes according to the value of the pairwise variation (Vn/Vn + 1). When Vn/Vn + 1 is less than 0.15, the optimal number of reference genes is n [43]. According to this standard, the V2/V3 values are lower than 0.15 in all groups, so it is considered that

The Analysis of Expression Stability of Candidate Reference Genes
In order to improve the accuracy of analysis, three frequently used algorithms, namely geNorm, NormFinder, and BestKeeper, were separately employed to further assess the stability of candidate reference genes under different treatments and tissues.

geNorm Analysis
In geNorm analysis, measurement (M) values of expression stability of each reference gene is generated for each pair of genes. This ranks the expression stability by M values, and the lower M value represents the higher stability of gene expression. Generally, if M > 1.5, it could be regarded as an unsuitable reference gene. As shown in Figure 4, the stability trend of the 15 candidate internal reference genes in the tissue group, MeJA group, and different abiotic stress groups were arranged from high to low by the geNorm analysis, and the most stable internal reference genes given by all groups had some differences. Concretely, the five most stable internal reference genes in each group were as follows: for the total group. More uniformly, for all samples, except for PTBP3 and F-box in the cold group, the two most unstable genes in other groups were TIP41 and GAPDH. In addition, geNorm can also recommend the required number of optimal internal reference genes according to the value of the pairwise variation (Vn/Vn + 1). When Vn/Vn + 1 is less than 0.15, the optimal number of reference genes is n [43]. According to this standard, the V2/V3 values are lower than 0.15 in all groups, so it is considered that the optimal number of internal reference genes is 2 ( Figure 5), indicating that 2 reference genes can be sufficient for gene normalization. the optimal number of internal reference genes is 2 ( Figure 5), indicating that 2 reference genes can be sufficient for gene normalization.    the optimal number of internal reference genes is 2 ( Figure 5), indicating that 2 reference genes can be sufficient for gene normalization.

NormFinder Analysis
The NormFinder can directly evaluate the stability of reference genes, based on the variance in intra-and inter-group, to calculate the normalization factors using ANOVA. Generally, the lower stability value also reflects better stability of the corresponding ref-erence gene expression. As shown in Table 2 from the top to the bottom, the stability values of reference genes are listed from the lowest to the highest with the NormFinder analysis. In detail, the most stable internal reference genes could be discovered clearly in all experimental conditions. The five most stable internal reference genes in each group were as follows: TBP > NCBP2 > eIF > RPL5 > PTBP3 for the H 2 O 2 group; NCBP2 > CYP71 > TBP > PP2A > eIF for the PEG group; TBP > CYP71 > PP2A > APT1 > UBCE for the cold group; EF-1α > CYP71 > NCBP2 > TBP > UBCE for the UV group; TBP > APT1 > PP2A > EF-1α > eIF for the MeJA group; MDH > NCBP2 > RPL5 > UBCE > TBP for the tissue group; NCBP2 > MDH > APT1 > UBCE > TBP for the total group. According to the comprehensive analysis, TBP can be regarded as a most stable reference gene, similar to the results with the geNorm analysis. In addition to the unstable expression of NCBP2 in the cold and MeJA group, it was also very stable in other groups. Furthermore, the two most unstable genes (GAPDH and TIP41) in each group were consistent with the results of the geNorm analysis apart from the cold group. Combined with the analysis results of geNorm, we can judge the best combinations of the two internal reference genes in each group as follows: TBP + NCBP2 for the H 2 O 2 group; NCBP2 + CYP71 for the PEG group; TBP + CYP71 for the cold group; CYP71 + EF-1α for the UV group; TBP + APT1 for the MeJA group; RPL5 + MDH for the tissue group; APT1 + NCBP2 for the total group.

BestKeeper Analysis
BestKeeper is an Excel-based program that directly uses the raw Ct values without conversion to calculate coefficient of variation (CV) and standard deviation (SD) of the reference gene from each experimental group, so as to compare their stability [45]. The smaller CV ± SD value, the more stable the reference gene. If the SD value is greater than 1, it is generally considered that this gene is unstable and should be discarded. All analysis results are shown in Table 3, and the stable arrangement of each gene is somewhat different from that of geNorm and NormFinder. Specifically, the five most stable internal reference genes in each group were as follows: TBP > RPL5 > NCBP2 > eIF > PTBP3 for the H 2 O 2 group; NCBP2 > CYP71 > TBP > eIF > PP2A for the PEG group; EF-1α > UBCE > CYP71 > APT1 > TBP for the cold group; EF-1α > CYP71 > TBP > UBCE > NCBP2 for the UV group; eIF > EF-1α > TBP > APT1 > NCBP2 for the MeJA group; PTBP3 > PEPC > UBCE > TBP > EF-1α for the tissue group; MDH > TBP > PTBP3 > UBCE > CYP71 for the total group. However, for most treatments, TBP or NCBP2 still showed good stability, basically consistent with the geNorm and NormFinder analyses. In addition, except for the Cold group, GAPDH and TIP41 were still more unstable than others, a result which was nearly close to the results of geNorm and NormFinder analysis.

Comprehensive Analysis and Validation of Reference Genes
Considering the differences in stability analysis from the first three programs, we further used the online comprehensive ranking platform RefFinder, which used geNorm, Normfinder, BestKeeper, and the comparative ∆Ct method to verify the rankings of candidate reference genes [49] (Table 4). The geometric mean of the attributed weight from the stable value of each gene was calculated by RefFinder, and the smaller the value, the more stable the internal parameter gene. The ranking order of the top five most stable and unstable candidate reference genes acquired by RefFinder were basically the same as the results provided by geNorm and Normfinder, which is slightly different from the results of BestKeeper. For example, in all treatment groups, the first two most unstable genes generated by RefFinder were basically the same as the first three most unstable genes calculated by geNorm, Normfinder, and BestKeeper. Moreover, the stability ranking order of the 15 reference genes was counted in Table S2 to better observe the rankings of the 4 software analysis results, which suggested that TBP, NCBP2, and CYP71 could be the 3 best stable genes in D. huoshanense under the different conditions.  Then, in order to further verify the reliability of the screened stable internal reference genes, the expression of two genes (CESA and GMPP) related to polysaccharide biosynthesis pathway in different tissues (root, stem, and leaf) and MeJA treatment for stem were analyzed by qPCR with the unstable reference gene (TIP41), three stable reference genes (RPL5, MDH, and TBP), and their combinations as reference genes based on the comprehensive rankings mentioned in Table S2. The results showed that, when RPL5, MDH and RPL5 + MDH were used as internal reference genes, the expression levels of CESA and GMPP were slightly different in the root, stem, and leaf samples. Among them, the expression levels of CESA were stem > root > leaf, and the expression levels of GMPP were stem > leaf > root ( Figure 6A). However, when the unstable internal reference gene TIP41 was used, the relative expression levels of CESA and GMPP in root and stem were significantly different from the quantitative results of RPL5, MDH, and RPL5 + MDH ( Figure 6B). Likewise, when the stable APT1 or TBP was selected as internal reference gene, it was obvious that the relative quantitative results of CESA or GMPP in the stem were similar, and while using the unstable internal reference gene TIP41, the relative quantitative results of CESA or GMPP in stem were significantly different ( Figure 6C). generated by RefFinder were basically the same as the first three most unstable genes calculated by geNorm, Normfinder, and BestKeeper. Moreover, the stability ranking order of the 15 reference genes was counted in Table S2 to better observe the rankings of the 4 software analysis results, which suggested that TBP, NCBP2, and CYP71 could be the 3 best stable genes in D. huoshanense under the different conditions. Then, in order to further verify the reliability of the screened stable internal reference genes, the expression of two genes (CESA and GMPP) related to polysaccharide biosynthesis pathway in different tissues (root, stem, and leaf) and MeJA treatment for stem were analyzed by qPCR with the unstable reference gene (TIP41), three stable reference genes (RPL5, MDH, and TBP), and their combinations as reference genes based on the comprehensive rankings mentioned in Table S2. The results showed that, when RPL5, MDH and RPL5 + MDH were used as internal reference genes, the expression levels of CESA and GMPP were slightly different in the root, stem, and leaf samples. Among them, the expression levels of CESA were stem > root > leaf, and the expression levels of GMPP were stem > leaf > root ( Figure 6A). However, when the unstable internal reference gene TIP41 was used, the relative expression levels of CESA and GMPP in root and stem were significantly different from the quantitative results of RPL5, MDH, and RPL5 + MDH ( Figure  6B). Likewise, when the stable APT1 or TBP was selected as internal reference gene, it was obvious that the relative quantitative results of CESA or GMPP in the stem were similar, and while using the unstable internal reference gene TIP41, the relative quantitative results of CESA or GMPP in stem were significantly different ( Figure 6C).

Discussion
D. huoshanense has attracted many researchers because of its broad pharmacological activity [9]. Our research group on D. huoshanense has carried out a series of studies on the field resource protection, tissue culture and cultivation, chemical component separation and identification, bioactivity analysis and action mechanism, and new drug and product development [3,5,10,39,[50][51][52][53], and also completed its whole genome sequencing [54]. In addition, the transcriptome sequencing of D. huoshanense had been reported several times in other research groups [13][14][15][16][17]. Since the main active components of similar plants have been gradually known, many studies also mainly focused on exploring and understanding the biological control [55] and biosynthetic pathways of polysaccharides, alkaloids and flavonoids, etc. [13][14][15][16][17].
The gene expression level has an inevitable and important role in the biosynthesis of secondary metabolites and related gene function mining. Due to its high sensitivity and specificity, qPCR is often used for high-throughput analysis at the level of gene transcription. To obtain accurate and reliable results for the data, a stable and suitable reference gene is required. In fact, in the quantitative studies of the genetic expression of D. huoshanense or its related species, multiple genes were randomly used as internal reference genes [14][15][16][17][18][19], which would cause errors in the analysis results. To our knowledge, there are still no systematical studies on reference genes in D. huoshanense. Moreover, there are three main cultivation methods of D. huoshanense including facility cultivation, under forest cultivation, and simulative habitat cultivation [3]. Among these, the simulative habitat cultivation of D. huoshanense has the comprehensive effects of "excellent environment", "excellent shape", "high quality", and "excellent effect" [3], which is similar to the report of Dendrobium officinale [56]. This cultivation method has almost no manual intervention in the growth process, and is mainly faced with various environmental stress factors, such as UV, temperature, and drought. In order to further study the molecular mechanism of the effects of environmental factors on the quality of D. huoshanense, it is necessary to screen the internal reference genes that can be stably expressed in different tissues, and under the main abiotic stress conditions (oxidation, drought, temperature, and ultraviolet radiation) and hormone stress. Thus, we set up six treatments in this study, and the results also confirmed the necessity of this configuration. Therefore, we comprehensively analyzed the expression levels and stability of 15 candidate reference genes under corresponding conditions.
All raw Ct values of D. huoshanense samples under different abiotic conditions, hormone treatment, and in different tissues were acquired from qPCR and processed using geNorm, NormFinder, and BestKeeper for ranking the reference genes. According to the results of this study, the 15 candidate reference genes firstly exhibited a relatively reliable range for expression profiles from 22.00 to 28.58, revealing that the selected candidate genes were able to offer an accurate normalization. Based on the relatively narrow Ct variation ranges, TBP, APT1, NCBP2, and MDH might be tentatively considered as the better stable reference genes ( Figure 3). Obviously, this results on the rankings of 15 candidate reference genes were somewhat different from the outcomes calculated by geNorm, NormFinder, and BestKeeper, which indicated that it is necessary to use multiple procedures to obtain the best results (Figure 4, Tables 2 and 3), and also revealed that none of the 15 reference genes could be stably expressed for meeting all conditions in D. huoshanense. Next, three Excel-based programs were further employed to estimate the stability of the candidate genes, and their results showed some differences in ranking order because their analysis principle and emphasis are different. In the analysis of geNorm and NormFinder, the five most stable genes given by them were basically the same, while the difference was their different stability ranking order. However, the analysis results of BestKeeper are different from those of geNorm and NormFinder, mainly due to the CV and SD values, which are the key factors for determining the stability ranking of reference genes obtained by BestKeeper. This discrepancy was acceptable because the results of the three methods showed some consistency in selecting the five most stable genes. Eventually, the above three analysis results were synthesized by RefFinder, and the five most stable genes were close to those given by geNorm and NormFinder. Combining the results of these four software analysis (Table S2), the five most stable genes were as follows: TBP > NCBP2 > RPL5 > eIF > PTBP3 in the H 2 O 2 group; NCBP2 > CYP71 > TBP > PP2A > eIF in the PEG group; TBP > CYP71 > PP2A > UBCE > APT1 in the cold group; EF-1α > CYP71 > TBP > NCBP2 > UBCE in the UV group; TBP > APT1 > eIF > EF-1α > PP2A in MeJA group; TBP > UBCE > MDH > RPL5 > NCBP2 in the tissue group; MDH > TBP > UBCE > NCBP2 > APT1 in the total group. Accordingly, these results demonstrated that TBP, NCBP2, and CYP71 could be comprehensively regarded as the top three stable reference genes in D. huoshanense. Unexpectedly, although GAPDH is a commonly used reference gene, its expression in D. huoshanense showed fairly poor expression stability under all experimental conditions in our research, which meant that it could not be used as a stable reference gene for qPCR analysis. However, it was regarded as the best reference gene in D. catenatum [57], Diabrotica undecimpunctata [58], and Corydalis yanhusuo [59]. Similarly, although TIP41 is not suitable as an internal reference gene in this study, its expression presented high stability in Momordica charantia [60]. These results further illustrated that there is no universally applicable reference gene with invariant expression, and that the selection of stable reference genes remain indispensable. Furthermore, in qPCR analysis, multiple reference genes are better than a single reference gene for normalization. According to the pairwise variation results (Figure 5), the combination of two reference genes should be enough to meet the demand for the normalization in different tissues and under all the experimental conditions ( Figure 6).
In order to further confirm the suitability of the selected reference genes, two pairs of the most stable genes (RPL5 and MDH, and APT1 and TBP) and their combinations and one least stable gene (TIP41) were selected for the normalization of CESA [46] and GMPP [20] involved in the polysaccharide biosynthesis of D. huoshanense under MeJA treatments and in various tissues. Noticeably, when the selected stable reference gene pairs or their combinations were standardized alone, there were only slight differences in the relative expression levels of CESA and GMPP between different tissues and MeJA treatment. However, when using the least stable reference gene, a significantly different result was found in the relative expression levels of CESA and GMPP between different tissues and MeJA treatments, indicating that the selection and confirmation of appropriate and stable reference genes is particularly critical before they are used for a set of samples.

Conclusions
To the best of our knowledge, our study is the first to systematically explore and evaluate the expression stability of 15 reference genes from D. huoshanense under different abiotic stress factors, hormone treatment, and in different tissues. The results indicated that TBP, NCBP2, and CYP71 could be regarded as the optimal reference genes based on their better stability in different conditions when analyzed by four commonly used programs (geNorm, NormFinder, BestKeeper, and RefFinder). Furthermore, the two unstable genes of each group were basically identical according to the comprehensive ranking order from the four programs at their relevant conditions. Finally, the validation experiments on the expression analysis of CESA and GMPP further accentuated that it is necessary to screen stable reference genes for the normalization of gene expression analysis by qPCR under different conditions. This study will be useful to increase the accuracy of gene expression analysis by qPCR and promote future research on gene functions in D. huoshanense and related Dendrobium species. Therefore, this research should also arouse great interest in researchers engaged in the mining of key genes in plant secondary metabolic pathways, inspiring further effort in plant molecular research and natural product biosynthesis pathway analysis.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/genes13081486/s1, Table S1: Information of 15 candidate reference genes and two target genes selected for evaluation and validation in D. huoshanense. Table S2: Expression stability rankings of the 15 candidate reference genes estimated by geNorm, NormFinder, BestKeeper and RefFinder. Figure S1: Standard curves of 15 candidate reference genes and two target genes were directly generated by StepOne TM Real-time PCR system.