Internal reference genes for normalizing quantitative real-time PCR in different tissues of Gelsemium elegans or under low temperature, MeJA, and SA stresses

Quantitative real-time PCR (RT-qPCR) is a technique for analyzing gene expression levels^[1−3], which is widely used due to its advantages of sensitivity, high efficiency, and throughput. In RT-qPCR experiments, there are differences in the quality and quantity of initial RNA, cDNA synthesis efficiency, PCR amplification efficiency, and other factors among samples, resulting in the deviation of target genes in relative expression analysis^[4]. It is necessary to use one or more internal reference genes to correct and standardize the target gene expression data^[5]. Secondary metabolites play a major role in medicinal plants, and it is very important for mining functional genes in their synthetic pathways^[6]. Therefore, appropriate reference genes are the prerequisite for the analysis of key gene expression in medicinal plants.

Generally, housekeeping genes have stable expression and are often used as internal reference genes for RT-qPCR analysis. There are various internal reference genes commonly used in plants, such as the actin gene (Actin), β-tubulin gene (TUB), eukaryotic elongation factor 1-alpha (eEF-1α), polyubiquitin gene (UBQ), etc^[7,8]. However, the expression of these genes has a certain specificity. Actin was the most stable reference gene in different tissues and in leaves treated with methyl jasmonate (MeJA) in Artemisia argyi^[9]. Serinethreonine-protein phosphatase (PP2A) can be used as an optimal internal reference gene for analyzing gene expression in different tissues of Polygonum multiflorum Thunb^[10]. eEF-1α was suitable to be used as the internal reference gene under low temperature (4 °C) treatment in Chenopodium quinoa Willd^[11]. Even in the same species, the expression of housekeeping genes in different tissues, at various developmental stages, or under different biotic and abiotic stresses is variable. Not all common housekeeping genes can be applied to all species or all experimental environments^[12]. We thus reasonably deduce that selecting appropriate internal reference genes for calibration always plays an important role in the expression pattern analysis of target genes.

Gelsemium elegans (Gardner and Champ.) Benth. (G. elegans) (National Center for Biotechnology Information Taxonomy ID: 427660), an evergreen woody vine of Gelsemium in Loganiaceae, is also known as Gou Wen, Hu Man Teng, Da Cha Yao, and Duan Chang Cao, etc^[13]. It was first recorded in 'Shen Nong's Materia Medica' and used as a traditional Chinese medicinal plant with a pungent and bitter taste^[14]. G. elegans has the functions of immune regulation, anti-inflammatory, anti-anxiety, anti-tumor, and neuropathic pain^[15−19]. As a short-day medicinal plant, G. elegans is not resistant to low temperature. When G. elegans suffers from cold stress, its young leaves turn from green to red, and thus affect its normal growth and metabolism. The root, stem, and leaf of G. elegans can be used as medicines, and its main medicinal components alkaloids are distributed in different tissues^[20]. As is well known, hormones such as JA and SA affect the synthesis and accumulation of secondary metabolites in plants. At present, most of the research on the medicinal chemistry of G. elegans focuses on its main active component indole alkaloids^[21,22]. There are no previous reports on functional gene mining in the synthesis pathway and molecular biology related to the gene expression level in G. elegans. So far, the internal reference genes of G. elegans have not yet been reported under different stress conditions.

