Comparison of Reliable Reference Genes Following Different Hormone Treatments by Various Algorithms for qRT-PCR Analysis of Metasequoia

Quantitative reverse transcription polymerase chain reaction (qRT-PCR) is the most sensitive technique for evaluating gene expression levels. Choosing appropriate reference genes for normalizing target gene expression is important for verifying expression changes. Metasequoia is a high-quality and economically important wood species. However, few systematic studies have examined reference genes in Metasequoia. Here, the expression stability of 14 candidate reference genes in different tissues and following different hormone treatments were analyzed using six algorithms. Candidate reference genes were used to normalize the expression pattern of FLOWERING LOCUS T and pyrabactin resistance-like 8. Analysis using the GrayNorm algorithm showed that ACT2 (Actin 2), HIS (histone superfamily protein H3) and TATA (TATA binding protein) were stably expressed in different tissues. ACT2, EF1α (elongation factor-1 alpha) and HIS were optimal for leaves treated with the flowering induction hormone solution, while Cpn60β (60-kDa chaperonin β-subunit), GAPDH (glyceraldehyde-3-phosphate dehydrogenase) and HIS were the best reference genes for treated buds. EF1α, HIS and TATA were useful reference genes for accurate normalization in abscisic acid-response signaling. Our results emphasize the importance of validating reference genes for qRT-PCR analysis in Metasequoia. To avoid errors, suitable reference genes should be used for different tissues and hormone treatments to increase normalization accuracy. Our study provides a foundation for reference gene normalization when analyzing gene expression in Metasequoia.


Introduction
Quantitative reverse transcription polymerase chain reaction (qRT-PCR) is an efficient, sensitive and reliable technique for quantifying the expression profiles of target genes in different tissues, following different hormone treatments and under various stresses [1]. qRT-PCR enables comparison of changes in the expression of a target gene to changes in reference genes, including abundant transcripts and low-abundance transcripts of the target gene. Thus, the qRT-PCR is widely preferred over classic transcriptome analysis methods [2]. Appropriate normalization is a key factor in quantifying transcript expression levels to avoid bias [3]. Normalization requires the selection of one or more reference genes showing constant and stable expression levels in different tissues and following different treatments [4,5]. To eliminate differences in the initial volume of cDNA template, quality of RNA from
Only one band of the expected size was detected by 2% agarose gel electrophoresis (Table 1, Figure S1). The specificity of PCR detection was confirmed as a single peak in melting-curve analysis ( Figure S2). Therefore, these primer pairs were used in the next qRT-PCR experiments. The PCR efficiency varied from 1.8518 for RP to 2.1419 for AP-2 (Table 1).
Only one band of the expected size was detected by 2% agarose gel electrophoresis (Table 1, Figure S1). The specificity of PCR detection was confirmed as a single peak in melting-curve analysis ( Figure S2). Therefore, these primer pairs were used in the next qRT-PCR experiments. The PCR efficiency varied from 1.8518 for RP to 2.1419 for AP-2 (Table 1).   The Ct values are listed in Supplementary Table S1 and reveal the different gene expression levels. The minimum Ct value of the 14 reference genes was 17.30 for RA, while the maximum value was 35.87 for TUB. The Ct value of each candidate reference gene is presented in Figure 2. In different Metasequoia tissues and organs, the Ct values varied from 18.80 to 35.87, following flowering induction hormone treatments, the Ct values ranged from 18.98 to 31.99 and the Ct values ranged from 17.36 to 34.77 following ABA treatments.
The expression of 14 genes in all 19 samples was calculated using the 2 ⁻ΔCT method to analyze expression stability; the expression of each gene in the leaves was measured as a control. The heatmap directly reveals the stability of gene expression in all samples ( Figure 3). The expression levels of ACT2, GAPDH, HIS and TATA were stable, while TUB and UBQ showed large variations in expression levels.   The Ct values are listed in Supplementary Table S1 and reveal the different gene expression levels. The minimum Ct value of the 14 reference genes was 17.30 for RA, while the maximum value was 35.87 for TUB. The Ct value of each candidate reference gene is presented in Figure 2. In different Metasequoia tissues and organs, the Ct values varied from 18.80 to 35.87, following flowering induction hormone treatments, the Ct values ranged from 18.98 to 31.99 and the Ct values ranged from 17.36 to 34.77 following ABA treatments.
The expression of 14 genes in all 19 samples was calculated using the 2 ⁻ΔCT method to analyze expression stability; the expression of each gene in the leaves was measured as a control. The heatmap directly reveals the stability of gene expression in all samples ( Figure 3). The expression levels of ACT2, GAPDH, HIS and TATA were stable, while TUB and UBQ showed large variations in expression levels.

