Selection of Suitable Reference Genes Based on Transcriptomic Data in Ginkgo biloba under Di ﬀ erent Experimental Conditions

: Ginkgo biloba , a deciduous tree species in the Ginkgo family, has a long history of cultivation in China and is widely used in garden landscapes, medicine, food, and health products. However, few reports have focused on the systematic selection of optimal reference genes based on transcriptomic data in G. biloba . The purpose of our research was to select an internal reference gene suitable for di ﬀ erent experimental conditions from thirteen candidate reference genes by the delta cycle threshold ( ∆ Ct) method, geNorm, BestKeeper, NormFinder, and RefFinder programs. The reference genes were used for gene expression analyses of Ginkgo biloba . These results showed that elongation factor 1( EF1 ) and ubiquitin ( UBI) were the best choices for samples of di ﬀ erent ginkgo genotypes. The expression of UBI and HAS28 presented the most stable at di ﬀ erent developmental stages of ginkgo, and EIF3I and RPII were considered as suitable reference genes in di ﬀ erent tissues of ginkgo. For methyl jasmonate (MeJA) treatment, ACA and ACT were identiﬁed as the optimal reference genes. For cold stress treatment, RPII and EIF4E were chosen for the gene expression normalizations. HAS28 and GAPDH presented the most stable expression for the heat treatment. To validate the above results, a chalcone synthase gene ( GbCHS ) in ginkgo was ampliﬁed by quantitative real-time polymerase chain reaction (qRT-PCR). Our results provide di ﬀ erent suitable reference genes for further gene expression studies in ginkgo.


Introduction
Quantitative real-time polymerase chain reaction (qRT-PCR) is a commonly used method to analyze the expression level of functional genes and has become the preferred method to quantify mRNA, as this method meets the requirements for quantitative analyses of data in molecular medical, biotechnological, microbiological, and diagnostic-type studies [1,2]. The benefits of this technique compared with traditional RNA measurement methods include its good specificity, high accuracy, wide detection range, simple operation, and safety. Some traditional reference genes (ribosomal 18S, glyceraldehyde-3-phosphate dehydrogenase, and elongation factor 1-α) have been used extensively for RT-qPCR-based analysis [3][4][5][6]. However, the most appropriate reference genes may not be the same for different experimental groups, as shown for Elymus sibiricus [7], Glycine max [8], Fragaria ananassa [9], Caragana korshinskii [10], and Daucus carota [11]. Inappropriate selection of internal reference genes may lead to quantitative errors. Therefore, the expression data may be misinterpreted. In the last decade, several methods, including the delta cycle threshold (∆Ct) method [12] and the use of software programs such as geNorm [13], BestKeeper [14], NormFinder [15], and RefFinder [16], have been commonly used to select the appropriate reference genes.

Materials and Treatments
Leaf samples (4)(5)(6) of six different ginkgo genotypes were collected from the Pizhou Ginkgo Germplasm Garden and Nanjing Forestry University (Supplementary Table S1). Leaf samples of genotype Nanjing (NJ) were collected every 30 days from May to October 2018 (M1-M6). Different tissue samples were collected as follows: branches, male flowers, and female flowers were collected at the flowering stage, and leaves and testae were collected during the fruiting period. Ginkgo seedlings of uniform size were subjected to hormone and abiotic stress treatments. For the methyl jasmonate (MeJA) treatments, the leaves of seedlings were sprayed with 1 mmol/L MeJA solution and were collected after 0, 12, and 24 h. For the heat and cold treatments, the seedlings were subjected to 42 and 4 • C for 0, 12, and 24 h, respectively. Under each treatment, 4-6 leaves were collected at each time point. A total of 87 samples (3 biological replicates each) were collected for this study. All the samples were collected and quickly frozen in liquid nitrogen and then stored at −80 • C for later use.

