Identification of Reference Genes for RT-qPCR Analysis in Gleditsia microphylla under Abiotic Stress and Hormone Treatment

Gleditsia microphylla is an important galactomannan gums source plant with characteristics of drought resistance, barren tolerance, and good adaptability. However, the underlying molecular mechanisms of the biological process are not yet fully understood. Real-time quantitative PCR (RT-qPCR) is an accurate and convenient method to quantify the gene expression level and transcription abundance of suitable reference genes. This study aimed to screen the best internal reference genes in G. microphylla under abiotic stresses, hormone treatments, and different tissues. Based on the transcriptome data, twelve candidate reference genes were selected, and ultimately, nine of them were further evaluated by the geNorm, NormFinder, BestKeeper, and RefFinder algorithms. These results show that TATA-binding protein 1 (TBP1)and Eukaryotic translation initiation factor 4A1 (EIF4A1)were the two most stable reference genes, and glyceraldehyde-3-phosphate dehydrogenase A subunit, chloroplastic (GAPA)and glyceraldehyde-3-phosphate dehydrogenase B subunit, chloroplastic (GAPB)were the two most unstable reference genes across all samples under the given experimental conditions. Meanwhile, the most stable reference genes varied among the different groups and tissues. Therefore, this study suggests that it is better to use a specific reference gene for a particular case rather than using a common reference gene.


Introduction
Real-time quantitative PCR (RT-qPCR) is a technique that can precisely quantify nucleic acid molecules by monitoring fluorescence signals during the entire PCR process. It has been widely utilized in gene expression and transcript abundance analysis due to its high sensitivity, good repeatability, and strong specificity [1][2][3][4][5]. Although RT-qPCR is a powerful tool for understanding gene roles in metabolic pathways, signaling pathways, and complex regulatory networks in organisms, the normalization accuracy depends on suitable and stable reference genes as internal standards [6][7][8][9]. For an ideal reference gene used for normalization in real-time PCR analysis, its expression should remain constant between the cells of different tissues and under various experimental conditions. Commonly, the housekeeping genes involved in fundamental cellular processes such as 18S rRNA, ACT (actin-related protein), TUB (β-1 tubulin), and GAPDH (glyceraldehyde-3-phosphate dehydrogenase) are used as reference genes in RT-qPCR for plants [2][3][4][5][6][7][8][9]. However, no universal reference genes with stable expression profiles in different tissues and organs, developmental stages, and experimental conditions have been discovered [10][11][12][13][14]. Therefore, it is vital to identify suitable internal reference genes to study the expression levels of target genes in different sample types, and under various experimental conditions. Several statistical algorithms such as geNorm [15,16], NormFinder [17], BestKeeper [18], and RefFinder [19] have (SA), while the control seedlings received only water. For heavy-metal and salinity stresses, the seedlings were treated with Hoagland solution supplemented with 200 µM copper sulfate (CuSO 4 ) or 100 µM sodium chloride (NaCl). Samples of the top third and fourth leaves, stem, and root were collected 12 h after hormone, heavy-metal, and salt treatments. For cold and heat shock treatments, the plants were grown in a growth chamber at 4 • C or 42 • C for 24 h. For simulated drought treatment, the seedlings were subjected to Hoagland nutrient solution with 10% PEG6000 for 7 d. The mock-treated seedlings with the same time interval served as the control. Male and female flower samples were collected from adult plants in Jinzhou, Hubei Province. All the samples were frozen in liquid nitrogen and preserved at −80 • C.

Extraction of Total RNA and cDNA Synthesis
Total RNA was isolated from the frozen tissues using the EasyPure ® Plant RNA Kit (TransGen Biotech, Beijing, China) and quantified by a NanoDrop 2000c spectrophotometer (Thermo Scientific, Waltham, MA, USA). Only RNA samples with A 260 /A 280 of 1.9-2.1 were used for cDNA synthesis. RNA integrity was analyzed by 1% agarose gel electrophoresis. The first-strand cDNA was synthesized with 1.0 µg total RNA in a 20 µL reaction system according to the TransScript ® All-in-One First-Strand cDNA Synthesis SuperMix for qPCR (One-Step gDNA Removal) Kit (TransGen Biotech, Beijing, China). All the cDNA samples were stored at −20 • C.