Analysis of Gene Expression Stability
The mean standard deviation (mSD) can be used to assess whether the expression of a gene is stable. A low mSD value of a candidate reference gene indicates that the gene is expressed stably in different samples. We used this ∆Ct approach to screen candidate reference genes ( Table 2). The ranking results showed in a box plot visually revealed whether the distribution of ∆Ct values between a gene and the other 13 genes was discrete or concentrated. A longer box indicates greater deviation in the ∆Ct value, suggesting that one or both genes is changing in the sample and that expression is relatively unstable. In contrast, a shorter box indicates smaller deviations, suggesting relatively stable expression of two genes ( Figure 4).  In BestKeeper analysis, the correlation coefficient (r), standard deviation (SD) and coefficient of variation (CV) indicate the variation among Ct values, which can determine candidate gene stability. In this study, the most stable gene in different tissues was AP-2, while the least stable gene was RA (Table 3). HIS was selected as the most stable gene in leaves following MglFlora treatment and Cpn60β was the most stable in the treated buds, while the least stable gene was AP-2 (Table 3). In contrast to In BestKeeper analysis, the correlation coefficient (r), standard deviation (SD) and coefficient of variation (CV) indicate the variation among Ct values, which can determine candidate gene stability. In this study, the most stable gene in different tissues was AP-2, while the least stable gene was RA (Table 3). HIS was selected as the most stable gene in leaves following MglFlora treatment and Cpn60β was the most stable in the treated buds, while the least stable gene was AP-2 (Table 3). In contrast to the ranking order identified by geNorm, TATA was the best choice for samples treated with ABA, while TUB is the most unstable gene. geNorm uses the mean expression stability values (M-values) to determine the reference gene expression stability, which is represented as the average pairwise variation (PV, V) of a specific gene with all other tested candidate genes. Lower M-values indicate greater expression stably, while higher M-values indicate lower expression stability. In our analysis, the most stable genes were AP-2 and Cpn60β with an M-value of 0.478 in different tissues (Table 4; Figure 5A). RPL17 and elF-5A were the most suitable candidate genes for normalization with an M-value of 0.105 in the leaves following MgFlora treatment (Table 4; Figure 5B), whereas EF1α and HIS were the most stable genes with an M-value of 0.19 in treated buds (Table 4; Figure 5C). For ABA treatment, EF1α and TATA showed stable expression with an M-value of 0.139 (Table 4; Figure 5D). The geNorm algorithm revealed the appropriate number of reference genes needed for standardization to calculate the pairwise variation between the normalization factors NF n and NF n+1 of two sequences (V n/n+1 ). When V n/n+1 is less than 0.15, n is the number of genes used in normalization.
In this study, all V 2/3 values were less than 0.15 in the four sample sets and thus the best number for normalization was two ( Figure 5E-H). Considering the accuracy of the experiment and previous reports, three reference genes were used in this study for target gene expression normalization in Metasequoia.  According to the NormFinder algorithm, the most stable candidate reference genes have low stability values (SV). When the stability value is close to zero, the gene displays stable expression. The NormFinder rankings of the 14 candidate genes are shown in Table 5. In the six different tissues, AP-2 ranked first with an SV of 0.197, while TUB ranked last with an SV of 3.839. TATA was found to be the most stable choice for the leaf and bud of with SVs of 0.152 and 0.233, whereas ACT2 and AP-2 were the most unstable genes in flowering induction hormone treatments, respectively. TATA was the most stable gene in response to ABA treatment with an SV was 0.07, which contrasted the results calculated by the NormFinder algorithm, showing that the gene with the lowest ranking was TUB with an SV of 4.255. Different algorithms yielded different ranking results in this study. RankAggreg was used to determine the final comprehensive ranking list of all genes based on the mSD value of ∆Ct, the SD of Bestkeeper, the SV of NormFinder and the M value of geNorm. The results shown as line graphs suggest that the most stable gene combination in different tissues was AP-2, Cpn60β and elF-5A ( Figure 6A). In the leaf treated with MglFlora, TATA, elF-5A and RPL17 were the most suitable choices for normalization ( Figure 6B), while TATA, HIS and ACT2 ranked highest in the treated bud ( Figure 6C). Additionally, EF1α, elF-5A and TATA were more stable compared to other genes following ABA treatment ( Figure 6D).
GrayNorm is a new algorithm useful for identifying combinations of reference genes. This program provides all combinations of all genes; the number of combinations is reflected as 2 n−1 (where n is the number of candidate genes). According the lowest number of reference genes required for reliable normalization provided by geNorm, four gene combinations in this study were selected based on the GrayNorm results ( Table 6, Table S2). Table 6 shows that a combination of ACT2, HIS and TATA should be used in different tissues. A combination ACT2, EFα and HIS can be used in leaves treated by flowering induction hormone solution, while a combination of Cpn60β, GAPDH and HIS should be used in treated buds. EF1α, HIS and TATA can be used to evaluate ABA-response signaling.

