Selection and Validation of Optimal RT-qPCR Reference Genes for the Normalization of Gene Expression under Different Experimental Conditions in Lindera megaphylla

Lindera megaphylla, a broad-leaved evergreen that is used as a landscape ornamental plant and medicinal plant, is an ecologically important and dominant tree species. However, little is known about the molecular mechanisms of its growth, development, and metabolism. The selection of suitable reference genes is critical for molecular biological analyses. To date, no research on reference genes as a foundation for gene expression analysis has been undertaken in L. megaphylla. In this study, 14 candidate genes were selected from the transcriptome database of L. megaphylla for RT-qPCR assay under different conditions. Results showed that helicase-15 and UBC28 were most stable in different tissues of seedlings and adult trees. For different leaf developmental stages, the best combination of reference genes was ACT7 and UBC36. UBC36 and TCTP were the best under cold treatment, while PAB2 and CYP20-2 were the best under heat treatment. Finally, a RT-qPCR assay of LmNAC83 and LmERF60 genes were used to further verify the reliability of selected reference genes above. This work is the first to select and evaluate the stability of reference genes for the normalization of gene expression analysis in L. megaphylla and will provide an important foundation for future genetic studies of this species.


Introduction
Lindera megaphylla is a predominant, broad-leaved, and aromatic evergreen tree species belonging to the Lauraceae family and is widely distributed in the subtropical and warmtemperate zones of China. L. megaphylla has not only ecological and ornamental value, but also medicinal and therapeutic value as a source of an essential oils, spices, and drugs [1][2][3]. For example, d-dicentrine, an aporphine alkaloid, is isolated from the root of L. megaphylla and has potential antitumor activity [4,5]. These trees are rich in terpenoids, alkaloids, and flavonoids, many of which could be used to make pesticides or industrial feedstocks, while its volatile compounds, mainly terpenoids, have strong bactericidal ability. These trees can also help to improve air quality [6], a property that could be improved through molecular breeding. To synthesize antibacterial compounds, it is necessary to first explore the related regulatory genes and to analyze their functions. To date, studies of L. megaphylla have largely focused on the cultivation of seedlings, various kinds of biotic and abiotic stress responses and the analysis of volatile substances and their potential applications [7][8][9][10][11][12], while few studies have focused on the molecular biology of L. megaphylla due to a lack of genomic information. Gene expression analysis circumvents that lack of a sequenced genome to explore the molecular mechanisms underlying transcriptional regulation of phenotype. Transcriptome datasets derived from different tissues and differently aged binding protein 2-like (PAB2). These 14 candidate genes were then assessed for stability of expression under specific conditions, including different tissues of one-year-old seedlings (roots, stems, and leaves), tissues of 10-year-old trees (leaf buds, young stems, young seeds, young leaves, and mature leaves), 16 different leaf developmental stages and under different temperature stresses (cold and heat). Finally, the NAC and ERF genes, which are from an important family of transcription factors in plants, were used to verify the reliability of the selected reference genes in different samples. Our research identified the best reference genes for RT-qPCR analysis of L. megaphylla tissues under different conditions, laying a basis for further studies of the molecular mechanisms regulating gene expression in this important tree species.