Real-Time PCR
A CFX96 Touch Real-Time PCR Detection System (Bio-Rad Laboratory, Inc., Hercules, CA, USA.) was used, and the program (3 min at 95 • C followed by 40 cycles at 94 • C for 15 s, and at 60.5 • C for 60 s) employed for RT-qPCR used a reaction mixture volume of 20 µL in an optical 96-well plate. Then, 10.0 µL of AceQ ® SYBR Green Master Mix (Vazyme Biotech, Nanjing, China), 0.3 µL of each final primer (150 nM), 2.0 µL of final cDNA (20 ng), and 7.4 µL RNase-free water were added to the reaction mixture. A control was also included in each plate with 2.0 µL of RNase-free water as a template. Two or three technical replicates were included for each biological replicate contained in each plate. The threshold cycle (Ct) values were calculated automatically by CFX Manager according to the overall expression levels of each gene analyzed.

Statistical Analysis of Gene Expression and Comparison of Normalization Methods
The Ct values obtained from RT-qPCR for each candidate reference gene were inputted into software or an online website and analyzed according to the corresponding manuals of geNorm [15,16], NormFinder [17], BestKeeper [18], and RefFinder [19]. Briefly, when analyzed by geNorm and NormFinder, raw Ct values of each candidate must be converted to relative quantitative values (2 −∆Ct ), and an expression stability measurement (M) value and pairwise variation (V) value must be calculated by geNorm. All candidates were ranked by the M values, where the smaller the M value, the better the stability. The V value was used to determine the optimal numbers of reference genes, and it is generally considered that when the value of Vn/n + 1 is more than 0.15, the (n + 1) th reference gene is in need. Compared with geNorm, NormFinder was used to compare the expression differences in candidate reference genes based on the calculated stability value. Similar to geNorm, the smaller the M value, the better the stability of the candidate; however, the NormFinder program can only select the most suitable gene as the internal reference gene. BestKeeper directly utilized the raw Ct value for stability analysis by calculating the coefficient of variance (CV) and the standard deviation (SD); the criteria for gene stability were smaller CV and SD values. Finally, the web-based tool RefFinder (http://blooge.cn/RefFinder/) (accessed on 1 May 2022) integrated all three methods mentioned above and raw Ct values to calculate the geometric mean for each reference gene and the comprehensive ranking index of stability. A lower index value indicates a higher stability of the candidate.

Validation of Reference Genes by Expression Analysis of Stress-Responding Genes
Under stress conditions, ethylene, ABA, and reactive oxygen species (ROS) are essential signal transduction molecules to activate the defense system in plants, and ACO (ethylene-forming enzyme) [28], PYLs (abscisic acid receptors) [29,30], and CSD2 (superoxide radicals detoxifying enzyme) [31,32] play key roles during the process. To validate the selected reference genes in this study, the expression levels of these three homologs were analyzed using the most and least stable reference genes under stress conditions and calculated using the 2 −∆∆Ct method [33]. Three biological replicates were included for each treatment, and three technical replicates were included for each biological sample.

Specificity of Primers for Candidate Reference Genes and Target Genes
To detect the specificity of the designed primers for the twelve candidate reference genes and three target genes, analyses of the gel electrophoresis of the PCR products and melting curves were performed. The gel electrophoresis showed a single band with the expected size of each primer pair (Figure 1), and the melting curves of each primer pair exhibited a single peak ( Figure S1), indicating the specificity of these primer pairs of candidate genes and target genes. The standard curves indicated that the RT-qPCR amplification efficiency of the candidate reference genes ranged from 86.90% (TIP41) to 118.60% (CBP20), and the correlation coefficients (R 2 ) varied from 0.982 (CBP20) to 1.000 (ACT7) (Table S1, Figure S2). Thus, all the primer pairs were specific for their respective genes and could be used in RT-qPCR analysis except TIP41 and ACT1 because of their lower amplification efficiencies. smaller CV and SD values. Finally, the web-based tool RefFinder (http://blooge.cn/RefFinder/) ( accessed on 1 May, 2022) integrated all three methods mentioned above and raw Ct values to calculate the geometric mean for each reference gene and the comprehensive ranking index of stability. A lower index value indicates a higher stability of the candidate.

Validation of Reference Genes by Expression Analysis of Stress-Responding Genes
Under stress conditions, ethylene, ABA, and reactive oxygen species (ROS) are essential signal transduction molecules to activate the defense system in plants, and ACO (ethylene-forming enzyme) [28], PYLs (abscisic acid receptors) [29,30], and CSD2 (superoxide radicals detoxifying enzyme) [31,32] play key roles during the process. To validate the selected reference genes in this study, the expression levels of these three homologs were analyzed using the most and least stable reference genes under stress conditions and calculated using the 2 −ΔΔCt method [33]. Three biological replicates were included for each treatment, and three technical replicates were included for each biological sample.

Specificity of Primers for Candidate Reference Genes and Target Genes
To detect the specificity of the designed primers for the twelve candidate reference genes and three target genes, analyses of the gel electrophoresis of the PCR products and melting curves were performed. The gel electrophoresis showed a single band with the expected size of each primer pair (Figure 1), and the melting curves of each primer pair exhibited a single peak ( Figure S1), indicating the specificity of these primer pairs of candidate genes and target genes. The standard curves indicated that the RT-qPCR amplification efficiency of the candidate reference genes ranged from 86.90% (TIP41) to 118.60% (CBP20), and the correlation coefficients (R 2 ) varied from 0.982 (CBP20) to 1.000 (ACT7) (Table S1, Figure S2). Thus, all the primer pairs were specific for their respective genes and could be used in RT-qPCR analysis except TIP41 and ACT1 because of their lower amplification efficiencies.

Expression Profiling of Candidate Reference Genes
Analyses of the expression levels of the remaining ten candidates were performed in all samples using RT-qPCR. The statistical results show that PTB1 had the highest average Ct value of 32.08 among all the candidates, and the maximum values were over 35.00 in some root samples, which implied it was expressed at a low level and unsuitable to act as a normalization gene. Therefore, the remaining nine candidates were further analyzed, and the raw Ct values obtained using RT-qPCR in various samples are shown in a boxplot ( Figure 2). In all samples, the Ct values varied from 21.08 (GAPA in leaves spritzed with ABA) to 33.00 (TBP1 in roots stressed by NaCl), and the average Ct values ranged from 25.89 (GAPA) to 30.76 (CBP20), indicating that these candidate genes present different expression levels under experimental conditions. Intra-and intergroup statistical analyses of Ct values of nine candidates were further conducted, the results demonstrated that significant differences (p < 0.05) were always observed, similar to the results in organs. For example, in the cold group, Ct values of TBP1 showed significant difference (p < 0.001) from that of GAPA and eIF2, but no significant difference from the other candidates. Among groups, Ct values of TBP1 in the ABA group had significant difference (p < 0.001) from that of the cold and heat groups. In organs, Ct values of TBP1 in roots showed significant difference (p < 0.001) from that of the stems and leaves. However, although GAPA exhibited a relatively higher expression, expression bias was evident in various organs as the average Ct values were 30.46, 24.65, and 22.90 in the total roots, stems, and leaves, respectively. Moreover, TBP1, EF-1α, CBP20, and EIF4A1 had a relatively narrower Ct value range than the other genes, implying these genes might be expressed more stably.

Expression Stability of Candidate Reference Genes
To evaluate the stability of the candidate reference genes, the raw Ct values obtained from all samples were analyzed by geNorm, NormFinder, and BestKeeper; and a comprehensive stability ranking was finally generated by RefFinder. In the geNorm analysis performed for each group, all the groups showed a V2/3 value that was more significant than the threshold value of 0.15 ( Figure S4), which suggests that it was difficult to find common reference genes for each treatment. It is better to identify suitable reference genes for individual tissues under a specific condition. Hence, we performed expression stability analysis on two sets of data: (1) a group set that combined the data of all organs under a specific condition; (2) data from individual organs under a specific condition. An investigation of all the stresses and stimuli was also performed.

geNorm Analysis
In geNorm analysis, a cut-off M value of 1.5 is recommended for evaluating all genes' stability. In this study, the M values of the candidate reference genes were all lower than 1.5 except GAPA in all groups and EF-1α in the PEG group when analyzed by the group set ( Figure 3). Still, they were all lower than 1.5 when analyzed by individual organs in each group ( Figure S3). As shown in Figure 3, the two most stable reference genes were not identical under different treatments; they were TBP1 + eIF2 for the cold and heat groups, TBP1 + EIF4A1 for the PEG and tissue groups, TBP1 + CBP20 for the NaCl and total groups, EF-1α + TUB1 for the CuSO 4 group, EF-1α + EIF4A1 for the ABA group, EIF4A1 + TUB1 for the MeJA group, and CBP20 + EIF4A1 for the SA group, according to their lowest M values. The most stable reference genes varied more when analyzed by organ; they were not identical among organs even under the same condition ( Figure S3). Interestingly, GAPA was the most unstable gene in all groups, but it was the most stable gene in the NaCl-and MeJA-treated roots, heat-, PEG-, and SA-treated and CK stems, and CK leaves, which suggests that it can be used as a reference gene under specific conditions. Moreover, the appropriate number of reference genes for normalization was investigated using a threshold value of 0.15 (Vn/n + 1). In our study, the V2/3 values in most organs under a specific condition were lower than 0.15, which meant two reference genes met the requirements for normalization. However, two reference genes were insufficient in CK roots, NaCl-treated stems, and heat-, CuSO 4 -, and ABA-treated leaves ( Figure 4).    . The pairwise variation (V) measures of the nine candidate reference genes using geNorm. Vn/n + 1 values were used to calculate the optimal number of reference genes (n). The letters R, S, L, and F represent root, stem, leaf, and flower, respectively. Vn/n + 1 values were used to calculate the optimal number of reference genes (n). The letters R, S, L, and F represent root, stem, leaf, and flower, respectively.

NormFinder Analysis
Similar to geNorm, the raw Ct values were transformed to relative quantitative values prior to NormFinder analysis. As a result, a variation value of each candidate was given, and all genes were ranked [17]. In this study, the most stable genes could be found directly in each organ under different experimental conditions, and the best combination of two genes was recommended when analyzed by the group set (Table 1). In most cases, the most stable gene varied among organs even under the same treatment, but occasionally it was the same. For example, the most stable gene was TBP1, EF-1α, and TUB1 in the roots, stems, and leaves under cold stress, respectively. However, for the ABA and MeJA stimuli, the most stable gene was the same in the roots, stems, and leaves. As for the best combination of two genes in a specific group, they were CBP20 + eIF2 for the total group, ACT7 + eIF2 for the PEG group, GAPB + TUB1 for the CuSO 4 group, TBP1 + eIF2 for the cold, heat, and ABA groups, and EIF4A1 + eIF2 for the NaCl, MeJA, SA, and tissue groups. Total

BestKeeper Analysis
Unlike geNorm and NormFinder, the raw Ct data could be used directly by BestKeeper. Finally, standard deviations (SDs) and coefficients of variation (CVs) between pairs of genes were obtained to evaluate the stability of reference genes in each experimental group. The most stable reference gene was considered to have the lowest CV and SD, and the SD values should be less than 1.0. As shown in Table 2, CV ± SD values were calculated and ranked by BestKeeper, and the stability of the reference genes decreased gradually from left to right in the table, similar to NormFinder. In the group's total and tissue groups, there were three genes with an SD value under the cut-off value of 1.0 in each group, but in the other groups, there were three to six genes. Three genes were in the cold group (TUB1, EIF4A1, and EF-1α), four in the heat group (EIF4A1, CBP20, TUB1, and TBP1), four in the PEG group (EF-1α, EIF4A1, TBP1, and ACT7), five in the CuSO 4 group (TUB1, EF-1α, CBP20, TBP1, and EIF4A1), four in the NaCl group (CBP20, TBP1, EF-1α, and EIF4A1), six in the ABA group (TBP1, TUB1, ACT7, EF-1α, CBP20, and EIF4A1), six in the MeJA group (TUB1, EF-1α, EIF4A1, ACT7, CBP20, and TBP1), and five in the SA group (ACT7, TUB1, TBP1, EF-1α, and EIF4A1). Due to the different analytical principles, the most stable reference gene evaluated by BestKeeper differed somewhat from the results of geNorm and NormFinder. However, the top three most stable reference genes in each group from BestKeeper also had some relatively consistent rankings with geNorm and Normfinder; in particular, TBP1, CBP20, and EIF4A1 appeared with the highest frequencies. Unexpectedly, the top two most unstable genes analyzed by the three software programs were the same in all group sets, namely, GAPA and GAPB. As for the results analyzed with BestKeeper by individual organs in each group, 55.56% (5/9) to 100.00% (9/9) of the candidate reference genes in each organ had an SD value below 1.0.

RefFinder Analysis
Due to the discrepancy in the stability of the candidate reference genes calculated by geNorm, Normfinder, and BestKeeper, the RefFinder program was applied to integrate these results to obtain a comprehensive ranking ( Table 3). The ranking order of the top three most stable and unstable reference genes obtained by RefFinder was broadly consistent with the results provided by geNorm and NormFinder but had slight differences with the results from BestKeeper. For instance, under all given experimental conditions, the top two most unstable genes generated by RefFinder appeared in the top three least stable genes yielded by geNorm, Normfinder, and BestKeeper. Ranking orders of the stability of nine reference genes were summarized to better observe the rankings of the analysis results from the four software programs (Table S3). * RefFinder exploits computational programs (such as BestKeeper, geNorm, NormFinder, or the comparative delta-Ct method) to rank and compare candidate reference genes. The values following the genes indicate the geometric mean of the attributed weights measured by this software for the overall final ranking.

Validation of Selected Reference Genes
In this study, PYL1 exhibited an average Ct value of 34.03 in all samples, which meant a relatively lower expression abundance compared with ACO1 (33.00) and CSD2 (26.34). To obtain a more reliable validation result, only the expression analysis of ACO1 and CSD2 was conducted. As depicted in Figure 5a, the expression levels of ACO1 and CSD2 in the roots and leaves under cold stress were overestimated when the least stable candidate gene was used as the internal control instead of the most stable gene. In contrast, their expression levels in the stems were underestimated, and the normalization differences between the most stable and unstable genes were significant. Moreover, in ABA-treated roots, ACO1 was enhanced at a similar level (no significant difference) when the two most stable reference genes (EF-1α and ACT7) and their combination (EF-1α + ACT7) were used as normalization factors, while a significant difference (p < 0.001) was observed when the least stable gene GAPA was used (Figure 5b). Therefore, incorrect results would be obtained when using an improper reference gene for the normalization of target genes, highlighting the importance of validating reference genes before conducting an RT-qPCR experiment, to obtain accurate results.

Validation of Selected Reference Genes
In this study, PYL1 exhibited an average Ct value of 34.03 in all samples, which meant a relatively lower expression abundance compared with ACO1 (33.00) and CSD2 (26.34). To obtain a more reliable validation result, only the expression analysis of ACO1 and CSD2 was conducted. As depicted in Figure 5a, the expression levels of ACO1 and CSD2 in the roots and leaves under cold stress were overestimated when the least stable candidate gene was used as the internal control instead of the most stable gene. In contrast, their expression levels in the stems were underestimated, and the normalization differences between the most stable and unstable genes were significant. Moreover, in ABAtreated roots, ACO1 was enhanced at a similar level (no significant difference) when the two most stable reference genes (EF-1α and ACT7) and their combination (EF-1α + ACT7) were used as normalization factors, while a significant difference (p < 0.001) was observed when the least stable gene GAPA was used (Figure 5b). Therefore, incorrect results would be obtained when using an improper reference gene for the normalization of target genes, highlighting the importance of validating reference genes before conducting an RT-qPCR experiment, to obtain accurate results.
(a) (b) Figure 5. Validation of the selected reference genes. Relative expression levels of ACO1 and CSD2 were normalized using candidate reference genes under different treatments. (a) The expression level was normalized using the most stable and unstable reference genes in various organs under cold treatment. The most stable and unstable reference genes in the roots, stems, and leaves were TBP1 and EF-1α, EF-1α and ACT7, and TUB1 and EF-1α, respectively. (b) The expression level of ACO1 was normalized using the most stable reference genes and their combination in roots under ABA treatment. EF-1α and ACT7 represent the two most stable reference genes, and GAPA is the most unstable gene. Data are displayed as means ± SEM (n = 3), and the statistical analyses were performed using multiple comparisons of one-way ANOVA to compare those between two reference genes or combinations for normalization. * p < 0.05; ** p < 0.01; *** p < 0.001, **** p < 0.0001. N.S.: no significant difference.

Discussion
G. microphylla is a member of the Fabaceae family. Its seeds are rich in galactomannan and proteins, which are broadly used in industries and animal feed. In addition, the plant has biological characteristics such as drought resistance, barren tolerance, and good adaptability which endows it with good ecological functions. Previous studies [23][24][25][26][27] mainly focused on the application of seeds and the investigation of biological traits, but not on the gene expression and underlying molecular mechanisms of biological processes. Therefore, de novo transcriptome sequencing of G. microphylla was conducted by our team for studies concerning sex determination, seed development, stress tolerance, and related gene functions. To achieve these objectives, stable and suitable internal reference genes for RT-qPCR analysis should be identified and evaluated first.
An ideal reference gene should be stably expressed regardless of the experimental condition and cell type. However, no universal reference genes for normalization using RT-qPCR have been discovered, making it particularly important to find proper reference genes when working with tissues of different histological origin or under different conditions. Owning to the absence of reports on reference genes in G. microphylla, twelve housekeeping genes (TBP1, CBP20, TIP41, EF-1α, EIF4A1, PTB1, ACT1, ACT7, GAPA, GAPB, TUB1, and eIF2) were selected from our transcriptome datasets for evaluation according to the literature [20,[34][35][36][37][38][39]. Among the designed primer pairs for all the candidate reference genes, ten presented high specificity and PCR efficiency, indicating that they could be used in RT-qPCR analysis (Table S1); and two, TIP41 (86.90%) and ACT1 (87.90%), were not further analyzed until higher-efficiency primer pairs were redesigned. In our study, non-coding genes such as rRNAs or miRNAs were not taken into consideration for controversial reasons, because some researchers thought that their high abundance compared with target mRNA transcripts made it difficult to normalize the expression of genes with low expression levels [5,13]. However, other researchers proved they were the best reference genes [7,40,41]. The remaining candidates showed a relatively wide range of expression profiles under the given experimental conditions, confirming again that no single gene could be used for normalization in all the samples of different tissues and under various experimental conditions, similar to the results in Osmanthus fragrans [42], Urochloa brizantha [43], Lycoris aurea [44], and Eleusine coracana [45].
The raw Ct values obtained using RT-qPCR in all samples under different conditions were the original data analyzed by the four most popular computational programs (geNorm, NormFinder, BestKeeper, and RefFinder) to rank the reference genes. Therefore, it is vital to ensure an optimal Ct value range for each reference gene through moderate dilution of cDNA. Compared with the other genes, PTB1, exhibiting the highest average Ct value of 32.08 and a wider range, was unsuitable as a reference gene. Finally, nine candidate reference genes were evaluated with the four algorithms in this study, and their mean Ct values ( Figure 3) and Ct value ranges were sufficient for experimental needs (Table S2). The stability of genes may be directly reflected in their Ct ranges. For instance, the expression of TBP1 might be more stable than others according to the relatively narrower Ct range. Importantly, these results were consistent with the outputs calculated by geNorm, NormFinder, BestKeeper, and RefFinder ( Figure 3 and Tables 1-3). Moreover, these results also reveal that none of the nine candidates could be expressed constantly in all tissues and under different conditions in G. microphylla.
As reference gene determination programs, geNorm evaluated the stability of reference genes with the average pairwise variation [16], and NormFinder exhibited the expression stability of reference genes by analyzing their intra-and intergroup variation [17]. However, the CV and SD values were the key factors in determining the stability rankings of the reference genes obtained by BestKeeper [18]. Thus, the ranking results were not identical and were reasonable. In this study, TBP1 and EIF4A1 were the most stable genes, and GAPA and GAPB were the most unstable genes under different conditions in geNorm analysis, which was relatively consistent with the results returned by NormFinder, and thus further ensured the accuracy and reliability of the analysis results. However, EF-1α, the most unstable reference gene in the PEG group from geNorm and NormFinder analyses, was ranked at a top position in the BestKeeper analysis. Similar research findings were reported in the literature on the selection and validation of reference genes in other plant species such as R. yunnanensis [20] and sorghum bicolor [46]. Fortunately, BestKeeper still showed some conformance to the top three most stable reference genes with geNorm and Normfinder (Figure 4; Tables 2 and 3). Finally, a widely used web-based tool, RefFinder, was used to obtain a consensus stability ranking of each candidate gene under the given conditions according to the geometric mean of the attributed weights of each gene [19]. For the RefFinder analysis, TBP1, EIF4A1, and CBP20 were ranked as the top three most stable reference genes, similar to geNorm, Normfinder, and BestKeeper. Synthetically, TBP1, EIF4A1, CBP20, and TUB1 were the four most frequent and stable genes ranked by the four programs, with TBP1 occurring 38 times, EIF4A1 36 times, CBP20 25 times, and TUB1 21 times. Among them, the top four most stable genes returned by BestKeeper were TBP1 nine times, EIF4A1 eight times, CBP20 six times, and TUB1 six times. Correspondingly, the top four most unstable genes in all groups were GAPA occurring 40 times, GAPB occurring 39 times, eIF2 occurring 23 times, and EF-1α occurring 22 times. The most unstable gene seemed to be determined from the top four most unstable genes in each group (Table S3).
According to the stability value, TBP1, EIF4A1, and CBP20 were considered as the most stable internal reference genes in all samples. However, the most stable reference genes were not identical in different groups; they differed even among organs in the same group. Therefore, we suggest that it is better to choose the best reference genes for a specific case rather than using a common one for normalization, although this is slightly laborious and time-consuming. Theoretically, the proteins encoded by housekeeping genes are used either to maintain the cell structure or participate in basic cellular metabolic processes, and they should be stably expressed regardless of the tissue type or the physiological state. However, some housekeeping genes have been proven to have poor expression stability in specific experimental conditions and could not act as internal reference genes for RT-qPCR analysis [7,9]. In this study, the traditional housekeeping genes TBP1, EIF4A1, and CBP20 exhibited reasonably good stability under the given conditions. To obtain a more accurate normalization result, an increasing number of studies have applied multiple reference genes rather than a single reference gene [46,47].
In response to environmental stimuli, the expression level of ACO1 (encoding ethyleneforming enzyme) and CSD2 (encoding superoxide radicals detoxifying enzyme) would change in a wide range in G. microphylla. Thus, these two stress-responding genes were chosen to validate the suitability of the reference genes selected in our study. The most stable and unstable reference genes recommended by RefFinder in the roots, stems, and leaves under cold treatment were used to normalize ACO1 and CSD2. The relative expression levels of ACO1 and CSD2 involved in stress responses had significant differences when normalized with the most stable and unstable genes (Figure 5a). In ABA treated roots, the relative expression levels of ACO1 showed a similar level (no significant difference) when normalized with the two most stable genes alone or in combination. Still, a lower expression level (significant difference) was observed when normalized with the least stable gene, GAPA (Figure 5b), suggesting that the selection and confirmation of suitable and stable reference genes were particularly critical for the proper normalization for the RT-qPCR data in G. microphylla.

Conclusions
Our study evaluated the expression stability of nine candidate reference genes selected from transcriptome data of G. microphylla under a wide range of experimental conditions, five types of abiotic stress (cold, heat, drought, salinity, and heavy-metal) and three types of hormone stimulus (MeJA, ABA, and SA). According to the comprehensive results analyzed by four widely used programs (geNorm, NormFinder, BestKeeper, and RefFinder), TBP1 and EIF4A1 were the most stable genes, and GAPA and GAPB were the most unstable genes, in all groups. Meanwhile, the most stable genes varied among different conditions and tissues (root, stem, leaf, and flower), suggesting that it is better to choose the best reference gene for a specific case rather than using a common one, and normalization with two-or multi-gene combinations is encouraged. In addition, the expression analysis of ACO1 and CSD2 emphasized the importance of selecting suitable and stable genes for the normalization of gene expression analysis using RT-qPCR. This study is the first report on the selection and validation of reference genes in G. microphylla and related species of Gleditsia. It provides an essential foundation for future research on gene expression analyses using RT-qPCR.
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/genes13071227/s1: Figure S1: Melting curves for the candidate reference genes and target genes; Figure S2: Standard curves of the candidate reference genes and target genes; Figure S3: Average expression stability M value of nine candidate reference genes calculated with GeNorm by organ; Figure S4: The pairwise variation (V) measure of the nine candidate reference genes using GeNorm by group. Table S1: Information of candidate reference genes, target genes, and primers for RT-qPCR; Table S2: Ct values for nine candidate reference genes in organs under different abiotic stresses or hormone treatments; Table S3: Expression stability ranking of the nine candidate reference genes estimated by geNorm, NormFinder, BestKeeper and RefFinder.