Validation of selected candidate reference genes
The relative expression data under MglFlora treatment are shown in Figure 7. The expression of MgFT was not increased sharply until 1 and 7 days after treatment but was significantly changed by 14 days of treatment ( Figure 7A). The MgFT expression pattern in treated buds was similar to that in the leaf, with expression after 1 day of treatment showing slight upregulation. In the 7 day bud sample set, expression was the same as that on day 1. However, the expression of MgFT in the bud was sharply increased after 14 days of treatment and was significantly higher than that under the same conditions used in the treated leaf ( Figure 7B). MgPYL8 expression is shown in Figure 7C. When the expression levels of these two genes were normalized to the different reference gene combinations recommended by the different algorithms, there was a large change in the maximum value of the gene expression and fold-change. Comparing normalized data to non-normalized data, a smaller difference leads to more accurate results. The expression levels of the two genes normalized using the combination recommended by six algorithms showed a small fold-change and the expression level was similar to that of non-normalized genes.

Validation of Selected Candidate Reference Genes
The relative expression data under MglFlora treatment are shown in Figure 7. The expression of MgFT was not increased sharply until 1 and 7 days after treatment but was significantly changed by 14 days of treatment ( Figure 7A). The MgFT expression pattern in treated buds was similar to that in the leaf, with expression after 1 day of treatment showing slight upregulation. In the 7 day bud sample set, expression was the same as that on day 1. However, the expression of MgFT in the bud was sharply increased after 14 days of treatment and was significantly higher than that under the same conditions used in the treated leaf ( Figure 7B). MgPYL8 expression is shown in Figure 7C. When the expression levels of these two genes were normalized to the different reference gene combinations recommended by the different algorithms, there was a large change in the maximum value of the gene expression and fold-change. Comparing normalized data to non-normalized data, a smaller difference leads to more accurate results. The expression levels of the two genes normalized using the combination recommended by six algorithms showed a small fold-change and the expression level was similar to that of non-normalized genes.  Error bars indicate standard errors of three replicates (bars ± SE).

Discussion
The selection of stably expressed genes in different tissues and treatments is a key step in qRT-PCR analysis. However, in gymnosperms, the selection of reference genes is limited to Cycas elongate [37], Abies alba [38] and marine pine [39]. Metasequoia originated in the Mesozoic Cretaceous and was declared extinct on Earth until 1941, when Chinese botanists discovered traces of this ancient and rare species in Hubei province, which had been preserved during the Quaternary Ice Age. Metasequoia is referred to as a "living fossil" and is a unique tree species in China that has not evolved and its growth process has not diverged. Few studies have examined reference genes in Metasequoia. Therefore, we evaluated the stability of the expression of candidate reference genes under specific conditions and then analyzed the usefulness of these genes as internal controls to normalize the expression of other genes.
In our study, 14 candidate reference genes were successfully identified and their expression characteristics were further analyzed using six statistical algorithms: ΔCt, Bestkeeper, NormFinder, geNorm, RankAggreg and GrayNorm. As shown in Figure 1 and Figure 2, each candidate reference gene showed different Ct values in all samples. The stability of each gene in different samples was clearly observed. ACT2 and TATA exhibited the most stable expression pattern. Our data demonstrate that TATA had the highest rank in six tissues and following different hormone treatments. The TATA-box binding protein interacts with the TATA-box, which was the first promoter identified in eukaryotes. The TATA-box is relatively fixed in eukaryotic promoters and plays an important role in regulating the initiation of gene transcription [40]. Notably, the TATA-box binding protein is preferable to other commonly used genes for qRT-PCR analysis in various species, such as, mouse [41], human [42] and Aphis glycines [43]. Recently, Izabela also recommended using this gene for normalization in molecular studies of primary and secondary dormancy in Avena fatua L [44]. An important feature of an ideal reference gene is that its expression should be constant