Total RNA Isolation and cDNA Synthesis
The total RNA of ginkgo samples was isolated by an RNA Simple Total RNA Kit (Omega Bio-tek, Inc., Norcross, GA, USA). RNA quality and concentration were determined by a NanoDrop 2000 spectrophotometer (Thermo, Wilmington, DE, USA). The OD260/280 nm and OD260/230 absorbance ratios of the total RNA were 1.8-2.0 and 2.0-2.6, respectively. The integrity of the RNA was evaluated by 1% agarose gel electrophoresis. Then, following the manufacturer's instructions, cDNA was synthesized in 20 µL of the reaction mixture with the PrimeScript™ RT master mix (TaKaRa, Beijing, China), and cDNA was stored at −20 • C until use.

Gene Sequence Search and Primer Design
Thirteen candidate reference genes and one target gene were selected based on the transcriptomic data of eight different tissues (root, immature leaf, mature leaf, microstrobilus, ovulate strobilus, immature fruit, stem, and mature fruit tissues) and six different leaf developmental stages in ginkgo (unpublished data) (Supplementary Table S2). All the RNA sequencing (RNA-seq) data of eight tissues were downloaded from the National Center for Biotechnology Information (NCBI) database (accession Nos. SRR7948405-SRR7948413 and SRP149113) [40]. Fragments per kilobase of exon model per million mapped reads (FPKM), the mean expression values (MV), standard deviations (SDs), and coefficients of variation (CVs) (calculated by the division of the SDs by the MVs), were calculated to estimate the expression stability of each gene, and stable genes with a CV < 0.5 were chosen (Supplementary Table S2) [41][42][43]. The specificity of the primers was confirmed by 1% agarose gel electrophoresis. The primers were designed with Primer Premier 5.0, and the specific sequences and related parameters are shown in Table 1.

qRT-PCR Assays
qRT-PCR amplification was performed with TB Green ® Fast qPCR Mix (TaKaRa, Beijing, China) on an Applied Biosystems 7300 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). All cDNA was diluted to 100 ng·µL −1 . The PCR mixtures and the procedures were set according to the instructions. The primer specificity was verified by evaluating the melting curve of each reaction. We used a chalcone synthase gene (GbCHS) to validate the accuracy of the final selected internal reference genes. All qRT-PCR experiments' reactions contained three technical replicates and three biological replicates.

