Screening and Validation of Internal Reference Genes for Quantitative Real-Time PCR Analysis of Leaf Color Mutants in Dendrobium officinale

Leaf color mutants (LCMs) are important resources for studying diverse metabolic processes such as chloroplast biogenesis and differentiation, pigments’ biosynthesis and accumulation, and photosynthesis. However, in Dendrobium officinale, LCMs are yet to be fully studied and exploited due to the unavailability of reliable RGs (reference genes) for qRT-PCR (quantitative real-time reverse transcription PCR) normalization. Hence, this study took advantage of previously released transcriptome data to select and evaluate the suitability of ten candidate RGs, including Actin (Actin), polyubiquitin (UBQ), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), elongation factor 1-α (EF1α), β-tubulin (β-TUB), α-tubulin (α-TUB), 60S ribosomal protein L13-1 (RPL13AD), aquaporin PIP1-2 (PIP1-2), Intima protein (ALB3) and Cyclin (CYCB1-2) for normalizing leaf color-related genes’ expression levels via qRT-PCR. Stability rankings analysis via common software Best-Keeper, GeNorm, and NormFinder disclosed that all ten genes met the requirements of RGs. Of them, EF1α exhibited the highest stability and was selected as the most reliable. The reliability and accuracy of EF1α were confirmed through qRT-PCR analysis of fifteen chlorophyll pathway-related genes. The expression patterns of these genes via EF1α normalization were consistent with the results by RNA-Seq. Our results offer key genetic resources for the functional characterization of leaf color-related genes and will pave the way for molecular dissection of leaf color mutations in D. officinale.


Introduction
Leaf color mutations are common and easily detectable morphological phenotype variations in higher plants [1,2]. They are usually expressed at the seedling stage and are generally categorized into many types, including greenish-white, albino, light green, greenish-yellow, white emerald, etiolation, striped, yellow-green, purple, dark green, and brown [3][4][5]. Besides, LCMs can also be classified into four types: total chlorophyll deficient type, total chlorophyll increased type, chlorophyll a deficient type, and chlorophyll b deficient type [6]. Genes involved in leaf color often play critical roles in the biosynthesis of photosynthetic pigments and chloroplast development [1,6]. Therefore, LCMs represent essential resources for insight into plant physiology and metabolism as color mutations mostly affect plants' photosynthetic efficiency, which results in poor growth and development performances and considerable economic losses [1,2]. Consequently, LCMs have become ideal materials for investigating pigments' metabolism, chloroplast formation and differentiation, photosynthesis, biotic and abiotic stress responses, etc. [1,2]. In rice, studies on LCMs have significantly contributed to the crop improvement [1]. For instance, LCMs and senescence-associated mutants have been identified, and the functional analysis of fore, it is essential to identify stable and reliable RG(s) for specific traits and experimental conditions in each plant species [18,31,32].
In this study, based on the previously released RNA-seq data [14], we selected ten candidate RGs and investigated their expression stability using common statistical algorithms (BestKeeper, GeNorm, and NormFinder). As a result, we identified and validated the most suitable RG for qRT-PCR normalization of leaf color-related genes in D. officinale. Our findings provide optimal RG for analyzing the expression of leaf color genes in D. officinale. Moreover, they will promote genomics studies on LCMs in D. officinale and contribute to understanding color mutation in plants.

Plant Materials and Sampling
The wild-type (green leaf, CK) and a leaf color-mutant (yellow) D. officinale tissue culture materials were used in this study (Figure 1). When the tissue cultures were at the seedlings stage, leaves were sampled from ten plants of each genotype, mixed, and immediately frozen in liquid nitrogen. Subsequently, the collected samples were stored at −80 • C in a refrigerator until further experimentations. Three biological and experimental repeats were achieved for each analysis.
Genes 2023, 14, x FOR PEER REVIEW 3 of 13 stable expressions at different developmental stages or under diverse experimental conditions in all plant organs [18,27,28]. But, other studies in many species have revealed that RGs do not always exhibit stable expression levels, specifically under changing environments [28][29][30]. Therefore, it is essential to identify stable and reliable RG(s) for specific traits and experimental conditions in each plant species [18,31,32].
In this study, based on the previously released RNA-seq data [14], we selected ten candidate RGs and investigated their expression stability using common statistical algorithms (BestKeeper, GeNorm, and NormFinder). As a result, we identified and validated the most suitable RG for qRT-PCR normalization of leaf color-related genes in D. officinale. Our findings provide optimal RG for analyzing the expression of leaf color genes in D. officinale. Moreover, they will promote genomics studies on LCMs in D. officinale and contribute to understanding color mutation in plants.