Discussion
The selection of stably expressed genes in different tissues and treatments is a key step in qRT-PCR analysis. However, in gymnosperms, the selection of reference genes is limited to Cycas elongate [37], Abies alba [38] and marine pine [39]. Metasequoia originated in the Mesozoic Cretaceous and was declared extinct on Earth until 1941, when Chinese botanists discovered traces of this ancient and rare species in Hubei province, which had been preserved during the Quaternary Ice Age. Metasequoia is referred to as a "living fossil" and is a unique tree species in China that has not evolved and its growth process has not diverged. Few studies have examined reference genes in Metasequoia. Therefore, we evaluated the stability of the expression of candidate reference genes under specific conditions and then analyzed the usefulness of these genes as internal controls to normalize the expression of other genes.
In our study, 14 candidate reference genes were successfully identified and their expression characteristics were further analyzed using six statistical algorithms: ∆Ct, Bestkeeper, NormFinder, geNorm, RankAggreg and GrayNorm. As shown in Figures 1 and 2, each candidate reference gene showed different Ct values in all samples. The stability of each gene in different samples was clearly observed. ACT2 and TATA exhibited the most stable expression pattern. Our data demonstrate that TATA had the highest rank in six tissues and following different hormone treatments. The TATA-box binding protein interacts with the TATA-box, which was the first promoter identified in eukaryotes. The TATA-box is relatively fixed in eukaryotic promoters and plays an important role in regulating the initiation of gene transcription [40]. Notably, the TATA-box binding protein is preferable to other commonly used genes for qRT-PCR analysis in various species, such as, mouse [41], human [42] and Aphis glycines [43]. Recently, Izabela also recommended using this gene for normalization in molecular studies of primary and secondary dormancy in Avena fatua L [44]. An important feature of an ideal reference gene is that its expression should be constant regardless of the environmental and experimental conditions. However, recent studies showed that widely used reference genes are only stably differentially expressed under certain conditions. TUB was predicted to be an unstable reference gene in this study. TUB is a member of the tubulin gene family, which are commonly used as reference genes [45] but in Metasequoia, its expression stability was low in different tissues and hormone treatment. TUB showed lower stable expression than several other reference genes under abiotic stresses in Seashore Paspalum [46]. The results show that not all common reference genes are suitable as reference genes without selection and validation and the same reference genes may have different expression characteristics in different tissue samples and experimental conditions. Notably, several Metasequoia reference genes exhibited different expression patterns compared to in other species treated with ABA. It was previously reported that UBQ is the most stable reference gene in NaCl-treated and ABA-treated Platycladus orientalis [45]; however, UBQ exhibited unstable expression following ABA treatment in this study. Several other widely used reference genes, including EF1α, elF-5A and TATA, were relatively better reference genes for ABA. In previous studies, TATA was used as a stable reference gene under heavy metal stress [47] but few reports have examined its expression under abiotic stress. EF1α is a commonly used reference gene and was shown to be one of the best reference genes for ABA in Populus euphraticarice [8]. Additionally, EF1α was selected as best choice for Petunia hybrida [48], soybean [49] and B. brizantha [50]. This indicates that the reference gene verified in Metasequoia is identical to previously reported reference genes. We focused on the influence of hormones on gene expression patterns. Evaluation of hormone treatment revealed the large impact of ABA on gene expression (Table 5, Figure 5D,H). ABA treatment is a type of abiotic stress that leads to activation of numerous stress-related genes and the synthesis of multiple functional proteins in plants. The expression patterns of these reference genes are likely to be greatly affected during this process. TATA and elF-5A are expressed stably during the flowering stage and are induced by hormone treatment in the leaves and in ABA-treated leaves, suggesting that the expression of these two genes in leaves is not affected by hormones, although different types of hormones have little effects on their expression in the leaves. In future studies, TATA and elF-5A can be used as valuable reference genes in qRT-PCR assays of hormone treatment in Metasequoia. Additionally, HIS and ACT2, which are also stably expressed in treated buds, did not show stable expression in the treated leaves. These results demonstrate that different tissues under the same treatment conditions exhibit variable gene expression patterns. Our research highlights the importance of choosing suitable reference genes for different tissues.
The expression levels of the target genes were significantly changed when different reference genes were used for normalization, leading to unreliable experimental results. To identify stably expressed reference genes and compare the advantages of these different algorithms, we analyzed and compared the relative expression levels of two functional genes following different hormone treatments. The results revealed numerous differences in the calculated expression patterns of these two target genes when different combinations of reference genes were used. Our results in Figure 7 support that a combination of reference genes is preferred rather than a single internal reference gene. The results of this study provide important information on ABA-response signaling in gymnosperms. Gymnosperms have large genomes but their genome-wide information is limited and the current results provide appropriate resources for qRT-PCR analysis in other species related to Metasequoia.