In the present study, seven commonly used housekeeping genes in plants (Actin, TUB, eEF-1α, UBQ, cullin (CUL), adenine phosphoribosyl transferase (APRT), and pseudo response regulator (PRR) which showed stable expression abundances in our previous transcriptome database of G. elegans under cold stress and in different tissues were selected and validated. GeNorm^[23,24], NormFinder^[25], BestKeeper^[26], Delta CT^[27] and RefFinder^[28] was used to analyze the expression level and expression stability of these candidate internal reference genes in different tissues (root, stem, and leaf) of G. elegans and under 4 °C stress. Moreover, the expression patterns of three candidate genes in G. elegans, including a key gene (geranyl diphosphate synthase, GPPS) in the alkaloid synthesis pathway, an ethylene-responsive transcription factor (ERF), and a structural gene 60S ribosomal protein (60S), were investigated by RT-qPCR with the selected reference genes or their combinations. This study aims to lay a foundation for more accurate quantification of gene expression in G. elegans.

G. elegans used in this study were collected from Yongding District, Longyan City, Fujian Province (altitude 602 m, 116.9° E, 24.4° N). The G. elegans plants were identified to be an evergreen woody vine of Gelsemium in Loganiaceae by Professor Zhongyi Zhang of Fujian Agriculture and Forestry University. This species is widely distributed in China, and no permits are required to collect it. The present work did not contain any studies with human participants or animals and did not involve any endangered or protected species. The experimental research on plants performed in this study complied with relevant institutional, national, and international guidelines and legislation. There were two experimental groups in this study. The first group was the three different tissue samples (root, stem, and young leaf) of G. elegans, which was referred to as DTS. The second group of G. elegans under low-temperature (4 °C) stress was referred to as LTS. Firstly, the G. elegans seedlings in Yongding District were transplanted into the greenhouse at the research institute of Fujian Agriculture and Forestry University. After 60 d of transplantation, the healthy seedlings with consistent growth were selected as experimental materials. The young leaves under low-temperature stress for 0 h (as a control), 12 h, and 24 h were collected. The leaves of G. elegans (60-day-old seedlings) were sprayed with 100 μM Methyl jasmonate (MeJA) and 5 mM salicylic acid (SA) in the greenhouse, respectively. They were the third MeJA-treated group and the fourth SA-treated group. The leaf samples were collected at 0, 6, 12, and 24 h. Each sample contained three biological replicates, and each replicate consisted of two mixed plants. All the collected samples were frozen in liquid nitrogen and stored in a −80 °C refrigerator.

Total RNA extraction and cDNA synthesis

Total RNA was extracted from different tissues and low-temperature treatment samples using TRIzol reagent (Yeasen, Shanghai, China). The concentration of RNA samples was detected by a UV spectrophotometer (BioTek, Vermont, USA). Agilent 2100 BioAnalyzer and 1% gel electrophoresis (Supplemental Fig. S1) were used to detect the purity, concentration, and integrity of the RNA samples. The first strand cDNA was reversed according to the instructions of Hifair® III 1st Strand cDNA Synthesis Kit (gDNA digester plus) (Yeasen, Shanghai, China) and detected by 1% gel electrophoresis (Supplemental Fig. S2).

Primer design for candidate internal reference genes and three G. elegans genes

Seven candidate internal reference genes (CUL, eEF-1α, Actin, TUB, UBQ, APRT, and PRR), which showed stable expression abundances in our previous transcriptome database of G. elegans under both cold and MeJA stress (unpublished) (Supplemental Table S1), were selected. In addition, the key alkaloid synthesis gene GPPS, a transcription factor ERF, and a structural gene 60S were also screened from the above transcriptome database. RT-qPCR primers were designed by Primer Premier 5.0 software (Table 1). Primer specificity was analyzed using Primer-BLAST (www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC=BlastHome). All primers were synthesized by Fuzhou Shangya Biotechnology Co., Ltd.

Primer-specific detection and amplification efficiency analysis of candidate internal reference genes

The RT-qPCR test was carried out using QuanStudio 3 real-time PCR instrument (ThermoFisher, Waltham, USA). Its reaction system (20 μL) included 10 μL SYBR Green Master Mix (Vazyme, Nanjing, China), 0.4 μL 10 μmol·L⁻¹ upstream primer, 0.4 μL 10 μmol·L⁻¹ downstream primer, 1.0 μL cDNA template, and 8.2 μL ddH₂O. RT-qPCR with distilled water as a template was performed as controls. The RT-qPCR reaction program was 50 °C for 2 min, 95 °C for 30 s, 40 cycles of 95 °C for 10 s, and 60 °C for 30 s, and the melting curve was collected at 95 °C for 15 s, 60 °C for 1 min and 95 °C for 1 s. Each sample was repeated three times. The amplification specificity of the ten candidate genes with a 5 × cDNA template of G. elegans leaves was analyzed according to the melting curve.

The cDNA template of G. elegans leaves was diluted five times successively, and five concentration gradients were set, which were 5⁻¹, 5⁻², 5⁻³, 5⁻⁴, and 5⁻⁵ times of the cDNA template, respectively. The CT values of each candidate gene under different gradients were obtained through RT-qPCR analysis, and the standard curve of the seven candidate internal reference genes (CUL, eEF-1α, Actin, TUB, UBQ, APRT, and PRR) was manufactured by Microsoft Excel 2016. Then the slope (k) and linear correlation coefficient (R²) were obtained, and the amplification efficiency was calculated following a formula of E = (10^−1/k − 1) × 100%.

Stability evaluation and validation of candidate internal reference genes

The expression stability of candidate internal reference genes was evaluated by BestKeeper, geNorm, NormFinder, Delta CT, and RefFinder software^[23−28]. For BestKeeper, the input data was the original CT value and the internal reference genes were evaluated by calculating the standard deviation (SD) of CT values of all samples and the variation coefficient (CV) of candidate genes. The input data for geNorm and NormFinder were the relative quantities converted from CT values with a formula of 2^−ΔCT, where ΔCT meant each corresponding CT value minus the minimum CT value. The data in geNorm and NormFinder were evaluated by the M value, and the lower the M value, the more stable the gene. Relative stability values for gene expression were generated in NormFinder, with smaller values indicating higher stability. The CT values of candidate internal reference genes were directly analyzed by BestKeeper and Delta CT. According to the analysis results of the four programs and the RefFinder online website (http://blooge.cn/RefFinder/?type=reference#), the stability of all seven candidates were evaluated by comprehensive comparison, and the optimal internal reference genes or their combinations of G. elegans for tissues specific expression and cold stress response was screened.

According to the optimal internal reference genes, the expression patterns of the three selected genes, including GPPS, ERF, and 60S were analyzed in different tissues of G. elegans and at different time points after low-temperature treatment. The relative expression of the gene was calculated using the 2^−ΔΔCT method^[29]. Quantitative data for internal reference genes validation were sorted to obtain the CT value of the gene and the CT value distribution map was drawn using GraphPad 6.0 software. The significance level of the data was analyzed using DPS 7.05 software, and graphs were plotted using GraphPad 6.0 and Adobe Illustrator CS6.

The total RNA of tissue samples of G. elegans (root, stem, and leaf) and G. elegans leaf samples under different stresses (cold, MeJA, and SA) for 0, 6, 12, and 24 h was extracted. The detection results showed that the concentration of RNA samples was 113.254−823.005 ng·μL⁻¹, and the ratio of absorbance at 260 and 280 nm ranged from 1.901 to 2.184 (Supplemental Table S2). These data indicate that the RNA qualities of the 18 G. elegans samples satisfy the subsequent test.

Primer design and amplification specificity of candidate internal reference genes

The melting curves for the seven candidate internal reference genes (CUL, eEF-1α, Actin, TUB, UBQ, APRT, and PRR) were all single-peak (Supplemental Fig. S3), indicating their high specificity. The cDNA templates of five concentration gradients were used to perform RT-qPCR analysis, and the obtained data was applied to draw the standard curves. The correlation coefficient (R²) values of these seven genes were all above 0.99, suggesting that the estimated primer amplification efficiency was accurate and reliable (Table 2). The slopes (k) of six genes ranged from −3.000 to −3.314 except for Actin which was −2.992. In addition, the amplification efficiency of these genes ranged from 100.32% to 115.87%, which was in the acceptable amplification efficiency between 90% and 120%^[30].

Expression abundance of candidate internal reference genes

RT-qPCR analysis results showed that the CT values of the seven candidate internal reference genes in all samples were different (Fig. 1). The average CT values ranged from 18.426 to 25.572. Among them, the UBQ gene had the highest expression abundance in the samples with a CT value lower than 20, and the expression abundance of the Actin gene was the lowest which had the highest CT value.

Figure 1. 

CT value distribution of the seven candidate internal reference genes in different tissue samples (root, stem, and leaf) and under different stresses (cold stress, MeJA stress, and SA stress) of G. elegans by GraphPad 6.0. The smaller the CT value, the higher the gene expression abundance. The change of the CT value of the same candidate internal reference gene in different samples reflects the gene expression stability. CUL, cullin; eEF-1α, eukaryotic elongation factor 1-alpha; APRT, adenine phosphoribosyl transferase; TUB, β-tubulin; Actin, actin; PRR, pseudo response regulator; UBQ, polyubiquitin.

Stability of candidate internal reference genes by geNorm, NormFinder, BestKeeper, and Delta CT

The geNorm software was used to analyze the expression stability of each candidate internal reference gene by calculating the gene expression stability value (M) (Fig. 2). It showed that the expression stability of the candidate internal reference genes in the DTS group was CUL = UBQ > eEF-1α > TUB > APRT > Actin > PRR, and was eEF-1α = APRT > CUL > Actin > TUB > PRR > UBQ in the LTS group. In the MeJA group, the order of expression stability of the candidate internal reference genes from high to low was CUL= eEF-1α > Actin > APRT > TUB > UBQ > PRR and was CUL = eEF-1α > UBQ > PRR > TUB >APRT > Actin in the SA group. The M values of the seven candidate internal reference genes in the three groups were all less than 1.5, which was regarded as relatively stable^[23,25]. The coefficient of variation (CV) of the CT value is one of the factors to measure the stability of the internal reference gene. Supplemental Table S3 showed that the CV values of CUL, APRT, and eEF-1α were relatively low, which were 3.53%, 3.61%, and 4.22%, respectively. However, the PRR gene had the highest CV value (7.05%). According to the NormFinder software instruction, NormFinder evaluates the best internal reference genes by calculating the variability of expression between and within groups^[25]. In Fig. 3, the expression stability of the candidate internal reference genes in the DTS group was CUL = UBQ > eEF-1α > TUB > Actin > APRT > PRR. CUL and UBQ showed the highest stability with the same M values (0.162), while that of the PRR was the lowest with an M value of 1.535. In the LTS group, the expression stability of the candidate internal reference genes was CUL > APRT >Actin > eEF-1α > TUB > PRR > UBQ. In the MeJA group, the whole rank of gene stability was eEF-1α > CUL > Actin > APRT > UBQ > TUB > PRR and was CUL > eEF-1α > PRR > TUB > UBQ >APRT > Actin in the SA group. The Bestkeeper program calculates the SD of the CT value of the candidate internal reference genes. As reported, reference genes with SD >1 are considered unstable and should be avoided^[26,31].

Figure 2. 

Stability of the seven candidate internal reference genes by geNorm. DTS, different G. elegans tissue samples; LTS, G. elegans leaves under low-temperature stress; CTS, both of the DTS and LTS samples. CUL, cullin; eEF-1α, eukaryotic elongation factor 1-alpha; APRT, adenine phosphoribosyl transferase; TUB, β-tubulin; Actin, actin; PRR, pseudo response regulator; UBQ, polyubiquitin.

Figure 3. 

Stability of the seven candidate internal reference genes by NormFinder. DTS, different G. elegans tissue samples; LTS, G. elegans leaves under low-temperature stress; CTS, both of the DTS and LTS samples. CUL, cullin; eEF-1α, eukaryotic elongation factor 1-alpha; APRT, adenine phosphoribosyl transferase; TUB, β-tubulin; Actin, actin; PRR, pseudo response regulator; UBQ, polyubiquitin.

As shown in Table 3, except for PRR, the SD values of the other six genes (APRT, TUB, eEF-1α, CUL, UBQ, and Actin) in the DTS group were all less than 1, indicating that these genes are more suitable as candidate internal reference genes. Regarding the expression stability, it was APRT > TUB > eEF-1α > CUL > UBQ > Actin > PRR, from high to low. In the LTS, MeJA, and SA groups, the SD values of the seven candidate genes were all less than 1, and the ranking of expression stability was CUL > eEF-1α >TUB > Actin > APRT > PRR >UBQ in the LTS group, and was CUL > eEF-1α >TUB > UBQ > APRT > Actin >PRR in the MeJA group, and was CUL > eEF-1α >PRR > UBQ > TUB > APRT > Actin in the SA group, of which CUL and eEF-1α were the most stable. The strategy of the Delta CT method is to calculate the average standard deviation of the CT value of each reference gene. The smaller the average standard deviation, the more stable the gene expression^[27]. In the DTS group, CUL, UBQ, and eEF-1α ranked among the top three candidate internal reference genes, followed by TUB, Actin, APRT, and PRR (Fig. 4). In the LTS group, the expression stability of each reference gene was CUL>Actin>eEF-1α>APRT>TUB>PRR>UBQ. In the MeJA group and SA group, eEF-1α and CUL were more stable than the other four genes (TUB, Actin, PRR, APRT, and UBQ).

Figure 4. 

Stability of the seven candidate internal reference genes by Delta CT. DTS, different G. elegans tissue samples; LTS, G. elegans leaves under low-temperature stress; CTS, both of the DTS and LTS samples. CUL, cullin; eEF-1α, eukaryotic elongation factor 1-alpha; APRT, adenine phosphoribosyl transferase; TUB, β-tubulin; Actin, actin; PRR, pseudo response regulator; UBQ, polyubiquitin.

In summary, the geNorm, NormFinder, Bestkeeper, and Delta CT analyses all revealed that CUL was the highest stable gene in different groups, but there were certain differences in other genes, which may be due to different software and the gene stability. To determine the comprehensive stability of genes, RefFinder analysis was further conducted for comprehensive stability in this study.

Comprehensive stability of candidate internal reference genes

RefFinder method avoids the one-sidedness of single algorithm, and the smaller the stable value, the better the stability of the internal gene^[28]. Figure 5 showed that the comprehensive stability order (from stable to unstable) of the candidate internal reference genes in the DTS group and LTS group was CUL > UBQ > eEF-1α > TUB > APRT > Actin > PRR and CUL> eEF-1α > APRT > Actin> TUB > PRR > UBQ, respectively. The ranking order of the comprehensive stability of the seven reference genes in the MeJA group and SA group was eEF-1α > CUL >Actin > APRT > TUB> UBQ >PRR and CUL > eEF-1α > PRR > UBQ > TUB > APRT >Actin, respectively. Therefore, CUL has the highest expression stability in different G. elegans tissues, as well as under low-temperature treatment and SA treatment, which is suitable as an RT-qPCR internal reference gene in G. elegans, while under MeJA treatment it was eEF-1α. Followed by the comprehensive stability of CUL, it was eEF-1α, which could also be used as an internal reference gene for G. elegans in the LTS group and SA group. CUL could also be used as a good candidate gene for internal reference in the MeJA group. In addition, PRR has the least stability among the seven candidate internal reference genes in different G. elegans tissues and under MeJA stress, and the same as UBQ under low-temperature stress, Actin under SA stress, suggesting that they may not be suitable as internal reference genes for G. elegans under these conditions.

Figure 5. 

Comprehensive stability of the seven candidate internal reference genes. DTS, different G. elegans tissue samples; LTS, G. elegans leaves under low-temperature stress; CTS, both of the DTS and LTS samples. CUL, cullin; eEF-1α, eukaryotic elongation factor 1-alpha; APRT, adenine phosphoribosyl transferase; TUB, β-tubulin; Actin, actin; PRR, pseudo response regulator; UBQ, polyubiquitin.

Stability verification for internal reference genes by the expression pattern of GPPS, ERF, and 60S genes

To verify the applicability of candidate internal reference genes, two reference genes (CUL and eEF-1α) with the highest comprehensive scores were obtained from the RefFinder software. They were composed of single (CUL and eEF-1α), and double (CUL + eEF-1α) internal reference genes to verify the relative expression of GPPS, ERF, and 60S genes in DTS, LTS, MeJA, and SA samples, respectively. RT-qPCR analysis showed (Fig. 6) that the expression trends of GPPS, ERF, and 60S in the DTS, LTS, MeJA, and SA groups with the single internal reference gene of CUL or eEF-1α were consistent with that using double (CUL + eEF-1α) reference genes. In the DTS group, the expression patterns of GPPS and ERF were similar. The highest expression level of GPPS and ERF was in the root, followed by that in the leaf, and the lowest in the stem. However, the 60S gene showed the highest expression level in the root, but the lowest in the leaf. In the LTS group, the expression level of GPPS was unchanged under low-temperature stress from 0 h to 12 h, but significantly increased at 24 h. The expression level of ERF was decreased at 12 h but increased at 24 h. For 60S, its expression level was decreased at 12 h and remained stable at 24 h. In the MeJA group, the expression patterns of GPPS, ERF, and 60S were similar. The expression level of GPPS, ERF, and 60S showed a tendency to increase and then decrease, with both GPPS and ERF with the highest expression at 6 h, followed by 12 h, and 60S with the highest expression at 6 h and 12 h. In the SA group, the expression level of GPPS was the highest at 12 h, followed by 24 h. The expression level of ERF showed a continuously increasing trend, with the highest at 24 h, and the next at 12 h. The expression level of the 60S was a significant difference from 0 to 24 h, among them, 12 h > 6 h > 24 h > 0 h. These results indicated that the single internal reference gene of CUL or eEF-1α, and the combinations of double (CUL + eEF-1α) internal reference genes were stable reference genes for RT-qPCR analysis in different G. elegans tissues and under low temperature, MeJA, and SA stresses.

Figure 6. 

Expression patterns of GPPS, ERF, and 60S in different tissues (DTS) of G. elegans and at different time points after low temperature (LTS), MeJA, and SA stresses. CUL, cullin; eEF-1α, eukaryotic elongation factor 1-alpha; GPPS, geranyl diphosphate synthase; ERF, ethylene-responsive transcription factor; 60S, 60S ribosomal protein. All data points were means ± standard error (n = 3). Data analysis was conducted using DPS 7.05 software, and Duncan's new multiplex range test was employed to assess the significance of differences. Lowercase letters were utilized to indicate significant distinctions between groups, with each letter corresponding to a unique group. (p-value < 0.05).

RT-qPCR is an indispensable experimental technique in the research on the synthesis and regulation of secondary metabolites in medicinal plants, which is an important method to detect gene expression^[32]. Selecting an appropriate internal reference gene for calibration can improve the accuracy of RT-qPCR, and the appropriate internal reference gene should be selected according to specific experimental conditions^[33,34]. geNorm, NormFinder, BestKeeper, and Delta CT are commonly used methods for assessing the stability of reference genes. geNorm can identify the most stable combination of genes, but it incurs higher computational costs for a large number of samples^[23]. NormFinder considers both within-sample and between-sample variations, suitable for small samples but with higher computational cost^[35]. BestKeeper provides a comprehensive stability index, applicable to small samples but sensitive to outliers^[36]. The Delta CT method is simple and practical but overlooks potential mutual influences among genes^[37]. The choice of method should consider experimental conditions, and integrating results from multiple methods enhances the reliability of the assessment.

In this study, seven candidate internal reference genes were screened from the G. elegans transcriptome database. Through geNorm, NormFinder, BestKeeper, and Delta CT analysis (Figs 2−4; Table 3), the expression stability of these seven candidates in the root, stem, and leaf of G. elegans and the leaf at different time points after low temperature, MeJA, and SA stresses was comprehensively evaluated. The ranking of candidate internal reference genes varied slightly among algorithms^[38]. In general, the results of this study showed that the stability of CUL and eEF-1α was essentially better with the above four analyzed platforms in the four experiments, which were in accordance with the previous reports that both CUL and eEF-1α can be used as internal reference genes for gene expression normalization in different plant tissues and under low temperature, MeJA, and SA stresses^[39−41]. To synthesize the results of the four algorithms, RefFinder was utilized to rank the identified candidate genes in the G. elegans, and the platform played an important role in integrating the results of the other algorithms for the screening of internal reference genes^[42,43]. Fortunately, it was found that RefFinder's results were similar to those of the different algorithms in the four experiments (Fig. 5; Table 3).

GPPS is a key enzyme in the synthesis pathway of iridoid compounds^[44] and in the generation of indole alkaloid precursors. ERF is a key hub for regulating hormone and stress signals^[45] and a downstream regulator in the ethylene signal transduction pathway^[46], which plays an important role in plant stress response^[44]. Ribosomal proteins are important structural components of ribosomes^[47]. In this study, the single (CUL or eEF-1α), and double (CUL + eEF-1α) internal reference genes were used as internal reference genes to detect the relative expression level of GPPS, ERF, and 60S genes in leaf, root, and stem of G. elegans and under low temperature, MeJA, and SA stresses. Our data indicated that GPPS, ERF, and 60S genes were constitutively expressed in G. elegans tissues and significantly enriched in the root (unpublished). The content of the important medicinal ingredient koumine (indole alkaloid) of G. elegans in the root is significantly higher than that in other tissues^[20]. It is thus speculated that GPPS may play an important role in the synthesis of indole alkaloids in G. elegans. Prolonged periods of sustained cold temperatures have serious effects on G. elegans growth^[48]. The expression level of ERF and 60S genes could be induced by low-temperature (Fig. 6). A previous study found that ERF plays an active role in the cold resistance of Citrus reticulata Blanco^[49], and it is inferred that ERF can be used as a candidate gene for G. elegans in response to low-temperature stress. MeJA is a volatile methyl ester of jasmonic acid, which has been identified as the main signaling molecule in abiotic and biological stress^[50]. SA is a well-known inducer of systemic acquired resistance and could also promote the synthesis of secondary metabolites in plants. MeJA and SA as exciters were able to promote the synthesis of different active compounds in a variety of medicinal plants in different cultures^[51]. GPPS, ERF, and 60S genes were induced to be expressed under MeJA and SA treatments, respectively. Active participation of ERF in hormone signaling pathways enhances the synthesis of medicinally active plant compounds in various medicinal plants^[52]. In addition, the expression trends of GPPS, ERF, and 60S in the DTS, LTS, MeJA, and SA groups with the internal reference gene of CUL or eEF-1α were consistent with that using double reference genes (Fig. 6). For economic reasons, the best quantitative results can be obtained by using fewer reference genes. Therefore, single reference genes (CUL or eEF-1α) or double reference genes (CUL + eEF-1α) can be selected as the best or perfect combination of reference genes in subsequent studies of G. elegans.

In the present study, CUL and eEF-1α were selected to be the optimal internal reference gene for RT-qPCR analysis in different G. elegans tissues and under low temperature, MeJA, and SA stresses according to the combination analysis of geNorm, NormFinder, BestKeeper, Delta CT, and RefFinder programs. Moreover, single reference gene (CUL or eEF-1α) or double reference genes (CUL + eEF-1α) can be selected as the best or perfect combination of reference genes for gene expression analysis in G. elegans. This study not only lays a foundation for the expression analysis of related genes in the growth and development of G. elegans under different conditions, but also provides a reference for the subsequent studies on the expression and functional verification of other genes in G. elegans.

The authors confirm contribution to the paper as follows: study conception and design: You C, Zang S, Que Y, Que W; data collection: You C, Cui T, Sun X; analysis and interpretation of results: Su Y, Lin Q, Lin H; draft manuscript preparation: You C, Que Y, Que W. All authors reviewed the results and approved the final version of the manuscript.

The data supporting the conclusions of this article are within the paper.