Plant Materials and Sampling
The wild-type (green leaf, CK) and a leaf color-mutant (yellow) D. officinale tissue culture materials were used in this study ( Figure 1). When the tissue cultures were at the seedlings stage, leaves were sampled from ten plants of each genotype, mixed, and immediately frozen in liquid nitrogen. Subsequently, the collected samples were stored at −80 °C in a refrigerator until further experimentations. Three biological and experimental repeats were achieved for each analysis.

Total RNA Extraction and cDNA Synthesis
Total RNA from all samples was extracted using the RNAprep pure polysaccharide polyphenol plant total RNA extraction kit according to the manufacturer's instructions. The extracted RNA concentration and quality (OD 260 /OD 280 value) were assessed with Nanodrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). Finally, the integrity of total RNA was confirmed through electrophoresis (1% agarose gel).
The GeneAmp ® PCR System 9700 (Applied Biosystems, Foster City, CA, USA) was used for cDNA construction. The reverse transcription was achieved using Hiscript II QRT super mix according to the manufacturer's instructions for users. The reaction conditions were 15 min at 42 • C and 5 s at 85 • C. Each synthesized cDNA was ten-time diluted (with nuclease-free water) and stored in the refrigerator at −20 • C.

qRT-PCR Analysis and Validation of the Most Reliable RG
The qRT-PCR was conducted on Lightcycle ® 480 type II fluorescent quantitative PCR (Roche, Indianapolis, USA) using the quantifast ® SYBR ® green PCR kit and following instructions by the manufacturer. The PCR reaction was a mixture of cDNA 1 µL; Upstream and downstream primers (10 µMol·L −1 ) 0.2 µL for each; 2× QuantiFast ® SYBR ® Green PCR Master Mix 5 µL; and Nuclease-free H 2 O 3.6 µL. The reaction procedure was as follows: pre-denaturation at 95 • C for 5 min; denaturation at 95 • C for 10 s (40 cycles); and annealing at 60 • C for 30 s. Three technical and biological repetitions were set for each gene. At the end of the cycle, the melting curve was used to detect the product specificity, and each gene's Ct (cycle threshold) values were computed automatically and saved [24,33].
The reliability of the most reliable RG was confirmed through qRT-PCR analysis of the relative expression of the fifteen selected chlorophyll pathway-related genes, followed by expression level calculation via the 2 −∆∆Ct method [34]. The internal control was the most stable RG.

Stability and Statistical Analyses
The Ct values were considered as relative quantities to perform gene expression stability analysis [23,33]. Three commonly used software, including NormFinder [16,18], BestKeeper [34], and geNorm [31]. The software were run following instructions in their respective manuals. Finally, we applied the GM (geometric mean) method to fit the results from the three software and generate a comprehensive stability ranking for the candidate RGs [35]. Data Processing System, GraphPad Prism v9.0.0121 (GraphPad 159 Software Inc., La Jolla, CA, USA), and Microsoft Excel 2016 were used for data analysis and graphing [35,36]. Statistical differences were obtained via an independent t-test at p < 0.05. Finally, we used TBtools software to construct the heatmap of genes [37].