Data Analysis
The ranking of the thirteen reference genes was calculated by the geNorm, NormFinder, BestKeeper, the ∆Ct method, and RefFinder (http://www.ciidirsinaloa.com.mx/RefFinder-master/). The stability value (M) and the pairwise variation (Vn/Vn + 1) are the main parameters used for evaluating the stability level of the reference gene within the geNorm [13]. A gene associated with the lower M value is more stable. The NormFinder program also screens the most suitable internal reference genes based on stable values, not only to rank the intragroup expression of the reference gene but also to calculate intergroup variations [15]. Through the BestKeeper, the standard deviation (SD) and coefficient of variation (CV) were calculated, and the reference gene having the lower CV and SD (CV ± SD) is considered the more stable [14]. RefFinder is a web-based comprehensive tool that integrates the above computational programs. Therefore, RefFinder was used for the overall final ranking of the tested genes [16]. The correlation coefficient (R 2 ) and the amplification efficiency (E) were calculated by the slope (k) of standard curves. The applied equation was E = (10 −1/k −1) × 100% [44]. The regression coefficient (R 2 ) of the standard curve was greater than 0.98, and the E value was between 90% and 110%. The expression profiles of GbCHS were validated using the 2 −∆∆ct method with three biological samples [45], and the significant differences were analyzed with SPSS statistical software 23.

Designing and Validation of Primers
Thirteen potential reference genes (Table 1) based on the transcriptomic data from ginkgo were selected, and the sequences were used as the basis for specific primer design. Primer specificity was evaluated by RT-PCR products. PCR products generated by each primer pair using cDNA from G. biloba as a template were visualized as single bands on 1% agarose gel. All amplification products were between 100 and 200 bp (Supplementary Figure S1). Additionally, the no-template controls showed no peaks, which indicated that there were no primer dimers present. The qRT-PCR E-values of the thirteen candidate reference genes ranged from 95.76% (HAS28) to 107.28% (ADSS), and the R 2 ranged from 0.9801 (RPII) to 0.999 (ACT).

Expression Profile of Candidate Reference Genes
For the transcription profile analysis, the primers of thirteen candidate reference genes were used for qRT-PCR to obtain the Ct values of all samples. The Ct values of thirteen reference genes ranged from 20.82 (H2A expressed in the YL samples) to 28.70 (GAPDH expressed in the testa samples) among all the samples (Figure 1).

geNorm Analysis
The geNorm algorithm was used to evaluate the stability values (Ms) of reference genes and calculate the average pairwise variations (Vs) to determine the optimal number of internal reference genes. According to the geNorm program, only candidate genes with an M value lower than 1.5 could be selected as suitable reference genes. The smaller the M value is, the higher the expression stability of the candidate reference gene. Most genes were stable, with M values lower than 1.5 ( Figure 2). Our results showed that RPII and UBI had the smallest M values (0.054), so these two genes were the most stable in the different genotypes, while HAS28 was the least stable (0.625) (Figure 2a). At the different developmental stages, H2A and UBI with an M value of 0.107 were the most suitable genes, while ADSS with the highest M value (0.621) was the most unsuitable ( Figure 2b). EIF3I and HAS28 were indicated as the most suitable genes in the different tissues, but the common internal reference gene GAPDH was considered as the least suitable ( Figure 2c). ACT and ACA (0.099) were identified as the most stably expressed reference genes for MeJA treatment (Figure 2d). After heat stress treatment, GAPDH and HAS28 (0.04) were the most suitable genes with a stability value of 0.04 (Figure 2e). For the above treatments, EF1 was the most unstable gene. Within the cold stress group, EIF4E and RPII (M = 0.097) were the most suitable genes, and H2A (0.557) ranked last (Figure 2f). ACT and HAS28, reflecting the highest stability, were the optimal reference genes, while GAPDH was also the worst choice in all the samples (Figure 2g). On the basis of our results, we found that the most appropriate genes under various experimental conditions were different.

geNorm Analysis
The geNorm algorithm was used to evaluate the stability values (Ms) of reference genes and calculate the average pairwise variations (Vs) to determine the optimal number of internal reference genes. According to the geNorm program, only candidate genes with an M value lower than 1.5 could be selected as suitable reference genes. The smaller the M value is, the higher the expression stability of the candidate reference gene. Most genes were stable, with M values lower than 1.5 ( Figure 2). Our results showed that RPII and UBI had the smallest M values (0.054), so these two genes were the most stable in the different genotypes, while HAS28 was the least stable (0.625) (Figure 2a). At the different developmental stages, H2A and UBI with an M value of 0.107 were the most suitable genes, while ADSS with the highest M value (0.621) was the most unsuitable ( Figure 2b). EIF3I and HAS28 were indicated as the most suitable genes in the different tissues, but the common internal reference gene GAPDH was considered as the least suitable (Figure 2c). ACT and ACA (0.099) were identified as the most stably expressed reference genes for MeJA treatment (Figure 2d). After heat stress treatment, GAPDH and HAS28 (0.04) were the most suitable genes with a stability value of 0.04 (Figure 2e). For the above treatments, EF1 was the most unstable gene. Within the cold stress group, EIF4E and RPII (M = 0.097) were the most suitable genes, and H2A (0.557) ranked last (Figure 2f). ACT and HAS28, reflecting the highest stability, were the optimal reference genes, while GAPDH was also the worst choice in all the samples (Figure 2g). On the basis of our results, we found that the most appropriate genes under various experimental conditions were different. The optimal number of reference genes to be used was determined by pairwise variation (V), which is an important parameter in the geNorm and the cut-off value was 0.150 ( Figure 3). As shown in Figure 3, all Vs were lower than 0.15 in the groups, except for V12/13 of the different tissues (0.185) and the cold stress group V12/13 (0.204), so two reference genes were adequate for normalization. In all the samples, V5/6 (0.137) and V4/5 (0.177) indicated that ACT, HAS28, HYP, ADSS, and TUB were needed for normalization. The optimal number of reference genes to be used was determined by pairwise variation (V), which is an important parameter in the geNorm and the cut-off value was 0.150 ( Figure 3). As shown in Figure 3, all Vs were lower than 0.15 in the groups, except for V12/13 of the different tissues (0.185) and the cold stress group V12/13 (0.204), so two reference genes were adequate for normalization. In all the samples, V5/6 (0.137) and V4/5 (0.177) indicated that ACT, HAS28, HYP, ADSS, and TUB were needed for normalization. . Pairwise variation (V) of thirteen candidate reference genes according to geNorm. The Vn/n + 1 is used to determine the optimal number of reference genes for accurate normalization. A cut-off value of Vn/n + 1 < 0.15 indicates that no additional reference genes need to be added. Otherwise, n + 1 genes need to be introduced.

NormFinder Analysis
In NormFinder analysis, reference genes were ranked according to their stability value (Figure  4a and Supplementary Table S4). The gene with the lowest stability value ranked first. The results indicated that EF1 was the optimal reference gene for the different genotypes, but was the most inappropriate choice for different developmental stages. For different tissues, EIF3I was the most stable gene, with a stability value of 0.248. The stability value of HAS28 was the most stable for different developmental stages and heat treatments. For cold stress, EIF4E was the most stable gene while H2A was the least stable gene. The most suitable reference genes indicated by the NormFinder algorithm in most experimental datasets differed from those obtained with the geNorm program. However, the results of the two methods for the most unstable reference gene in all groups were consistent. For the MeJA treatment, RPII was the most suitable gene within NormFinder analysis, whereas RPII was ranked seventh by geNorm. In all the samples, HYP was the most stable gene, but HYP was ranked third with the geNorm program. Pairwise variation (V) of thirteen candidate reference genes according to geNorm. The Vn/n + 1 is used to determine the optimal number of reference genes for accurate normalization. A cut-off value of Vn/n + 1 < 0.15 indicates that no additional reference genes need to be added. Otherwise, n + 1 genes need to be introduced.

NormFinder Analysis
In NormFinder analysis, reference genes were ranked according to their stability value (Figure 4a  and Supplementary Table S4). The gene with the lowest stability value ranked first. The results indicated that EF1 was the optimal reference gene for the different genotypes, but was the most inappropriate choice for different developmental stages. For different tissues, EIF3I was the most stable gene, with a stability value of 0.248. The stability value of HAS28 was the most stable for different developmental stages and heat treatments. For cold stress, EIF4E was the most stable gene while H2A was the least stable gene. The most suitable reference genes indicated by the NormFinder algorithm in most experimental datasets differed from those obtained with the geNorm program. However, the results of the two methods for the most unstable reference gene in all groups were consistent. For the MeJA treatment, RPII was the most suitable gene within NormFinder analysis, whereas RPII was ranked seventh by geNorm. In all the samples, HYP was the most stable gene, but HYP was ranked third with the geNorm program.

BestKeeper Analysis
The BestKeeper program ranks reference genes according to SD, CV, and r values. A lower SD value corresponded to a more stable gene (Figure 4b and Supplementary Table S4). Therefore, the reference gene EF1 was the most suitable in the different genotypes. Across the different developmental stages, the expression of H2A was the most stable, so H2A was considered as the most appropriate reference gene. For the different tissues, EIF3I was the optimal reference gene with the lowest SD value of 0.23. In the MeJA treatment group, EIF3I was also the most suitable gene according to BestKeeper. For cold stress, GAPDH was considered to be the most stable gene with the lowest SD value (0.17). For heat stress, the SD value of ACA was the lowest, so ACA was the most stable gene. In all the samples, HAS28 was the best choice for normalization of target gene expression. RPII, EIF4E, H2A, and GAPDH were considered to be the least stable genes, as their SD values were the highest.

Analysis by the ∆Ct Method
The value of the standard deviation (SD) was used as the indicator for evaluating the expression stability of reference genes. A reference gene with a lower SD value is more stable (Figure 4c and Supplementary Table S5). In the groups of different genotypes, developmental stages, and tissues, UBI, HYP, and EIF3I presented the most stable expression, while HAS28, ADSS, and GAPDH presented the least stable expression. In the MeJA, cold stress, and heat stress treatment groups, ACA, UBI, and HAS28 were the most appropriate reference genes. EF1 was the worst in both the MeJA treatment and heat stress groups. For cold stress, H2A was the worst choice. These results were the same as the results obtained with the geNorm, NormFinder, and BestKeeper analyses.

RefFinder Analysis
RefFinder performed a comprehensive evaluation to obtain the overall final ranking based on the ranking of the other four methods ( Table 2). Among the different genotypes, UBI, EF1, and RPII were the three most stable genes, while HYP, HAS28, and UBI ranked in the top three across the different developmental stages. In the different issues, EIF3I, RPII, and HAS28 were the three most suitable choices. Under MeJA treatment, ACA, ACT, and HAS28 ranked first to third respectively, and the expression stabilities of UBI, RPII, and ACA ranked as the top three genes in terms of expression stability for cold stress. In the heat stress group, HAS28, GAPDH, and ACA were identified as the three most suitable genes. In all samples, HAS28, HYP, and ADSS were the most suitable choices. However, HAS28 was the most unstable gene in the different genotypes and ADSS was the worst gene across the different developmental stages. Under heat stress and MeJA treatment, EF1 presented the lowest expression stability. Across the developmental stages, GAPDH had the lowest expression stability. GAPDH was also the worst choice for all the samples. Finally, H2A was ranked as the worst choice under cold stress.
A comparison of the different experimental sample sets showed that the most suitable genes differed. Therefore, it is necessary to select the appropriate reference gene for each experimental condition. Through RefFinder analysis, we chose the two most suitable reference genes and the least suitable gene under each experimental group to verify the accuracy of the expression of the reference genes.

Validation of the Stability of Reference Genes
The relative expression of GbCHS was calculated using the selected reference genes to validate the stability of the reference genes. Figure 5 shows that similar expression patterns of GbCHS were obtained when using the two suitable reference genes. In contrast, when an unstable reference gene was selected for expression normalization, the relative expression level of GbCHS showed an abnormal trend. The expression of GbCHS in Yezi ginkgo (YZ) and NJ significantly changed compared to that in GL when using UBI and EF1 as the reference genes. However, the expression of GbCHS did not change significantly when HAS28 was chosen as the reference gene. When HAS28 and UBI were used, the expression of GbCHS in M2-M6 was significantly changed compared with that in M1. The expression of GbCHS in M5 was 7.47 times higher than that in M1 when HAS28 was used as the reference gene. Similarly, in different tissues, when the EIF3I and RPII genes were used as references, the expression change trends were similar, and the expression of GbCHS in the female flowers was 2.57 and 2.36 times higher than that in the leaves, respectively. The expression in the female flowers was 10.72 times higher than that in the leaves when using GAPDH as a reference gene. Taking GAPDH and HAS28 as reference genes, the expression of GbCHS at 8 h was 1.73 and 1.63 times higher than that at 0 h after heat stress. However, the expression was 8.4 times higher than at 0 h when EF1 was used as the reference gene. There was no significant change in the expression of GbCHS after cold stress. Upon choosing H2A as a reference gene, the expression of GbCHS decreased significantly, and its expression significantly decreased after MeJA treatment. When using EF1 as the reference gene, we found that the expression of GbCHS at 12 h was significantly higher than that at 0 h.

Validation of the Stability of Reference Genes
The relative expression of GbCHS was calculated using the selected reference genes to validate the stability of the reference genes. Figure 5 shows that similar expression patterns of GbCHS were obtained when using the two suitable reference genes. In contrast, when an unstable reference gene was selected for expression normalization, the relative expression level of GbCHS showed an abnormal trend. The expression of GbCHS in Yezi ginkgo (YZ) and NJ significantly changed compared to that in GL when using UBI and EF1 as the reference genes. However, the expression of GbCHS did not change significantly when HAS28 was chosen as the reference gene. When HAS28 and UBI were used, the expression of GbCHS in M2-M6 was significantly changed compared with that in M1. The expression of GbCHS in M5 was 7.47 times higher than that in M1 when HAS28 was used as the reference gene. Similarly, in different tissues, when the EIF3I and RPII genes were used as references, the expression change trends were similar, and the expression of GbCHS in the female flowers was 2.57 and 2.36 times higher than that in the leaves, respectively. The expression in the female flowers was 10.72 times higher than that in the leaves when using GAPDH as a reference gene. Taking GAPDH and HAS28 as reference genes, the expression of GbCHS at 8 h was 1.73 and 1.63 times higher than that at 0 h after heat stress. However, the expression was 8.4 times higher than at 0 h when EF1 was used as the reference gene. There was no significant change in the expression of GbCHS after cold stress. Upon choosing H2A as a reference gene, the expression of GbCHS decreased significantly, and its expression significantly decreased after MeJA treatment. When using EF1 as the reference gene, we found that the expression of GbCHS at 12 h was significantly higher than that at 0 h.

Discussion
The qRT-PCR method has become a prevalent tool for analysis of gene expression. Reports have indicated that inappropriate reference genes cause contrasting conclusions in gene expression analyses [46,47]; therefore, appropriate reference genes are needed for normalization. To the best of our knowledge, this is the first study concerning the systematic selection and evaluation of reference genes based on transcriptomic data for RT-qPCR normalization in G. biloba.
Plants are challenged with numerous abiotic stresses, such as cold, heat, drought stresses, etc. These stresses have adverse effects on the production, growth, and development of plant seeds [48][49][50][51]. Plants have different cold responses, and cold-tolerant species can develop effective strategies to adapt to the cold environment [48]. The freezing tolerance of most plants from temperate regions can be improved by exposure to low, nonfreezing temperatures (0 to 10 • C) via cold acclimation [52,53]. Previous studies have shown that the expression of flavonoid synthesis genes is closely related to the response of plants to cold and heat stresses. Chalcone synthase (CHS) catalyzes the first and key regulatory step of the flavonoid biosynthesis pathway [54], and we selected one GbCHS as a target gene. For example, in Citrus sinensis, the expression of PAL, CHS, DFR, and UFGT is significantly upregulated under a low-temperature environment [55]. After heat stress, the expression of three of the four SmCHS genes (SmCHS1, SmCHS2, and SmCHS3) was continuously downregulated, while that of SmCHS4 was upregulated [56]. Similarly, the expression of PAL, 4CL, CHS, ANR, FLS, and LAR was upregulated after cold stress in Tetrastigma hemsleyanum [57]. The expression of GbCHS was also upregulated in response to cold and heat stresses in our study. MeJA has been proposed to be an important signal transduction compound involved in the process of elicitation, leading to overproduction of a variety of secondary metabolites and several environmental stress factors [58]. Several recent studies have shown that MeJA treatment can increase the contents of flavonoids and terpenes in several different species, including Hypericum perforatum [59], Norway spruce [60], Camellia sinensis [61], and G. biloba [62]. Moreover, MeJA treatment was shown to induce CHS gene expression in Coleus forskohlii [63], Carthamus tinctorius [64], and Nicotiana tabacum [65], among others. After the citrus leaves were treated with MeJA, the expression of the two CHS genes decreased but the expression increased after 24 h [66]. In Physalis angulata, CHS1 was significantly downregulated by treatment with MeJA [67]. The expression of GbCHS was also downregulated under MeJA treatment in this study.
To better normalize the expression data, selection of reference genes based on transcriptomic data has been performed in many plant species, such as Eucalyptus globulus [68], Oxytropis ochrocephala [42], Arabidopsis pumila [69], Sinocalycanthus chinensis [70], and Polygonum cuspidatum [71]. In our study, the thirteen candidate reference genes were chosen based on transcriptome sequencing data. The expression stability was evaluated in all samples. The range of amplification efficiencies was 95.76-107.28%, the range of R 2 values was 0.9801-0.999, and Cq values ranged from 20.82 to 28.70, which indicates that the qRT-PCR data can be used for subsequent analysis. We used five common methods, the geNorm, NormFinder, BestKeeper, and RefFinder programs and the ∆Ct method, to evaluate the expression stability of the thirteen candidate reference genes. RefFinder integrates the above statistical algorithms for the overall final ranking. The above results showed that the most appropriate reference genes indicated by the four algorithms were inconsistent, even within the same experimental samples, which is consistent with previous research results [7]. For example, EF1 was the most suitable reference gene identified by the BestKeeper and NormFinder algorithms in the different genotypes but ranked fourth according to the geNorm analysis. Across the different developmental stages, H2A was the most suitable choice according to the BestKeeper and geNorm analyses, while it ranked sixth and fifth according to the ∆CT and NormFinder analyses, respectively. In the MeJA treatment group, EIF3I was also the most suitable gene according to BestKeeper, but it was ranked eighth and third according to the NormFinder and geNorm analyses, respectively. For heat stress, the SD value of ACA was the lowest, so it was the most stable gene, whereas it ranked third and fourth according to the NormFinder and geNorm analyses, respectively. The most stable genes were different under the other conditions. Therefore, RefFinder should be used for a comprehensive evaluation to select the most suitable gene.
In all the samples, HAS28 and HYP were the most stable genes, and they were sufficient for normalization according to their pairwise variation (Figure 3). The stability of reference genes varied greatly with different conditions. For example, EF1 was the most suitable reference gene in different genotypes, while it was the least stable reference gene after heat stress and MeJA treatment. Across different developmental stages, UBI was the most appropriate reference gene but was not an eligible option across different tissues. GAPDH is a frequently used reference gene and is widely applied for normalization in the qRT-PCR-based analysis. For heat stress, GAPDH was the best option. In contrast, it was ranked least in the different tissues, which was consistent with the results in different issues of Sorghum bicolor [72] and Salix matsudana [38]. For MeJA treatment, ACT was the most appropriate reference gene; moreover, ACT7 was one of the best-ranked reference genes in Brassica napus [73]. ACT was also highly stably expressed in Lagerstroemia indica, and GAPDH was shown to be an eligible reference gene for Lagerstroemia speciosa at different flower developmental stages [74]; however, these reference genes were unstable in ginkgo across different leaf developmental stages. Similarly, TUB was shown to be an unsuitable gene in L. indica [74] and G. biloba. These results indicated that different experimental conditions require different reference genes, thus demonstrating the importance of selecting appropriate reference genes.

Conclusions
In our study, the expression levels of thirteen candidate reference genes were analyzed with five methods and the appropriate reference genes were indeed found to differ under the different experimental conditions. Appropriate reference genes must be selected for specific experimental conditions to obtain accurate gene expression analysis results. We selected stably expressed genes on the basis of transcriptomic data, which is more targeted, convenient, and efficient than other methods for identifying new and more stably expressed genes. The geNorm, BestKeeper, NormFinder, and RefFinder programs and the ∆Ct method were the best-suited strategies for analyzing the suitability of the reference genes, and the expression of the target gene GbCHS was normalized across all conditions. Moreover, we not only screened several eligible reference genes of G. biloba under specific conditions but also offered a theoretical basis for the selection of reference genes in other plants.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/11/1217/s1, Table S1: Information of the six Ginkgo biloba accessions, Table S2: FPKM values of 12 reference genes covering 14 transcriptomes data of G. biloba, Table S3: Ranking of the potential reference genes of G. biloba by NormFinder (The tabulated data of Figure 4a), Table S4: Ranking of the potential reference genes of G. biloba by Bestkeeper (The tabulated data of Figure 4b), Table S5: Ranking of the potential reference genes of G. biloba by Delta CT (The tabulated data of Figure 4c), Figure S1: The primer specificity and amplification products of thirteen candidate reference genes.