Plant Materials and Treatments
L. megaphylla trees were grown in a field at the Zhengzhou Botanical Garden located in Zhengzhou city, Henan Province, China. Young seedlings were cultivated in plant growth chambers with LED lighting for temperature treatments. In total, 6 experimental sets were cultivated for RNA extraction. The first experimental set consisted of roots, stems, and leaves from 1-year-old seedlings that displayed robust and consistent growth. The second experimental set consisted of leaf buds, young stems, young seeds, young leaves, and mature leaves from adult trees of L. megaphylla that had been growing in the natural environment for approximately 10 years. The third experimental set included leaves at different developmental stages. Growing leaves were collected roughly every 3-7 or 3-15 days from the beginning of leaf bud growth in late March until the leaves were fully mature by the end of July. The specific sampling dates were 24 March, 29 March, 1 April, 4 April, 7 April, 15 April, 22 April, 29 April, 14 May, 23 May, 4 June, 16 June, 1 July, 15 July, and 31 July. The fourth set was exposed to cold stress [41]. Furthermore, 1-year-old seedlings were grown at 65% relative humidity and under 16 h/8 h light/dark conditions in an LED plant growth incubator (Shengyuan Instrument Co., Ltd., Zhengzhou, China). The seedlings were first treated at 25 • C for 7 days as control, and leaves were collected. Then, the seedlings were transferred to 4 • C for 7 days for long-term chilling acclimation (CA). Next, the seedlings were shifted to 0 • C for an additional 7 days for longterm freezing acclimation (FA). Then, the seedlings were again moved to control conditions (25 • C) for 7 days for long-term de-acclimation (DA). The fifth and sixth experimental groups consisted of leaves collected from 1-year-old seedlings that were treated with cold and heat, respectively. For cold treatments, the seedlings were cultivated in an LED plant growth incubator at 25 • C, 4 • C, 0 • C, −4 • C, or −6 • C for 24 h. For heat treatments, the seedlings were cultivated in an incubator at 25 • C, 30 • C, 35 • C, 40 • C, or 45 • C for 24 h. Data regarding all six sample sets described above are summarized in Table 1. All collected samples were immediately frozen in liquid nitrogen and stored at −80 • C. Three independent biological replicates were collected for each sample.

Total RNA Extraction and cDNA Synthesis
Total RNA was extracted from 0.05 g samples using a Quick RNA isolation kit (HUAYUEYANG Biotechnology, Beijing, China) according to the manufacturer's protocol [42]. The RNA integrity, purity and concentration were assessed using 2% (w/v) agarose gel electrophoresis and a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). Total RNA (1 µg) with A 260 /A 280 and A 260 /A 230 ratios greater than 1.8 was used for first-strand cDNA synthesis using the Evo M-MLV RT Kit with gDNA Clean for qPCR according to the manufacturer's instructions (Accurate Biotechnology, Changsha, China). Specifically, a 10 µL reaction system (1 µg total RNA, 2 µL 5 × gDNA Clean Reaction Mix, to a total volume of 10 µL with RNase-free water) was subjected to 42 • C for 2 min. Then, a 20 µL reaction system (10 µL of first reaction solution, 4 µL 5 × Evo M-MLV RT Reaction Mix and 6 µL RNase-free water) was subjected to 37 • C for 15 min and 85 • C for 5 s. Five-fold diluted cDNA was used for subsequent RT-qPCR experiments. All cDNA samples were stored at −20 • C until use.

RT-PCR and RT-qPCR Data Analysis
To verify the accuracy of the designed primers, each pair was used for RT-PCR amplification. Each 20 µL reaction system contained 2 µL of 10 µM forward and reverse primers, 1 µL cDNA, 10 µL 2 × Rapid Taq Master Mix (Vazyme, Nanjing, China) and 7 µL ddH 2 O. The reaction was carried out at 95 • C for 3 min, followed by 35 cycles of 95 • C for 30 s, 60 • C for 30 s and 72 • C for 15 s, with a final extension at 72 • C for 2 min. The PCR products were visualized by 1% (w/v) agarose gel electrophoresis.
To monitor the E-value, the cDNA templates from all samples were serially diluted five-fold (cDNA:water, v:v). RT-qPCR was performed for each pair of primers to obtain Ct values and to establish a standard curve; the R 2 , slope and E-values were calculated with Microsoft Office Excel 2019 using the following formula: E = (5 −1/slope − 1) × 100% [43].
All RT-qPCRs were performed with an Applied Biosystems TM (ABI) QuantStudio TM 5 real-time PCR system (ABI, Los Angeles, CA, USA) using the following amplification procedure: 95 • C for 30 s, followed by 40 cycles of 95 • C for 5 s and 60 • C for 30 s. A melting curve was generated at 60-95 • C. Each 20 µL RT-qPCR reaction (10 µL 2 × SYBR ® Green Pro Taq HS Premix, 2 µL 5-fold diluted cDNA, 0.4 µL 10 µM forward primer, 0.4 µL 10 µM reverse primer, 0.4 µL 4 µM ROX Reference Dye, and 6.8 µL RNase-free water) was prepared according to the instructions for SYBR ® Green Premix Pro Taq HS qPCR Kit (Rox Plus) (Accurate Biology, Shanghai, China). A negative control without the addition of cDNA was used to test for background amplification. Three technical replicates were performed for each sample, and the mean was used for RT-qPCR analysis.