Primers Specificity and Expression Profiles Analyses of Candidate RGs
To uncover a suitable RG for normalizing the relative expression levels of leaf colorrelated genes in D. officinale via qRT-PCR, we selected ten candidate RGs (Table S1) from transcriptome data and subjected them to qRT-PCR analysis. The melting curves showed that the primers of each potential RG were highly specific ( Figure 3). All ten genes exhibited a single peak with no primer dimer, and all amplicons showed good repeatability, confirming the accuracy, reliability, and specificity of primers.
In order to assess the expression levels of each potential RGs, we automatically computed their respective Ct values. The results showed that the Ct values of the ten genes ranged from 18.858 (GAPDH) to 26.915 (RPL13AD) on average ( Figure 4) Figure 4). The summary of the variation in expression levels (Ct values) of all candidate RGs is presented in Table S4. The coefficient of variation (CV) of Actin, UBQ10, GAPDH, EF1α, β-TUB, α-TUB, RPL13AD, PIP1-2, ALB3, CYCB1-2 was 2.33%, 2.13%, 2.18%, 1.57%, 1.54%, 1.45%, 1.30%, 2.82%, 3.11%, 1.14%, respectively (Table S4). Genes 2023, 14, x FOR PEER REVIEW 6 of 13 In order to assess the expression levels of each potential RGs, we automatically computed their respective Ct values. The results showed that the Ct values of the ten genes ranged from 18.858 (GAPDH) to 26

Stability Analysis of Candidate Internal RGs
To determine the stability rankings of the ten potential RGs, three independent software (BestKeeper, GeNorm, and NormFinde) were used, and the results were integrated via the GM method to obtain the comprehensive stability ranking. Firstly, we used geNorm to compute the expression stability measured (M) of each potential RG based on average pairwise expression ratios. A gene can be selected as an internal RG if its M value is below 1.5, and the lower the M value, the more the gene is stable [24]. The analysis by geNorm showed that the M values of all ten genes were >0.9, indicating they both met the basic requirements for an internal RG ( Figure 5A). EF1α was the most stable gene, with an M value of 0.395, followed by CYCB1-2 (M value of 0.439) ( Figure 5A). PIP1-2 and ALB3 recorded the highest stability M values of 0.753 and 0.861, respectively. Usually, a single RG is insufficient to achieve high accuracy. Accordingly, sometimes using two or more RGs for precise and reliable normalization is recommended. The optimal number of RGs is determined by pairwise variation value (Vn/n + 1, V value). If Vn/n + 1 ≥ 0.15, the adequate number of RGs equals n + 1. In contrast, if Vn/n + 1 < 0.15, the number is n [23,31]. The V value of all pairwise comparisons was less than 0.15 ( Figure 5D), indicating two RGs might be required. Thus, the M values indicate EF-1α + CYCB1-2 is the best combination for accurate normalization of leaf color-related genes' expression levels in D. officinale.
NormFinder serves to compute the SV (stability value) of RGs and unveil the optimal number of RGs for precise normalization through analysis of intra-and inter-group variations [24]. The most stable gene is one that exhibits a lower expression level than the average SVs. The stability ranking by NormFinder and GeNorm analyses were somewhat similar, with a little variation in the classification of the top four genes ( Figure 5B). Interestingly, EF1α (SV = 0.064) occupied the first rank, confirming it was the most stable RG. β-TUB (0.115), RPL13AD (0.117), and CYCB1-2 (0.143) ranked second, third, and fourth in terms of stability by NormFinder, respectively ( Figure 5B). PIP1-2 (SV = 0.481) and ALB3 (SV = 0.57) occupy the last ranking positions as per the GeNorm results.