Hormone Treatment and Sample Collection
Metasequoia trees located on the campus of Beijing Forestry University (latitude 40.0 N, longitude 116.2 E, 57 m above sea level) were growing in the natural environment. Different tissues were collected at 17:00, including leaves, stems, roots, buds, male cores and female cores. After collections, all samples were immediately frozen in liquid nitrogen and stored at −80 • C until RNA extraction.
Metasequoia flowering induction hormone solution contained indole butyric acid, Zeatin ribonucleoside (Sigma-Aldrich, St. Louis, MO, USA) and 0.5% Tween-20 [51]. Current-year leaves and buds were treated with hormone solution, which we named as MglFlora. The treatment was conducted on three Metasequoia trees and each represented a biological replicate. The treated shoots were evenly sprayed with MglFlora until liquid drops began to drip down from the leaves at 17:00 on 11 July 2017. In this treatment, mature leaves and buds from the same position were collected at 1, 7 and 14 days after treatment. For ABA treatment, the leaves were sprayed with 200 µM ABA (Sigma-Aldrich, St. Louis, MO, USA) solution at 17:00 on 11 July 2017. Leaves from the same position were collected at 0, 0.5, 1, 2, 4, 8 and 12 h after treatment. After collection, all samples were immediately frozen in liquid nitrogen and stored at −80 • C until RNA extraction. The average temperature in Beijing Forestry University from July 21 to 25, 2017 is 28.28 • C; the average rainfall during this period is about 170 mm.

Selection of Candidate Reference Genes and Primer Design
In a previous study, the first large-scale dataset of expressed sequence tags (ESTs) of Metasequoia from vegetative and reproductive tissues was reported using 454 pyrosequencing technology [52]. Based on the ESTs, 14 putative reference genes were selected as candidates for qRT-PCR normalization and two target genes were selected to investigate the effect of the choice of reference genes on normalization.

Total RNA Isolation and cDNA Synthesis
Frozen samples were ground into a fine powder under liquid nitrogen with a pestle and mortar. Total RNA was isolated from all samples using cetyltrimethyl ammonium bromide (CTAB) buffer according to the instructions [53]. Dithiothreitol (DTT, Thermo Fisher Scientific, Waltham, MA, USA) is a strong denaturant that impacts RNase A (Thermo Fisher Scientific, Waltham, MA, USA) enzyme activity and free thiol, thereby reducing the RNase A enzyme content. DTT was added to the CTAB buffer to improve the extraction efficiency. After chloroform/isoamyl alcohol extraction, RNA was precipitated with 10 M LiCl 2 , followed by extraction with SSTE buffer (1 M NaCl, 0.5% SDS, 10 mM Tris-HCl, 1 mM EDTA) and precipitated with ethanol again. Finally, 75% alcohol were used to clean the RNA precipitate. After the remaining alcohol was evaporated, all RNA samples were dissolved in RNA-free water for subsequent cDNA synthesis experiments. RNA quantity and purity were confirmed by measuring the optical density at OD260/OD280 and OD260/OD230 absorption ratio using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The integrity of the total RNA was checked on 0.8% agarose gels. After adjusting each of the RNA samples to the same concentration, 1000 ng RNA samples were prepared for the synthesis of cDNA strands. cDNA synthesis was performed using PrimerScript RT Enzyme (TaKaRa, Shiga, Japan) following the manufacturer's instructions.