Candidate Reference Gene Expression Stability Analysis
CT values were used to assess the expression levels of candidate reference genes in all samples by RT-qPCR. Four common algorithms, namely, delta Ct (∆Ct), geNorm (version 3.5), NormFinder (version 0.953) and BestKeeper (version 1.0), were used to evaluate the stability of the expression of the candidate reference genes in the different experimental groups. In geNorm and NormFinder, the M value reflects the stability of each candidate reference gene [28,29], with a smaller M value indicating higher stability. The geNorm package also determines the number of optimal reference genes based on the ratio V n/n+1 by calculating pairwise variations in the normalized factor after introducing a new internal reference gene [28]. A cut-off value of 0.15 was used for pairwise variation. If the value of V n/n+1 was less than 0.15, n was selected for the number of optimal internal reference genes; if the value of V n/n+1 was greater than 0.15, n + 1 was selected. For BestKeeper, the values of CV and SD were used to evaluate the relative expression stability of each candidate gene [29]. The smaller the CV and SD values are, the more stable the gene. Finally, the RefFinder program (http://blooge.cn/RefFinder/, accessed on 25 September, 2022) was used to comprehensively rank the candidate reference genes by ∆Ct, geNorm, NormFinder, and BestKeeper as previously described [31].

Validation of Candidate Reference Genes by RT-qPCR
The NAC and ERF transcription factors function as central switches of growth, development and various abiotic/biotic stress responses in plants [44][45][46][47][48]. Hence, LmNAC83 and LmERF60 genes were selected as targets to determine the reliability of the most stable and unstable reference genes. The primers for these two genes are shown in Supplementary  Table S2. The expression levels of these genes were calculated using the 2 −∆∆Ct method [35]. The RT-qPCR analysis was carried out using RNA from three biological replicates.

Verification of Amplicon Size, Primers Specificity and PCR Amplification Efficiency
A total of 20 candidate genes from L. megaphylla were selected as potential reference genes for the normalization of target gene transcript levels using RT-qPCR. Standard PCR amplification with primers ( Table 2, Supplementary Table S1) targeting the candidate genes was performed with reverse-transcribed cDNA from each sample as templates. Agarose gel electrophoresis indicated that all PCR products were single bands of the expected sizes, indicating that the primers were specific (Supplementary Figure S1). The melting curves, obtained after 40 cycles of amplification by RT-qPCR, showed single peaks, which also verified that the primers for the 20 candidate reference genes had strong specificity (Supplementary Figure S2). Standard curves, generated using a five-fold serial dilution for each candidate gene, had linear correlation coefficients (R 2 ) greater than 0.99 for each specific primer pair, and the amplification efficiencies of the RT-qPCR reactions were 64.42-107.17% (Table 2). Because the amplification efficiencies of ICIn (64.42%), UBQ (85.95%), SDE2 (86.75%), and EF1α (83.49%) were less than 90%, and the R 2 values of CYP95 and RHA2A were <0.99, these six candidate genes, based on the designed primers, were not appropriate for validating expression. Therefore, the remaining 14 candidate genes were used for subsequent experiments.

Transcript Abundance of Candidate Reference Genes
The transcript abundance of 14 candidate reference genes was estimated using the average cycle threshold (Ct) values for RNA extracted from different tissues or plants grown under different experimental conditions. All candidate reference genes were expressed at a wide range of transcript levels under the different experimental conditions. The gene expression level is negatively correlated with the Ct value, which means that a gene with a higher transcript level has a smaller Ct value. As shown in Figure 1

Estimation of the Stability of the Reference Genes under Different Experimental Conditions
To identify the optimal reference genes for the normalization of gene expression analysis in L. megaphylla, the stability of the 14 candidate genes was assessed by four different algorithms. The RefFinder software was used for overall ranking.

Delta Ct Method Analysis
The delta Ct (ΔCt) method ranks the stability of candidate reference genes based on the relative expression levels of "gene pairs" in each group of sample comparisons, while the mean standard deviation of gene expression differences (STDEV) is inversely proportional to its stability using the raw Ct value [34]. The stability of the transcript levels of each candidate reference gene was evaluated based on the STDEV value. The gene with the minimum STDEV value was regarded as the most stably expressed gene. The results demonstrated that the optimal reference genes were different in the different experimental sets ( Figure 2). In seedling samples, helicase-15 and PAB2 showed the lowest ΔCt values (0.43), indicating the most stability (Figure 2A), while UBC28 (0.31) was the most stably expressed gene in adult trees ( Figure 2B). In other tissues, helicase-15 (0.51) and UBC28 (0.57) were the most stable reference genes ( Figure 2C), which is consistent with the results

Estimation of the Stability of the Reference Genes under Different Experimental Conditions
To identify the optimal reference genes for the normalization of gene expression analysis in L. megaphylla, the stability of the 14 candidate genes was assessed by four different algorithms. The RefFinder software was used for overall ranking.

Delta Ct Method Analysis
The delta Ct (∆Ct) method ranks the stability of candidate reference genes based on the relative expression levels of "gene pairs" in each group of sample comparisons, while the mean standard deviation of gene expression differences (STDEV) is inversely proportional to its stability using the raw Ct value [34]. The stability of the transcript levels of each candidate reference gene was evaluated based on the STDEV value. The gene with the minimum STDEV value was regarded as the most stably expressed gene. The results demonstrated that the optimal reference genes were different in the different experimental sets (Figure 2). In seedling samples, helicase-15 and PAB2 showed the lowest ∆Ct values (0.43), indicating the most stability (Figure 2A), while UBC28 (0.31) was the most stably expressed gene in adult trees ( Figure 2B). In other tissues, helicase-15 (0.51) and UBC28 (0.57) were the most stable reference genes ( Figure 2C), which is consistent with the results for seedlings and adult trees. ACT7 (0.58) and UBC36 (0.59) were more stable across different leaf developmental stages ( Figure 2D). Data analyses from the entire growth cycle indicated that ubiquinone (0.67) and UBC36 (0.68) were the most stable ( Figure 2E). UBC36 showed good stability in all three cold stress sets ( Figure 2F-H). PAB2 (0.69) had the highest stability under heat treatment ( Figure 2I). Across all temperature stresses, PAB2 (0.6) was the most stably expressed gene ( Figure 2J). For all samples, ubiquinone (0.74), EF2 (0.75) and PAB2 (0.76) showed the most stability ( Figure 2K). TUA had relatively higher Ct values, indicating that it was the least stable reference gene in most of the experimental sets (Figure 2A-K). for seedlings and adult trees. ACT7 (0.58) and UBC36 (0.59) were more stable across different leaf developmental stages ( Figure 2D). Data analyses from the entire growth cycle indicated that ubiquinone (0.67) and UBC36 (0.68) were the most stable ( Figure 2E). UBC36 showed good stability in all three cold stress sets ( Figure 2F-H). PAB2 (0.69) had the highest stability under heat treatment ( Figure 2I). Across all temperature stresses, PAB2 (0.6) was the most stably expressed gene ( Figure 2J). For all samples, ubiquinone (0.74), EF2 (0.75) and PAB2 (0.76) showed the most stability ( Figure 2K). TUA had relatively higher Ct values, indicating that it was the least stable reference gene in most of the experimental sets (Figure 2A-K).

geNorm Analysis
The stability of the reference genes was ranked by calculating the average expression stability values (M value) using the geNorm program, taking into account only similar intergroup variation [35].
Genes with an M value less than 1.5 were considered stably expressed, with smaller M values indicating a more stable gene [28]. The geNorm analysis results for the 14 candidate genes in the different experimental sets are shown in Figure 3A-K. The M values of the 14 candidate reference genes were less than 1.5 under the different experimental conditions (Figure 3). EF2 and PAB2 had the highest stability in seedlings, with M values of 0.012 ( Figure 3A), while CYP26-2 (0.076) and helicase-15 (0.076) were most stably expressed in the adult tree ( Figure 3B). Between the different tissue sets, PAB2 and helicase15 were the optimal candidate genes, which is similar to the results in seedlings and adult trees ( Figure 3C). The genes PAB2 and helicase-15 had the lowest stability values (0.285), which is consistent with the results of the ∆Ct analysis. The two most stably expressed genes among the different leaf developmental stages were similar to those for the entire growth cycle, namely, UBC36 and UBC7 ( Figure 3D,E). These results are also similar to those of the ∆Ct analysis. UBC36 was more stable than the other candidate reference genes under the cold treatments ( Figure 3F-H) and was the same regardless of whether the cold treatment lasted for 7 days or 24 h. GAPDH (0.269), CYP20-2 (0.269), and PPR (0.306) were more stable than the other candidate genes under heat treatment ( Figure 3I). However, PPR and PAB2 (0.343) exhibited the strongest stability under temperature stress ( Figure 3J). These results are consistent with the results of the ∆Ct analysis. UBC36 and UBC7 (0.406) showed the strongest stability in all samples ( Figure 3K). In contrast, TUA and PPR were the least stable across most sets.
Best practices include using multiple reference genes as internal controls for standardization to improve the accuracy of RT-qPCR data [28,49]. The number of optimal genes for standardization of the different datasets from L. megaphylla was calculated using the V n/n+1 function of geNorm, with a threshold of 0.15 ( Figure 4). Interestingly, the values of V 2/3 were less than 0.15 for most experimental groups (0.033, 0.045, 0.111, 0.112, 0.126, 0.034, 0.044, 0.052, 0.098, and 0.129) except for the 'all samples' group, as shown in Figure 4. This suggested that two was the optimal number of reference genes for each type of samples. However, for the 'all samples' set, the V 2/3 and V 3/4 values were greater than 0.15, and it is not until four reference genes are used (V 4/5 value of 0.117) that the values is less than 0.15. Thus, at least four genes are required to obtain accurate results across many tissues and treatments.

NormFinder Analysis
To further determine the reliability of the results obtained by the geNorm algorithm, the NormFinder application was used to evaluate both the intra-and inter-group variation to calculate the stability of the candidate reference genes, with lower values of inter-and intra-group variation corresponding to increased stability of the candidate gene [48]. The results of the NormFinder analysis are shown in Figure 5. UBC28 (0.055) and helicase-15 (0.138) were the most stable genes in seedlings, while UBC28 (0.063) and UBC7 (0.123) were the most stable genes in adult trees ( Figure 5A,B). The best combination of reference genes for different tissues was helicase-15 (0.134) and ubiquinone (0.267) ( Figure 5C). In the combined leaf developmental stages set, the optimal combination of candidate reference genes was UBC7 (0.190) and ACT7 (0.245), while over the entire growth cycle, the most stable candidate genes were ubiqunone (0.267) and UBC36 (0.295) ( Figure 5D,E). TCTP (0.030) and ubiqunone (0.035) were the most stable genes under cold treatment for 7 d, while GAPDH (0.062) and UBC36 (0.088) were the most stable genes under cold treatment for 24 h (Figure 5F,G). Including both cold stress treatments, the candidate reference genes EF2 (0.168) and UBC36 (0.193) were the most stable ( Figure 5H). Under heat treatment for 24 h, PAB2 (0.089) and CYP20-2 (0.414) were the most stable reference genes ( Figure 5I). Under all temperature stresses, PAB2 (0.185) and EF2 (0.349) showed the highest stability ( Figure 5J). The best combination of reference genes to compare all sample sets was ubiquinone (0.375), EF2 (0.408), PAB2 (0.413), and GAPDH (0.483) ( Figure 5K). Notably, TUA and PPR were the least stable genes in most sets, similar to the results calculated by geNorm. In general, the results of the NormFinder analysis for the 14 candidate reference genes under different experimental conditions were similar to the results of the ∆Ct and geNorm analyses.

NormFinder Analysis
To further determine the reliability of the results obtained by the geNorm algorithm, the NormFinder application was used to evaluate both the intra-and inter-group variation to calculate the stability of the candidate reference genes, with lower values of inter-and intra-group variation corresponding to increased stability of the candidate gene [48]. The results of the NormFinder analysis are shown in Figure 5. UBC28 (0.055) and helicase-15 (0.138) were the most stable genes in seedlings, while UBC28 (0.063) and UBC7 (0.123) were the most stable genes in adult trees ( Figure 5A, B). The best combination of reference genes for different tissues was helicase-15 (0.134) and ubiquinone (0.267) ( Figure 5C). In the combined leaf developmental stages set, the optimal combination of candidate reference genes was UBC7 (0.190) and ACT7 (0.245), while over the entire growth cycle, the most stable candidate genes were ubiqunone (0.267) and UBC36 (0.295) ( Figure 5D, E). TCTP (0.030) and ubiqunone (0.035) were the most stable genes under cold treatment for 7 d, while GAPDH (0.062) and UBC36 (0.088) were the most stable genes under cold treatment for 24 h (Figure 5F, G). Including both cold stress treatments, the candidate reference genes EF2 (0.168) and UBC36 (0.193) were the most stable ( Figure 5H). Under heat treatment for 24 h, PAB2 (0.089) and CYP20-2 (0.414) were the most stable reference genes (Figure 5I). Under all temperature stresses, PAB2 (0.185) and EF2 (0.349) showed the highest stability ( Figure 5J). The best combination of reference genes to compare all sample sets was ubiquinone (0.375), EF2 (0.408), PAB2 (0.413), and GAPDH (0.483) ( Figure 5K). Notably, TUA and PPR were the least stable genes in most sets, similar to the results calculated by geNorm. In general, the results of the NormFinder analysis for the 14 candidate reference genes under different experimental conditions were similar to the results of the ΔCt and geNorm analyses.

BestKeeper Algorithm
The BestKeeper algorithm evaluates the most stable candidate reference genes based on the standard deviation (SD) and coefficient of variance (CV) of the average cycle threshold (Ct) values [16,19,32]. In general, the smaller the SD value is, the more stable the gene is. As shown in Table 3, the genes UBC7 (0.06 ± 0.24) and helicase-15 (0.14 ± 0.57) had the

Comprehensive Stability Analysis Using RefFinder
Since the different algorithms determined different stability rankings for the candidate genes, the program RefFinder was used to calculate the geometric mean of the ranking results from the four methods [16]. The comprehensively ranked candidate genes (Table 4) did not present one or two universal reference genes for the normalization of gene expression data for all samples. However, the consensus for the top two genes in adult trees during leaf development, under cold treatment for 7 days or 24 h, and under heat treatment for 24 h was consistent with the results of the ∆Ct and NormFinder analyses. The top two reference genes in seedlings or under temperature stress were consistent with the results of the ∆Ct or geNorm analysis, respectively. In the other experimental conditions, the top two most stable genes in the overall ranking appeared in either the top two or top three positions in one of the other four algorithms. In addition, all analyses revealed that TUA was the most unstable gene. We then analyzed the top five most stable genes as determined using the ∆Ct, geNorm, NormFinder, BestKeeper, and RefFinder algorithms ( Figure 6). The top two most stable reference genes selected using RefFinder for the various experimental sets were ranked highly by three or four of the other software programs.
Plants 2022, 11, x FOR PEER REVIEW 20 of 26 Figure 6. Venn diagrams displaying the overlap between the top five most stable reference genes as ranked by ΔCt, geNorm, NormFinder, and BestKeeper for the different datasets. Genes in the overlapping area were confirmed by more than one algorithm. The light purple, purple, light green and pink ellipses contain the top five most stable reference genes selected by ΔCt, geNorm, NormFinder, and BestKeeper, respectively.

Validation of Reference Genes
To verify the reliability of the selected reference genes, the transcripts of two genes were quantified using either different combinations of the two most stable genes, using a single reference gene or using the relatively unstable reference gene TUA under six different experimental conditions. The test genes were NAC and ERF, both of which show relatively high abundance levels, with fragments per kilobase of transcript per million fragments mapped (FPKM) values of approximately 22-223. The relative expression patterns and levels of the LmNAC83 gene showed similar trends when the most stable genes were used alone or in combination as reference genes for standardization ( Figure 7A-F). In contrast, the relative expression levels of LmNAC83 showed significant fluctuations when the relatively unstable gene TUA was used for relative quantification. For example, under heat treatment for 24 h, the expression levels of LmNAC83 in leaves was the highest at 30 °C Figure 6. Venn diagrams displaying the overlap between the top five most stable reference genes as ranked by ∆Ct, geNorm, NormFinder, and BestKeeper for the different datasets. Genes in the overlapping area were confirmed by more than one algorithm. The light purple, purple, light green and pink ellipses contain the top five most stable reference genes selected by ∆Ct, geNorm, NormFinder, and BestKeeper, respectively.

Validation of Reference Genes
To verify the reliability of the selected reference genes, the transcripts of two genes were quantified using either different combinations of the two most stable genes, using a single reference gene or using the relatively unstable reference gene TUA under six different experimental conditions. The test genes were NAC and ERF, both of which show relatively high abundance levels, with fragments per kilobase of transcript per million fragments mapped (FPKM) values of approximately 22-223. The relative expression patterns and levels of the LmNAC83 gene showed similar trends when the most stable genes were used alone or in combination as reference genes for standardization ( Figure 7A-F). In contrast, the relative expression levels of LmNAC83 showed significant fluctuations when the relatively unstable gene TUA was used for relative quantification. For example, under heat treatment for 24 h, the expression levels of LmNAC83 in leaves was the highest at 30 • C when using the most stable genes as the reference genes. However, the relative expression of LmNAC83 was low when using the unstable gene (TUA) as the reference gene ( Figure 7F). In addition, the expression levels and trends of the LmERF60 gene ( Figure 8) were very similar to those found in the analysis of LmNAC83. It is evident that the use of unstable internal reference genes for gene expression analysis in L. megaphylla can lead to unreliable results. This test illustrates the importance of screening stable internal reference genes under different experimental conditions.

Discussion
Changes in plant secondary metabolism are closely related to the transcriptional activities of key genes, and gene expression analysis is a key technique for understanding the mechanisms involved in these processes [32,49]. RT-qPCR is the most accurate technique to obtain gene expression profiles that relate to biological function and regulatory networks [50][51][52]. However, the accuracy of RT-qPCR results mainly depends on using optimal internal reference genes that are stably expressed in samples across different experimental conditions. Reference genes are crucial for normalization of gene expression data and avoiding experimental errors by minimizing non-biological variation between different samples [24,33,35]. To ensure the accuracy of experiments, it is important to select suitable reference genes for each species that similar transcript levels under different experimental conditions. Lindera megaphylla is an ecologically important and dominant broad-leaved evergreen tree species that is naturally distributed in the warm-temperate and subtropical zones of China [1]. This tree contributes to the seasonal landscape and produces volatile compounds with strong effects on bacteria and toxic gases. For example, the terpenes produced have a strong antibacterial effect [6,10,53]. L. megaphylla is also used as a medicinal plant [2]. However, few studies have focused on the molecular biology of L. megaphylla due to limited genomic information, and to date, no reference genes have been reported. In this study, we obtained transcriptome databases of different tissues of L. megaphylla and identified appropriate internal reference genes for use when studying the expression of genes. We initially used 40 internal reference genes from published model plants such as Arabidopsis thaliana and other species in the Lauraceae family to screen for constitutively expressed reference genes from the transcriptome database of L. megaphylla. We obtained 20 candidate genes based on an FPKM value > 50. As shown in Table 2, the final 14 candidate genes were selected based on their E-values in the range of 91.04-107.17% and R2 values in the range of 0.991-0.999, which indicated that the primer pairs for standardized evaluation by RT-qPCR had high sensitivity and accuracy. In addition, the average Ct values of the candidate genes ranged from 17.94 (PAB2) to 27.51 (UBC36), indicating different expression levels ( Figure 1). The results obtained are similar to those of many previous studies, such as on Cryptomeria fortune [33], Gerbera hybrid [19], and Piper species [48]. The results indicate that none of the reference genes had constant expression levels under all tested experimental conditions or in different species. Thus, it is necessary to carefully select the most appropriate reference gene to ensure the accuracy of RT-qPCR analysis for specific conditions and with specific materials. In this study, we combined four statistical algorithms (∆Ct, geNorm, NormFinder, and BestKeeper) to assess the expression stability of 14 candidate genes (TUA, PPR, EIF4A-3, CYP26-2, helicase-15, TCTP, ACT7, PAB2, GAPDH, UBC28, EF2, UBC7, UBC36, and ubiquinone) in different tissues across 16 different leaf developmental stages, and under different temperature stresses. The results demonstrated that the optimal reference genes were not the same under different conditions (Table 4).
There were slight differences in the rankings of candidate reference genes between the different algorithms. However, analysis by the ∆Ct and NormFinder algorithms consistently identified the most stable or unstable candidate reference genes for most experimental sets, while in a few experimental groups, the expression patterns of similar genes were the most or least stable in geNorm and NormFinder. The ranking of candidate genes by BestKeeper suggested some differences compared to the other algorithms. For example, for different leaf developmental stages, the ∆Ct and NormFinder platforms indicated that the ACT7 and UBC36 genes were the most stable. geNorm placed these genes in fourth and first place, while the BestKeeper program placed these genes in sixth and seventh place. Although the rankings of candidate genes produced by the different algorithms were slightly different, the top five stable candidate genes selected by the algorithms were similar for each group of experimental conditions (Table 3). For instance, ACT7, UBC36, TCTP, UBC7, and ubiquinone were the top five most stable genes based on the geNorm and NormFinder analyses across the leaf development stages, and the ∆Ct analysis showed similar results, except for ubiquinone. BestKeeper analysis identified two stable genes: UBC7 and TCTP. Numerous other studies have found similar differences between the outputs of geNorm and NormFinder [32,54], and many studies also demonstrated that these subtle differences result from the use of different algorithm models [33,34,55].
To comprehensively synthesize the results of the four algorithms, RefFinder was utilized to rank the identified candidate genes in L. megaphylla. This analysis plays an important role in integrating the screening results of reference genes from other algorithms by assigning an appropriate weight to each gene and calculating the geometric mean of its weights to produce a final ranking [32,34,35]. Fortunately, we found that the results from RefFinder were similar to those of the different algorithms in each experimental set, proving that RefFinder can assess and screen the optimal reference genes [56], as shown in Figure 6 and Table 4.
The results also indicate that we screened and identified the optimal reference gene combinations for use in L. megaphylla samples generated under different experimental conditions. To compare expression in different seedling tissues, helicase-15 and PAB2 were the most suitable, whereas UBC28 and UBC7 were most stably expressed in different adult trees tissues. When two different groups of tissues were analyzed, helicase-15 and UBC28 emerged as the most stable gene combination. For different leaf developmental stages, ACT7 and UBC36 were best, whereas ubiquinone and UBC7 were best when analyzing samples over the entire growth cycle. Interestingly, for cold stress of 7 days and 24 h, the optimal reference genes were TCTP + ubiquinone and GAPDH + UBC36, respectively. For heat treatment, the best reference genes were PAB2 and CYP20-2, while for overall temperature stress, PAB2 and PPR were most stable. When all samples were tested, ubiquinone, EF2, UBC7, and GAPDH were the optimal candidate reference genes overall for the normalization of gene expression in L. megaphylla (Table 4). These analyses are sufficient to