Stability Analysis of Candidate Internal RGs
To determine the stability rankings of the ten potential RGs, three independent software (BestKeeper, GeNorm, and NormFinde) were used, and the results were integrated via the GM method to obtain the comprehensive stability ranking. Firstly, we used geNorm to compute the expression stability measured (M) of each potential RG based on average pairwise expression ratios. A gene can be selected as an internal RG if its M value is below 1.5, and the lower the M value, the more the gene is stable [24]. The analysis by geNorm showed that the M values of all ten genes were >0.9, indicating they both met the basic requirements for an internal RG ( Figure 5A). EF1α was the most stable gene, with an M value of 0.395, followed by CYCB1-2 (M value of 0.439) ( Figure 5A). PIP1-2 and ALB3 recorded the highest stability M values of 0.753 and 0.861, respectively. Usually, a single RG is insufficient to achieve high accuracy. Accordingly, sometimes using two or more RGs for precise and reliable normalization is recommended. The optimal number of RGs is determined by pairwise variation value (Vn/n + 1, V value). If Vn/n + 1 ≥ 0.15, the adequate number of RGs equals n + 1. In contrast, if Vn/n + 1 < 0.15, the number is n [23,31]. The V value of all pairwise comparisons was less than 0.15 ( Figure 5D), indicating two RGs might be required. Thus, the M values indicate EF-1α + CYCB1-2 is the best combination for accurate normalization of leaf color-related genes' expression levels in D. officinale.
NormFinder serves to compute the SV (stability value) of RGs and unveil the optimal number of RGs for precise normalization through analysis of intra-and inter-group variations [24]. The most stable gene is one that exhibits a lower expression level than the average SVs. The stability ranking by NormFinder and GeNorm analyses were somewhat similar, with a little variation in the classification of the top four genes ( Figure 5B). Interestingly, EF1α (SV = 0.064) occupied the first rank, confirming it was the most stable RG. β-TUB (0.115), RPL13AD (0.117), and CYCB1-2 (0.143) ranked second, third, and fourth in terms of stability by NormFinder, respectively ( Figure 5B). PIP1-2 (SV = 0.481) and ALB3 (SV = 0.57) occupy the last ranking positions as per the GeNorm results. BestKeeper helps to evaluate the stability index of RGs by calculating mainly the SD (standard deviation) [38]. A gene is stable if its SD is lower than 1.0 [24,38]. The results by BestKeeper indicated that the SD of all ten genes was lower than 0.7, supporting that they are suitable for use as internal RGs ( Figure 5C). CYCB1-2 (SV = 0.2002) and EF1α (SV = 0.2386) ranked first and second in terms of stability, respectively ( Figure 5C). Both three software revealed that PIP1-2 and ALB3 were the least stable genes.
We applied the GM method to integrate the results from the three software and determine the comprehensive stability ranking of the ten RGs. As presented in Table 2, GM values confirmed that EF1α was the most suitable and reliable RG, followed by CYCB1-2, RPL13AD, β-TUB, UBQ10, α-TUB, GAPDH, Actin, PIP1-2, and ALB3.  BestKeeper helps to evaluate the stability index of RGs by calculating mainly the SD (standard deviation) [38]. A gene is stable if its SD is lower than 1.0 [24,38]. The results by BestKeeper indicated that the SD of all ten genes was lower than 0.7, supporting that they are suitable for use as internal RGs ( Figure 5C). CYCB1-2 (SV = 0.2002) and EF1α (SV = 0.2386) ranked first and second in terms of stability, respectively ( Figure 5C). Both three software revealed that PIP1-2 and ALB3 were the least stable genes.
We applied the GM method to integrate the results from the three software and determine the comprehensive stability ranking of the ten RGs. As presented in Table 2, GM values confirmed that EF1α was the most suitable and reliable RG, followed by CYCB1-2, RPL13AD, β-TUB, UBQ10, α-TUB, GAPDH, Actin, PIP1-2, and ALB3.

Validation of EF1α as the Most Reliable and Stable Internal RG
To confirm the reliability and stability of EF1α for accurate normalization of relative expression levels of leaf color-related genes, we investigated the expression of fifteen selected genes from the chlorophyll biosynthesis pathway ( Figure 2B, Table S2) via qRT-PCR. As shown in Figure 6, the expression patterns of the targeted genes via qRT-PCR and RNA-Seq were consistent with identical regulation patterns. Dca001861, Dca002994, Dca004000, Dca010603, Dca013248, Dca020330, Dca020854, and Dca024845 were significantly up-regulated as per the results of transcriptomics analysis ( Figure 6). Among them, Dca024845 and Dca020330 were the most significantly up-regulated, suggesting they might play crucial roles in yellow leaf development in D. officinale.