qRT-PCR Conditions and Amplification Efficiency
qRT-PCR was performed on a Roche 480 light Cycler instrument (Roche, Basel, Switzerland) using TB Green ™ Premix Ex Taq ™ II (Tli RNaseH Plus) (TaKaRa, Shiga, Japan). The reaction mixture was a total of 20 µL and contained 1 µL cDNA, 10 µL 2× TB Green ™ Premix Ex Taq ™ II (Tli RNaseH Plus), 0.6 µL 10 µM forward primer, 0.6 µL 10 µM reverse primer and 7.8 µL ddH 2 O. The amplification conditions were as follows: 95 • C for 10 s; 45 cycles at 95 • C for 5 s; 60 • C for 30 s; and 72 • C for 30 s. Each assay included three technical replicates. A standard curve was used to calculate the PCR efficiency (E) for each gene by qRT-PCR using a five-fold dilution series of the mixed cDNA template. The PCR efficiency was calculated by the equation: E = 10 −1/slope , where slope indicates the slope of a standard curve was given by LightCycler software (Roche, Basel, Switzerland).
The principle of the ∆Ct algorithm is to compare the relative expression of the 'gene pair' in each sample to identify the reference gene that can be stably expressed. The difference in Ct (∆Ct) between the two genes is first calculated. If the ∆Ct of the two genes in different cDNA samples remains the same, these two genes were considered as stably expressed. The standard deviation (sd) reflects the change in the ∆Ct of the two genes in all samples. Next, the ∆Ct value between one gene and the other 13 genes was calculated and the 13 SD values were averaged to obtain the mean standard deviation (mSD).
BestKeeper is a tool that determines the stable level of reference genes based on pairwise correlations calculated from the percentage standard deviation (SD) and covariance (CV) among the raw Ct value data. Candidate reference genes showing the lowest CV values were considered as the most stable genes. Another function of the BestKeeper program is that it applies the decision coefficient (r) as an index of the program's stability characterization. r was calculated to determine a credible normalization factor (NF), rather than evaluating the excellence of each reference gene's expression stability separately. In this study, we used r to rank the stability values of the 14 candidate genes.
NormFinder values the gene stability based on the expression variations of these candidate reference genes, with a lower stability value indicating more stable reference genes. The algorithm calculates inter-group and intra-group changes in response to the NF, with the algorithm requiring at least three candidate genes and more than two samples per group.
geNorm analyzes the expression stability of a candidate reference gene according the stability value (M, M-value), with smaller M values indicating better stability. The standard of selected reference genes was based on the M value, with values of less than 1.5 considered to indicate stable expression of an ideal reference gene. Additionally, geNorm also valued the pairwise variation values (V), indicating the lowest number of reference genes required for reliable normalization.
GrayNorm is an algorithm that can identify a combination of reference genes with minimal deviation from non-normalized data. The principle of GrayNorm is to calculate the NF for each combination of treatment group and each possible reference gene. The closer the mean value of 1/NF for each treatment group is to 1.0, the smaller the biological variation and the more accurate the expression level of gene of interest can be calculated.
The R program v3.0.1 can load the RankAggreg v0.6.4 package. RankAggreg meets the needs of complex rank aggregation easily and expediently, integrating the stability measurements obtained from the four methods and then establishes a comprehensive ranking of reference genes. RankAggreg contains various algorithms, such as Cross-Entropy Monte Carlo algorithm, Genetic algorithm and a brute force algorithm. The detailed algorithm process can be referenced at https://cran.r-project.org/ web/packages/RankAggreg/RankAggreg.pdf. According to the ranking orders of these four methods, the Cross-Entropy Monte Carlo algorithm (CE) was used to demonstrate hierarchical aggregation. The result was intuitively output in the form of a line graph.

Conclusions
Stable reference genes are used to determine the relative gene expression levels and show stable expression under different treatment conditions for qRT-PCR. Genes utilized as qRT-PCR candidate reference genes have not been reported for Metasequoia. This study was conducted to select reference genes after evaluating their expression stability in various tissues and hormone conditions. The results of this study emphasize the importance of validating reference genes for qRT-PCR analysis in Metasequoia.