Validation of EF1α as the Most Reliable and Stable Internal RG
To confirm the reliability and stability of EF1α for accurate normalization of relative expression levels of leaf color-related genes, we investigated the expression of fifteen selected genes from the chlorophyll biosynthesis pathway ( Figure 2B, Table S2) via qRT-PCR. As shown in Figure 6, the expression patterns of the targeted genes via qRT-PCR and RNA-Seq were consistent with identical regulation patterns. Dca001861, Dca002994, Dca004000, Dca010603, Dca013248, Dca020330, Dca020854, and Dca024845 were significantly up-regulated as per the results of transcriptomics analysis ( Figure 6). Among them, Dca024845 and Dca020330 were the most significantly up-regulated, suggesting they might play crucial roles in yellow leaf development in D. officinale. Figure 6. The relative expression levels of fifteen chlorophyll pathway-related genes via normalization with the most suitable RG (reference gene) EF1α. *, **, ***, **** indicate statistical differences at p < 0.05, 0.01, 0.001, and 0.0001, respectively.

Discussion
Leaf color is a critical agronomic trait, and its variation significantly affects global plant metabolism. Principally, leaf color mutations induce less photosynthetic efficiency, causing poor growth, reduced yield, and important economic losses [1,2]. Accordingly, leaf color mutants have become materials of considerable research interest. They are ideal for studying variations in plant metabolism and physiology, especially the molecular mechanisms governing chloroplast biogenesis and differentiation, pigment synthesis and accumulation, photosynthesis, diverse stress response, etc. [1,2]. In D. officinale, previous studies have shown that leaf color mutation is associated with a reduction in plant height, stem diameter, and the number of chloroplasts; low chlorophyll content; and altered chloroplast structure [14,15]. In addition, they have provided insight into the underlying molecular mechanisms through comparative transcriptome analysis and unveiled DEGs [14]. However, the identified potential candidate genes have not been validated, and functional characterizations are yet to be conducted due to the unavailability of a reliable RG for qRT-PCR analysis. Thus, the present study took advantage of the RNA-seq data to comprehensively select candidate RGs, examine their stability, and validate the most stable to promote genomics studies on leaf color mutants in D. officinale. This approach has been

Discussion
Leaf color is a critical agronomic trait, and its variation significantly affects global plant metabolism. Principally, leaf color mutations induce less photosynthetic efficiency, causing poor growth, reduced yield, and important economic losses [1,2]. Accordingly, leaf color mutants have become materials of considerable research interest. They are ideal for studying variations in plant metabolism and physiology, especially the molecular mechanisms governing chloroplast biogenesis and differentiation, pigment synthesis and accumulation, photosynthesis, diverse stress response, etc. [1,2]. In D. officinale, previous studies have shown that leaf color mutation is associated with a reduction in plant height, stem diameter, and the number of chloroplasts; low chlorophyll content; and altered chloroplast structure [14,15]. In addition, they have provided insight into the underlying molecular mechanisms through comparative transcriptome analysis and unveiled DEGs [14]. However, the identified potential candidate genes have not been validated, and functional characterizations are yet to be conducted due to the unavailability of a reliable RG for qRT-PCR analysis. Thus, the present study took advantage of the RNA-seq data to comprehensively select candidate RGs, examine their stability, and validate the most stable to promote genomics studies on leaf color mutants in D. officinale. This approach has been widely used in several plants to identify suitable RGs for specific traits and experimental conditions [20,24,[39][40][41].
The efficiency of qRT PCR results depends on the reliability of internal RGs. Only a stable expression of internal RG can guarantee the accuracy of the results. Previous studies have demonstrated that internal RGs are species-specific, and their stability varies according to the traits and experimental conditions [20,24,[39][40][41]. For instance, different RGs have been identified for studying traits related to pigment synthesis in different plant species. In Catalpa fargesii, CfMADH and CfEF-1 have been identified as the most reliable RGs to be used individually or in combination for normalizing the expression levels of leaf color-related genes [40]. In Lagerstroemia indica, LiEF1α-2 and LiEF1α-1 were identified as the most suitable for analyzing the expression levels of leaf color-related genes [42]. In wheat, 18S rRNA is used as the internal RG to calibrate the expression of leaf color-related candidate genes [26]. In Lilium regale, LrActin2 was identified as the best RG for normalizing photosynthesis-related candidate genes' expression levels [25]. Herein, we investigated ten genes, including EF1α, CYCB1-2, RPL13AD, β-TUB, UBQ10, α-TUB, GAPDH, Actin, PIP1-2, and ALB3. The melting curves showed that all ten genes exhibited a single peak with no primer dimer, and all amplicons showed good repeatability. These results indicate that the designed primers were highly specific, accurate, and reliable. These primers could be directly employed in future studies. The coefficient of variation (CV) of the Ct values of all ten genes was lower than 3%, confirming the repeatability of the experiment. Expression stability analysis via GeNorm, NormFinder, and BestKeeper revealed that these ten genes could be used as internal RGs to calibrate leaf color-related genes' expression levels in D. officinale. The comprehensive ranking analysis showed that EF1α was the most stable.These results suggest that EF1α is the most reliable RG for calibrating the expression levels of leaf color-related genes in D. officinale. Using this RG may promote the molecular dissection of the regulatory network of leaf color mutations in D. officinale. We did not investigate the expression stability of these RGs in other plant organs since this study aimed to identify a reliable RG for leaf color specifically.
To verify the suitability and reliability of EF1α, we analyzed fifteen chlorophyll pathway-related DEGs' expression levels via qRT qPCR. The expression patterns of the targeted genes via EF1α normalization were consistent with the transcriptome sequencing results. This result further confirms the accuracy of using EF1α as the primary RG for qRT-PCR calibration of leaf color-related genes' expression levels in D. officinale. Furthermore, of the fifteen analyzed genes, Dca024845 and Dca020330 were the most significantly upregulated in the yellow leaf mutant genotype. Accordingly, these genes may represent key candidate genes for deciphering chlorophyll synthesis-related molecular mechanisms in D. offinale. In rice, most leaf color mutants are directly or indirectly involved in chlorophyll metabolism [1]. Any defects in the development of chloroplasts also impair the stability of chlorophyll and other photosynthetic pigments, eventually changing the green color of the leaf [1]. In general, leaf color mutation is caused by complex mechanisms, including abnormal chlorophyll metabolism pathways (mutations of genes related to chlorophyll biosynthesis or degradation pathways and mutations of genes related to heme metabolism pathway), abnormal chloroplast development and differentiation, defects in carotenoid metabolism, abnormal anthocyanin biosynthesis and degradation (mutations of structural and regulation genes) [2]. Therefore, functional characterization of these genes is required to reveal their roles and gain insight into leaf coloring molecular mechanisms in D. officinale. With the rapid development of genomics tools, the discovery of this RG may accelerate the molecular dissection of leaf coloring mechanisms in D. officinale.

Conclusions
In summary, this study investigates ten RGs for qRT-PCR normalization of expression levels of leaf color-related candidate genes in D. officinale. Stability rankings analysis showed that all ten genes, including Actin (Actin), polyubiquitin (UBQ), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), elongation factor 1-α (EF1α), β-tubulin (β-TUB), α-tubulin (α-TUB), 60S ribosomal protein L13-1 (RPL13AD), aquaporin PIP1-2 (PIP1-2), Intima protein (ALB3) and Cyclin (CYCB1-2) met the requirement of an RG. Of them, EF1α was the most stable, and its accurability was confirmed through the qRT-PCR calibration of the expression levels of fifteen chlorophyll pathway-related DEGs. The regulation patterns of the fifteen chlorophyll pathway-related genes were identical to the RNA-seq results. Dca024845 and Dca020330 were the most significantly up-regulated genes in the yellow leaf mutant genotype. Therefore, they were proposed as candidate genes for dissecting leaf coloring molecular mechanisms in D. officinale. Our results provide key genetic resources for validating leaf color mutants-related candidate genes and perform functional genomics studies. Furthermore, they may provide comprehensive guidelines for uncovering suitable RGs in other plant species.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/genes14051112/s1, Table S1: List of the ten candidate reference genes with their FPKM expression values; Table S2: List of the fifteen chlorophyll pathwayrelated genes with their FPKM expression values; Table S3: Primers of the fifteen chlorophyll pathway-related genes for qRT-PCR analysis; Table S4: Summary of the qRT-PCR Ct value of ten candidate genes.