DOI: https://doi.org/10.21203/rs.3.rs-35251/v1
Background: Before studying the gene expression of different organisms, it is important to determine the best reference gene. At present, the most accurate method for detecting gene expression is quantitative real-time PCR(RT-qPCR). By using this method, reference genes that are stable in different biological systems and under different conditions can be obtained. Toona ciliata Roem (T. ciliata). is a valuable and fast-growing timber species. In this study, 20 reference genes were identified through RT-qPCR, as a primary prerequisite for future gene expression analysis. Four different methods, geNorm, NormFinder, BestKeeper, and RankAggreg were used to evaluate the stability of expression of 20 candidate reference genes in various tissues under different conditions.
Results: The experimental results showed that TUB-α was the most stably expressed reference gene across all samples; UBC17 was found to be the most stable in leaves & young stems under Hypsipyla robusta (H. robusta) and methyl jasmonate (MeJA) treatment. In addition, under H. robusta treatment, PP2C59 and UBC5B were the best-performing genes in leaves, while HIS1 and ACT7 were the best reference genes in young stems. Under low temperature (4℃) treatment, the two best reference genes were 60S-18 and TUB-α. The expression of HIS6 and MUB1 was the most stable under PEG6000 treatment. The accuracy of the selected reference genes was verified using transcription factor MYB3(TcMYB3) genes.
Conclusions: This is the first report to verify the best reference genes for normalizing gene expression in T. ciliata under different conditions, which will facilitate the future elucidation of gene regulations in this species.
Toona ciliata Roem. belongs to the Meliaceae family, which is widely distributed in China, Australia and India. Because of its straight trunks and russet wood, T. ciliata has the title of "Chinese Mahogany" [1]. But due to environmental degradation and destruction by humans, its population has decreased sharply, therefore, has been listed as one of the "National Class II Key Protected Endangered Plants" in China. T. ciliata has great economic value that its wood is often used to produce high-end furniture, high-end musical instruments and high-quality crafts [2], and more importantly, it is often used as a medicinal plant [3]. Its roots, stems, and leaves have rich chemical substances. Compounds that have been isolated from T. ciliata include ketones, steroids and coumarins, many of which have antifungal, anti-glycation or anti-tumor activity [4–7], and the flower extract has a therapeutic effect on gastric ulcers [8]. However, the yield of compounds isolated from T. ciliata is low. In addition, in previous research, it has been found that T. ciliata is very susceptible to the moth pest Hypsipyla robusta Moore [9] which eats mainly the young stems and causes the hollow branches to fail to grow and die in some cases. This pest is not only a regional issue in China, but also a worldwide problem. In some of the main areas where H. robusta is distributed, such as Australia and Brazil, T. ciliata also faces serious damage from H. robusta [10–12]. At present, there are no chemical and physical methods for effectively preventing and controlling H. robusta, and those that are available usually consume a lot of labor work and material resources, therefore, are not applicable to large-scale forest plantations [13]. Obtaining insect-resistant plants through molecular breeding may help to solve the problem of damage from the pest. In order to synthesize a desired compound related to the resistance mechanism, it is necessary to first explore the pathway and its related regulatory genes in those biological pathways [14, 15]. In exploring biosynthetic and insect-resistance mechanisms in T. ciliata, gene expression analysis is one of the most important steps. So far, the knowledge base ICG (http://icg.big.ac.cn) has collected reference genes from more than 120 plant species including Arabidopsis [16], peanut [17], cucumber [18], soybean [19], but it does not include T. ciliata. Nor are there any literature references about T. ciliata housekeeping genes which can be used for the standardization of gene expression.
RT-qPCR has advantages, including good repeatability, high sensitivity, accurate quantification, and fast reaction, making it a powerful tool for carrying out the entire PCR process and the most commonly used method for determining gene expression levels [20]. However, RT-qPCR can be affected by multiple sources of error, such as the amount of starting material, the integrity of the RNA, and the efficiency of the enzymatic reaction. It is therefore necessary to introduce a relatively stable expressed housekeeping gene as a reference for correction and standardization, so as to control the unnecessary errors generated within and between samples [21].
The commonly used housekeeping genes are those encoding actin (ACT), glyceraldehyde-3 phosphate dehydrogenase (GAPDH), and tubulin (TUB), because they are considered to be consistently expressed under all conditions [22]. However, more and more studies are now reporting that genes that are stably expressed across different tissues, different experimental conditions, and different species do not exist. In order to ensure the accuracy of an experiment, it is important to evaluate the reference genes to select those suitable reference genes for specific experimental conditions [23]. Software packages are widely used to assess the stability of expression of candidate reference genes, and determine the one or more best choices include geNorm [24], NormFinder [25], BestKeeper [26]. Many researchers have used these algorithms to successfully identify reference genes in different species [27, 28]. In species for which a reference genome sequence is available, which have generally been studied in depth, the use of reference genes in expression analysis has greatly facilitated research into plant development and evolutionary mechanisms [29].
In this study, 20 candidate genes from T. ciliata transcriptome database generated by our group were screened: actin 7 (ACT7), phosphoglycerate kinase (PGK), 60S ribosomal protein L13 (60S-13) and L18 (60S-18), histone deacetylase 1 (HIS1) and 6 (HIS6), protein phosphatase 2 C57 (PP2C57) and C59 (PP2C59), ubiquitin-conjugating enzyme E2 5B (UBC5B) and 17 (UBC17), S-adenosylmethionine decarboxylase proenzyme (SAMDC), elongation factor 1 (EF1) and 2 (EF2), peptidyl-prolyl cis-trans isomerase CYP95 (PPIA95) and CYP26-2 (PPIA26), 18S rRNA (18S), tubulin alpha-3 chain (TUB-α), tubulin beta-5 chain (TUB-β), membrane-anchored ubiquitin-fold protein 1 (MUB1), and TIP41-like protein (TIP41). The study was carried out to determine the most suitable T. ciliata candidate gene(s) as the reference(s) for gene expression analysis using RT-qPCR technique under specific conditions including different tissues (mature leaves, young leaves, flowers, shoots and young stems) and treatments (4℃, MeJA, PEG6000 and H. robusta). In addition, TcMYB3 gene was used to confirm the reliability and validity of the reference genes screened. MYB proteins, which constitute one of the largest family of transcription factors in plants, play important roles in plant growth and development, biotic and abiotic stress responses, and circadian rhythm regulation [30, 31]. For example, the R2R3 MYB transcription factor MdMYB30 modulates plant resistance against pathogens, and Arabidopsis transcription factor MYB102 increases plant susceptibility to aphids [32, 33]. Our research provided the best reference genes for use in RT-qPCR analysis of T. ciliata under different conditions, laying a foundation for studying molecular mechanisms in T. ciliata through gene expression analysis.
Reverse-transcribed cDNA from each sample was used as a template with primers for standard PCR amplification. After detection by 1.5% agarose nucleic acid electrophoresis, it was proved that all PCR products were accurate with single bands in the gel (Fig. S1). The melting profiles of all primers amplified candidate reference genes using RT-qPCR showed single peaks (Fig. S2). A standard curve for each candidate was obtained by stepwise dilution, and their linear correlation coefficients were all greater than 0.99 (R2 > 0.99). The amplification efficiency for the 20 candidate reference genes was from 90.41% for PPIA95 to 102.44% for PGK. Further details of primers are given in Table 2.
Experimental design |
tissue |
Biological repetition |
Sampling points |
Number of samples |
---|---|---|---|---|
Different tissues |
mature leaves, young leaves, flowers, shoots, young stems, roots |
3 |
1 |
18 |
H. robusta treatment |
leaves, young stems |
3 |
6 |
36 |
4 °C treatment |
leaves |
3 |
5 |
15 |
MeJA treatment |
leaves |
3 |
5 |
15 |
PEG6000 treatment |
leaves |
3 |
5 |
15 |
The expression levels of all candidate reference genes were determined by RT-qPCR under all of the following conditions: different tissues, H. robusta treatment, 4℃ treatment, MeJA treatment, and PEG6000 treatment. The expression levels of candidate genes were very different across the samples. The maximum cycle threshold (CT) value was 31.66, and the minimum was 13.18 (Fig. 1). Among them, PPIA26 showed the highest expression abundance, the maximum, minimum and median of the CT values being 31.66, 20.07 and 23.36, respectively. EF1 showed the lowest expression abundance, with maximum, minimum and median CT values of 22.16, 13.18, and 17.17 respectively. In addition, candidate genes exhibited significant variability in expression. MUB1 and UBC5B had a relatively narrow range of CT values compared with other genes, indicating that they are more stably expressed. However, more importantly, these results show that none of the genes is expressed stably across all conditions, so it is necessary to screen reference genes for T. ciliata under specific conditions.
The software packages geNorm, NormFinder and BestKeeper were used to evaluate the stability of expression of the 20 candidate reference genes under different experimental conditions. The R software RankAggreg package was used for overall ranking [34].
In geNorm analysis, M value is calculated for each pair of genes. The stability of gene expression is evaluated based on the M value; the genes with threshold M value below 1.5 are considered as stably expressed, and the gene with the lowest M value is regarded as the most stably expressed reference gene. The results of geNorm analysis of 20 candidate reference genes under different conditions are shown in Fig. 2A-H. In our processed sample groups, the M values of all candidate genes were below 1.5 (Fig. 3). Under H. robusta treatment, UBC17, PP2C59 and UBC5B were most stably expressed in leaves (Fig. 2A), while HIS1, UBC5B and ACT7 were with few expression fluctuations in young stems (Fig. 2B). Analyses of data from the two tissues under H. robusta treatment showed that the M values of UBC5B, HIS1 and ACT7 were the smallest, and therefore, are considered with the most stable expression (Fig. 2C). The most stable genes expressed across different tissues were 18S and TUB-α, with M values around 0.2 (Fig. 2D). PPIA95 showed good stability under both 4℃ and MeJA treatments; another gene with stable expression under 4℃ treatment was 60S-18, and UBC17 is another gene with the most stable expression under MeJA treatment (Fig. 2E and F). Under drought stress which was simulated by PEG6000 treatment, the two genes with the lowest M values were PP2C57 and EF1 (Fig. 2G). When data from all the above experimental conditions were analyzed, it was found that EF2 and EF1 had the highest stability with M value 0.49 (Fig. 2H).
In general, multiple reference genes are more reliable than a single reference gene when used for quantitative gene analysis. For this reason, geNorm calculates the pairwise variation (Vn/n+1) of the normalization factor after the introduction of a new reference gene, and determines the optimal number of reference genes based on this ratio. The default Vn/n+1 value for the software is 0.15. If the ratio is less than 0.15, it indicates that the number of internal gene combinations that can meet the requirements for relative quantification is n, otherwise the next reference gene needs to be introduced. In our study, the values of pairwise variation V2/3 under conditions with H. robusta treatment, 4℃ treatment, MeJA treatment, PEG6000 treatment, and different tissues, were all less than 0.15, which means that the optimal number of reference gene combinations is two (Fig. 3). Across all samples, pairwise variation (V2/3) was 0.180, V3/4 was 0.15, and V4/5 was 0.126, indicating that the addition of the third and fourth reference genes has different impacts on the results. However, less reference genes are always less time and cost consuming, which is desired by researchers. For all samples, the best reference gene combination was found to be EF2, EF1, ACT7 and UBC5B.
In order to further determine the stability of candidate reference genes, NormFinder was used to re-analyze the data. The results are shown in Fig. 4A-H. Under H. robusta treatment, the top three genes for stability of expression in leaves were TUB-α (stability value = 0.038), HIS1 (0.105), and PP2C59 (0.161) (Fig. 4A), while the most stable in young stems were ACT7 (0.042), UBC5B (0.042) and TIP41 (0.109) (Fig. 4B). The top three ranked reference candidate two tissues (leaves and stems) were TUB-α (0.170), UBC5B (0.206) and PPIA95 (0.250) (Fig. 4C). The genes 18S (0.088) and TUB-α (0.112) had lower stability values across different tissues since they showed the most stable expression, which is consistent with the results of geNorm analysis (Fig. 4D). However, ACT7 (0.014), TUB-α (0.082) and PGK (0.099) were found to be the most stable candidate genes under 4℃ treatment (Fig. 4E), which is inconsistent with the results of geNorm analysis; this may be because the two software packages use different algorithms. Under MeJA treatment, the two most stable reference genes were 18S (0.078) and UBC17 (0.100) (Fig. 4F), while under PEG6000 treatment, the two most stable were 18S (0.055) and HIS6 (0.082) (Fig. 4G). For all samples, TUB-α (0.281), 18S (0.316), and PGK (0.335) were the three genes with the lowest stability values (Fig. 4H) that was different from geNorm analysis, because PPIA26, PP2C59 and HIS6 were the most unstable genes in geNorm.
BestKeeper takes as input the CT data for each gene-primer pair combination and calculates the coefficient of variation (CV) and standard deviation (SD), as shown in Table 3. The stability of genes is judged by the value of CV ± SD. More stable genes have a lower value of CV ± SD. Under H. robusta treatment, UBC17 (1.16 ± 0.23) and 18S (1.62 ± 0.36) were the most stable genes in leaves, the expression of HIS1 (1.40 ± 0.31) and UBC17 (1.70 ± 0.33) were the most stable in young stems, and genes with the most stable expression across two tissues was UBC17 (1.73 ± 0.34) and 18S (2.12 ± 0.47), as was the case in leaves. HIS6 (5.41 ± 1.38) and MUB1 (5.48 ± 1.16) were the most stably expressed genes in the different tissues, 18S (0.73 ± 0.16) and 60S-18 (0.83 ± 0.18) in 4℃ treatment, 60S-18 (0.67 ± 0.14) and EF1 (0.83 ± 0.14) in MeJA treatment, PPIA26 (3.89 ± 0.87) and 60S-18 (4.48 ± 1.02) in PEG6000 treatment. For all samples, MUB1 (3.95 ± 0.87) showed the highest value for expression stability.
Gene symbol |
Gene Name |
Primer: Forward/reverse |
Amplification product size (bp) |
standard curve |
En |
R2 |
|
---|---|---|---|---|---|---|---|
ACT7 |
Actin-7 |
F: TGATTGGGATGGAAGCAGCA R: GAACATGGTTGAACCGCCAC |
122 |
y = -3.5133x + 29.711 |
0.9259 |
0.9931 |
|
PGK |
Phosphoglycerate kinase |
F: CCGCAAGCTTCTTTGCGATT R: GGCTTGGATATTGGACCCGA |
145 |
y = -3.2649x + 26.93 |
1.0244 |
0.9985 |
|
60S-18 |
60S ribosomal protein L18a-1 |
F: GCCTGGATGCCTTGTATGTTG R: GGGAAAGCACCAAGCAGTTTC |
108 |
y=-3.5672x + 27.862 |
0.9332 |
0.9993 |
|
60S-13 |
60S ribosomal protein L13-1 |
F: CCAACATGGCACTCATTCGC R: TTCCCAAGATGTGCTCGCAA |
200 |
y = -3.4076x + 29.173 |
0.9654 |
0.9929 |
|
HIS6 |
Histone deacetylase 6 |
F: ATTGTCCGGTGATAGGTTGGG R: GTCTCGTAGCACCAACAACG |
153 |
y=-3.4932x + 29.484 |
0.9332 |
0.9965 |
|
HIS1 |
Histone acetyltransferase MCC1 |
F: CTGCACGAATTGTGCTGGTC R: ACTGCACGACATGTTGGGAT |
193 |
y = -3.5229x + 30.411 |
0.9225 |
0.9959 |
|
PP2C57 |
Protein phosphatase 2C 57 |
F: TGTTGCAGCTTTACAAGGCG R: TGAACAAATCACCGCCTCCA |
185 |
y = -3.3057x + 32.193 |
1.0068 |
0.9949 |
|
PP2C59 |
Probable protein phosphatase 2C 59 |
F: TAAGCGATCGCCAACAAGGA R: CACGAGCTGCTGAGTATGTGA |
194 |
y = -3.3119x + 27.157 |
1.0042 |
0.9974 |
|
UBC5B |
Ubiquitin-conjugating enzyme E2 5B |
F: GGAGGACCCATGATTGTTGC R: TCGAAGCGGATCTTGAAGGAG |
116 |
y = -3.3156x + 25.529 |
1.0026 |
0.9987 |
|
UBC17 |
Ubiquitin-conjugating enzyme E2-17 kDa |
F: GCGTCGAAACGCATCTTGAA R: GAAACACCCCTCCCGCATAA |
148 |
y = -3.4489x + 26.966 |
0.9496 |
0.999 |
|
EF1 |
elongation factor 1-alpha-like |
F: CCGACCTTCTTCAGGTAGGAA R: TCCAAGGATGGTCAGACTCG |
164 |
y = -3.4295x + 23.45 |
0.9570 |
0.9915 |
|
EF2 |
elongation factor 1-alpha-like |
F: CACCCTTGGTGTGAAGCAAA R: GGTTGGTGGACCTCTCAATCA |
200 |
y = -3.4128x + 27.077 |
0.9561 |
0.9992 |
|
PPIA26 |
Peptidyl-prolyl cis-trans isomerase CYP26-2 |
F: GAAGCTGAAGTTGGTTGCCC R: GACGACCAGGGCTGAAACAT |
147 |
y = -3.4393x + 30.33 |
0.9532 |
0.9952 |
|
PPIA95 |
Peptidyl-prolyl cis-trans isomerase CYP95 |
F: ACCCGGCCTCTTATCTATGC R: ACAAGCTCCCCGAATACCAC |
117 |
y = -3.4295x + 23.45 |
0.9041 |
0.9915 |
|
SAMDC |
S-adenosylmethionine decarboxylase proenzyme |
F: AGCGATCTGCTATGACCCTG R: CCCGCAGAACCTGATTGGTC |
102 |
y = -3.3179x + 29.683 |
1.0017 |
0.9994 |
|
18S |
18S rRNA factor 2 |
F: GCTGCTAAGAGAGAGCGGG R: GGGAGCTCAGAATGGGTTCG |
128 |
y = -3.4317x + 30.095 |
0.9561 |
0.9987 |
|
TUB-α |
Tubulin alpha-3 chain |
F: TACAACAGTTGGCGGCTGAT R: TGTACCGCGGAGATGTTGTT |
137 |
y = -3.3568x + 29.206 |
0.9857 |
0.9998 |
|
TUB-β |
Tubulin beta-5 chain |
F: ACACACGCTGGACTTGACAT R: TCGCTACCTAACTGCTTCGG |
139 |
y = -3.3349x + 32.62 |
0.9946 |
0.9996 |
|
MUB1 |
Membrane-anchored ubiquitin-fold protein 1 |
F: GCATTCTTGCTCAATGGCCT R: GGTTGTAACTCCACCAGGGA |
152 |
y = -3.3956x + 28.572 |
0.9701 |
0.9984 |
|
TIP41 |
TIP41-like protein |
F: TGGTTGGAAGCAGGAAGGTT R: TTCACTTCCGCAGTATGGTG |
133 |
y = -3.3332x + 31.127 |
0.9953 |
0.9992 |
|
MYB3 |
Transcription factor MYB3 |
F: CGCACCCATAACAACTCCCA R: TCTTTCACTTACTCCCTCTTCAGC |
178 |
y = -3.4246x + 32.543 |
0.9589 |
0.9968 |
In this study, three algorithms were used to analyze the stability of expression stability of 20 candidate reference genes. However, because they use different algorithms, the gene ranking tables generated by them are also different. RankAggreg is an algorithm designed to aggregate large ranking lists. It performs aggregation of ordered lists based on the rankings via the Cross-Entropy Monte Carlo algorithm or a Genetic Algorithm [34]. To provide a consensus ranking, we used RankAggreg to calculate the overall gene ranking for each experimental condition, as shown in Table 4. The consensus for the top two genes in H. robusta treatment on leaves and under MeJA treatment was consistent with the results of geNorm analysis. HIS1 ranked first for young stem tissue under H. robusta treatment. The first-placed genes under 4℃ treatment and PEG6000 treatment were 60S-18 and HIS6. TUB-α was the most stable gene in different tissues and all samples. The expression of PPIA26 was the most unstable under all experimental conditions except PEG6000 treatment.
Ranking |
H. robusta -leaves |
H. robusta -young stems |
H. robusta -leaves & young stems |
Different tissues |
4 °C treatment |
MeJA treatment |
PEG6000 treatment |
All samples |
||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
gene |
CV ± SD |
gene |
CV ± SD |
gene |
CV ± SD |
gene |
CV ± SD |
gene |
CV ± SD |
gene |
CV ± SD |
gene |
CV ± SD |
gene |
CV ± SD |
|
1 |
UBC17 |
1.16 ± 0.23 |
HIS1 |
1.40 ± 0.31 |
UBC17 |
1.73 ± 0.34 |
HIS6 |
5.41 ± 1.38 |
18S |
0.73 ± 0.16 |
60S-18 |
0.67 ± 0.14 |
PPIA26 |
3.89 ± 0.87 |
MUB1 |
3.95 ± 0.87 |
2 |
18S |
1.62 ± 0.36 |
UBC17 |
1.70 ± 0.33 |
18S |
2.12 ± 0.47 |
MUB1 |
5.48 ± 1.16 |
60S-18 |
0.83 ± 0.18 |
EF1 |
0.83 ± 0.14 |
60S-18 |
4.48 ± 1.02 |
PPIA95 |
4.16 ± 0.96 |
3 |
HIS6 |
2.05 ± 0.48 |
TUB-β |
1.84 ± 0.44 |
MUB1 |
2.54 ± 0.56 |
PPIA95 |
6.24 ± 1.36 |
PPIA95 |
0.90 ± 0.21 |
HIS1 |
0.85 ± 0.19 |
PPIA95 |
5.10 ± 1.24 |
18S |
4.58 ± 1.02 |
4 |
UBC5B |
2.06 ± 0.37 |
18S |
1.87 ± 0.41 |
HIS6 |
2.69 ± 0.62 |
SAMDC |
6.42 ± 1.34 |
EF1 |
1.13 ± 0.21 |
TUB-β |
0.87 ± 0.20 |
ACT7 |
5.18 ± 1.18 |
PGK |
4.90 ± 0.99 |
5 |
PP2C59 |
2.31 ± 0.45 |
PGK |
2.05 ± 0.40 |
HIS1 |
3.09 ± 0.70 |
60S-18 |
6.57 ± 1.34 |
PGK |
1.19 ± 0.24 |
PPIA95 |
0.88 ± 0.20 |
HIS6 |
5.26 ± 1.42 |
PP2C59 |
5.13 ± 1.02 |
6 |
SAMDC |
2.84 ± 0.63 |
ACT7 |
2.25 ± 0.44 |
UBC5B |
3.09 ± 0.54 |
UBC5B |
7.32 ± 1.30 |
TUB-α |
1.32 ± 0.27 |
TIP41 |
0.89 ± 0.20 |
TUB-β |
5.32 ± 1.36 |
TUB-α |
5.14 ± 1.09 |
7 |
MUB1 |
2.95 ± 0.65 |
MUB1 |
2.26 ± 0.50 |
PPIA95 |
3.17 ± 0.74 |
EF2 |
7.47 ± 1.48 |
ACT7 |
1.51 ± 0.33 |
UBC17 |
0.97 ± 0.19 |
PGK |
5.43 ± 1.17 |
TIP41 |
5.28 ± 1.17 |
8 |
HIS1 |
3.06 ± 0.71 |
TIP41 |
2.46 ± 0.52 |
PP2C59 |
3.31 ± 0.66 |
HIS1 |
7.54 ± 1.77 |
SAMDC |
1.58 ± 0.34 |
60S-13 |
0.97 ± 0.21 |
SAMDC |
5.63 ± 1.32 |
SAMDC |
5.33 ± 1.15 |
9 |
ACT7 |
3.12 ± 0.65 |
PPIA95 |
2.75 ± 0.64 |
TUB-α |
3.42 ± 0.72 |
PGK |
27.62 ± 1.51 |
EF2 |
1.65 ± 0.35 |
SAMDC |
1.34 ± 0.28 |
PP2C57 |
5.76 ± 1.24 |
UBC5B |
5.36 ± 0.97 |
10 |
TUB-α |
3.17 ± 0.68 |
UBC5B |
2.81 ± 0.48 |
ACT7 |
3.66 ± 0.74 |
18S |
7.75 ± 1.66 |
TUB-β |
1.67 ± 0.43 |
HIS6 |
1.49 ± 0.34 |
18S |
6.07 ± 1.44 |
ACT7 |
5.37 ± 1.14 |
11 |
TIP41 |
3.46 ± 0.78 |
HIS6 |
2.87 ± 0.65 |
SAMDC |
3.85 ± 0.84 |
PP2C59 |
7.95 ± 1.59 |
HIS6 |
1.92 ± 0.44 |
18S |
1.57 ± 0.34 |
MUB1 |
6.34 ± 1.43 |
UBC17 |
5.48 ± 1.11 |
12 |
PPIA95 |
3.61 ± 0.84 |
60S-13 |
3.05 ± 0.62 |
PGK |
3.96 ± 0.80 |
ACT7 |
8.07 ± 1.68 |
60S-13 |
2.09 ± 0.47 |
EF2 |
1.64 ± 0.32 |
UBC5B |
6.37 ± 1.27 |
60S-18 |
5.56 ± 1.19 |
13 |
PGK |
4.16 ± 0.86 |
TUB-α |
3.36 ± 0.70 |
TUB-β |
3.98 ± 0.93 |
TIP41 |
8.08 ± 1.77 |
MUB1 |
2.11 ± 0.47 |
ACT7 |
1.68 ± 0.35 |
EF2 |
6.42 ± 1.43 |
HIS1 |
5.86 ± 1.37 |
14 |
PP2C57 |
4.35 ± 0.88 |
EF2 |
3.52 ± 0.65 |
TIP41 |
4.34 ± 0.94 |
TUB-α |
8.09 ± 1.70 |
PP2C57 |
2.21 ± 0.44 |
PP2C57 |
1.84 ± 0.38 |
EF1 |
6.61 ± 1.32 |
TUB-β |
6.09 ± 1.47 |
15 |
60S-13 |
4.45 ± 0.97 |
PP2C59 |
3.99 ± 0.81 |
60S-18 |
4.79 ± 1.00 |
60S-13 |
8.83 ± 2.03 |
UBC5B |
2.23 ± 0.41 |
UBC5B |
1.93 ± 0.35 |
UBC17 |
7.03 ± 1.55 |
HIS6 |
6.61 ± 1.59 |
16 |
TUB-β |
4.75 ± 1.08 |
EF1 |
4.21 ± 0.69 |
60S-13 |
5.07 ± 1.07 |
UBC17 |
8.97 ± 1.81 |
UBC17 |
2.31 ± 0.47 |
PGK |
2.00 ± 0.40 |
PP2C59 |
7.04 ± 1.48 |
EF2 |
6.82 ± 1.38 |
17 |
EF2 |
4.83 ± 0.97 |
60S-18 |
4.48 ± 0.91 |
EF1 |
5.21 ± 0.87 |
TUB-β |
9.33 ± 2.19 |
PP2C59 |
2.51 ± 0.48 |
TUB-α |
2.00 ± 0.42 |
TUB-α |
7.13 ± 1.62 |
60S-13 |
7.61 ± 1.70 |
18 |
60S-18 |
4.84 ± 1.04 |
SAMDC |
4.85 ± 1.04 |
EF2 |
5.38 ± 1.04 |
EF1 |
9.44 ± 1.57 |
HIS1 |
2.52 ± 0.59 |
MUB1 |
2.28 ± 0.49 |
TIP41 |
7.51 ± 1.75 |
PP2C57 |
7.77 ± 1.66 |
19 |
EF1 |
6.03 ± 1.04 |
PP2C57 |
5.49 ± 1.30 |
PP2C57 |
8.84 ± 1.94 |
PP2C57 |
12.30 ± 2.71 |
TIP41 |
3.11 ± 0.71 |
PP2C59 |
3.27 ± 0.64 |
HIS1 |
7.70 ± 1.97 |
EF1 |
8.11 ± 1.43 |
20 |
PPIA26 |
8.36 ± 2.02 |
PPIA26 |
5.99 ± 1.67 |
PPIA26 |
10.94 ± 2.84 |
PPIA26 |
15.21 ± 3.73 |
PPIA26 |
4.58 ± 1.03 |
PPIA26 |
4.46 ± 1.04 |
60S-13 |
10.75 ± 2.70 |
PPIA26 |
9.83 ± 2.37 |
Ranking |
H. robusta -leaves |
H. robusta -young stems |
H. robusta -leaves & young stems |
Different tissues |
4 °C treatment |
MeJA treatment |
PEG6000 treatment |
All samples |
---|---|---|---|---|---|---|---|---|
1 |
PP2C59 |
HIS1 |
UBC17 |
TUB-α |
60S-18 |
UBC17 |
HIS6 |
TUB-α |
2 |
UBC5B |
ACT7 |
UBC5B |
PPIA95 |
TUB-α |
PPIA95 |
MUB1 |
PGK |
3 |
HIS6 |
TIP41 |
ACT7 |
TIP41 |
PPIA95 |
EF1 |
PP2C57 |
UBC5B |
4 |
UBC17 |
PGK |
TUB-α |
EF2 |
ACT7 |
HIS1 |
UBC5B |
UBC17 |
5 |
HIS1 |
PPIA95 |
HIS1 |
PP2C59 |
PGK |
TIP41 |
EF1 |
TIP41 |
6 |
ACT7 |
UBC5B |
PPIA95 |
UBC5B |
PP2C57 |
ACT7 |
PGK |
PP2C59 |
7 |
TUB-α |
MUB1 |
HIS6 |
HIS1 |
UBC17 |
HIS6 |
EF2 |
ACT7 |
8 |
PPIA95 |
EF2 |
PP2C59 |
ACT7 |
EF2 |
EF2 |
SAMDC |
PPIA95 |
9 |
PGK |
UBC17 |
PGK |
EF1 |
SAMDC |
UBC5B |
ACT7 |
HIS1 |
10 |
EF2 |
TUB-β |
TIP41 |
60S-18 |
UBC5B |
PP2C57 |
60S-18 |
EF2 |
11 |
TIP41 |
TUB-α |
MUB1 |
SAMDC |
60S-13 |
PGK |
PPIA95 |
SAMDC |
12 |
MUB1 |
MUB1 |
EF2 |
UBC17 |
HIS1 |
TUB-β |
PPIA26 |
60S-18 |
13 |
PP2C57 |
60S-13 |
60S-13 |
MUB1 |
PP2C59 |
60S-18 |
UBC17 |
MUB1 |
14 |
SAMDC |
EF1 |
60S-18 |
60S-13 |
18S |
TUB-α |
PP2C59 |
TUB-β |
15 |
60S-13 |
PP2C59 |
SAMDC |
TUB-β |
TUB-β |
60S-13 |
18S |
HIS6 |
16 |
18S |
18S |
EF1 |
PGK |
TIP41 |
MUB1 |
TUB-α |
18S |
17 |
EF1 |
HIS6 |
TUB-β |
18S |
EF1 |
18S |
TIP41 |
60S-13 |
18 |
60S-18 |
60S-18 |
18S |
PP2C57 |
HIS6 |
SAMDC |
TUB-β |
EF1 |
19 |
TUB-β |
PP2C57 |
PP2C57 |
HIS6 |
MUB1 |
PP2C59 |
HIS1 |
PP2C57 |
20 |
PPIA26 |
PPIA26 |
PPIA26 |
PPIA26 |
PPIA26 |
PPIA26 |
60S-13 |
PPIA26 |
In order to verify the expression stability of the selected reference genes by the software, expression of TcMYB3 gene was relatively quantified by using the two most stable genes, alone and in combination, and the two most unstable genes in the consensus ranking. Under H. robusta treatment, the relative expression of TcMYB3 in leaves and young stems reached a peak at 12 h when the most stable genes and their combinations were standardized. But when standardized with the most unstable gene, the relative expression of TcMYB3 was abnormally increased (Fig. 5A, B). As shown in Fig. 5C, when using the most stable gene as the reference gene, the expression level of TcMYB3 was increased 1-1.5 times in young leaves compared with the expression level in mature leaves and the expression of other tissues (shoots, young stems, roots and flowers) was down-regulated and the depression multiple is basically the same. But the expression level of TcMYB3 in the flower was the highest using the most unstable gene (PPIA26) as the reference gene. In addition, under other stresses, including 4℃ treatment (Fig. 5D), MeJA treatment (Fig. 5E), and PEG6000 treatment (Fig. 5F), it was found that when the most stable two reference genes and their combination were used for relative quantification of the TcMYB3, the expression level and trends were very similar, whereas when the two most unstable internal reference genes were used for relative quantification, neither the expression level nor trend was consistent. It is evident that the use of unstable references for gene expression analysis in T. ciliata can result in biased results.
It is ideal that reference genes are stably expressed under all experimental conditions and show stable expression levels across various tissues and growth stages of the organism, but such genes are almost non-existent [35]. More and more studies are showing that the genes that are stably expressed in different species and under different conditions change [36–39]. Selection of the most suitable reference gene for specific conditions using RT-PCR is therefore very important. This study was dedicated to discovering the best reference genes for gene expression analysis in T. ciliata under different conditions. There were 20 candidate genes from the T. ciliata transcriptome database screened and analyzed by RT-qPCR. It was found that the best reference genes were not consistent across different conditions. For examples, PP2C59 and UBC5B were most suitable for leaves under H. robusta treatment, whereas HIS1 and ACT7 were more optimal for young stems under H. robusta treatment, TUB-α and PPIA95 for comparing different tissues, and 60S-18 and TUB-α for leaves under 4℃ treatment.
In this study, we used four methods, geNorm, NormFinder, BestKeeper, and RankAggreg, to evaluate the stability of expression of 20 candidate genes. The first three algorithms were used to evaluate the expression stability of candidate genes. Our results demonstrated that reference values and calculation methods used by the three algorithms were very different [40]. NormFinder calculates stability values based on intra- and inter-group differences [25], while geNorm compares a reference gene with all genes in a given sample to evaluate the best reference gene [24]. In BestKeeper, CV and SD values determine the ranking of stability of candidate genes [41]. Due to the difference of the algorithms in the three software packages, they generated different rankings for the same set of experimental data, although the results of analysis with geNorm and NormFinder had few variations in this study. For example, according to geNorm and NormFinder, 18S was the best reference gene to use across different tissue conditions, whereas the best reference gene identified by BestKeeper analysis was HIS6. For young stem tissue exposed to H. robusta stress, ACT7 and UBC5B were put forward by geNorm and NormFinder, while ranking behind BestKeeper. In order to consolidate the results from the three algorithms, RankAggreg was used for overall ranking [34]. Many researchers use ReFinder to calculate the final ranking [42–44]. ReFinder assigns an appropriate weight to each gene and calculates the geometric mean of its weights to give the final ranking [45]. RankAggreg uses a cross-entropy Monte Carlo algorithm or genetic algorithm to produce aggregated ordered lists based on rankings [34]. Both these tools play very important roles in the consolidation the screening results for internal reference genes from other software.
Under H. robusta stress, the reference genes that performed best in leaves and young stem tissues were different in our study. PP2C59 and UBC5B showed high stability of expression in leaves, while only PP2C59 was ranked highly for young stems. Other researchers have also studied the best reference genes for plants under pest stress, STP4 was found to be the best reference gene for use in Brassica juncea under biotic stress caused by aphid infestation [46]. ABCT and FBOX were found to be the most stable in soybean under soybean aphid (SBA) stress; TUB4 and TUA4 were stable under two-spotted spider mite (TSSM) stress [47]. Miranda indicated that both GmELF1A and GmTUA5 were stable reference genes for normalization of expression data obtained from soybean roots infected with Meloidogyne incognita, and GmCYP2 and GmELF1A were the best reference genes in soybean leaves infested by Anticarsia gemmatalis [48]. Once again, the appropriate reference genes for different species under different conditions and in different tissues vary. It is necessary to identify the best reference genes for specific conditions by RT-qPCR.
Under most experimental conditions in this study (all except for PEG6000 treatment), the reference gene with the worst performance was PPIA26, which was the best reference gene recommended by BestKeeper for use under PEG6000 stress. Another gene in this family, PPIA95, it ranked first in geNorm analysis for both 4℃ cold stress and MeJA treatment. For leaves and young stems under H. robusta treatment, and MeJA stress, PPIA95 ranked third in NormFinder analysis. In BestKeeper analysis, PPIA95 ranked third under 4℃ cold stress and PEG6000 induced drought stress, as for all samples it ranked second. Overall, the PPIA gene family is potential for use as a reference gene set in T. ciliate. The PPIA gene family encodes proteins with functions in immune responses, and resistance to cancer, autoimmune diseases, protozoan and viral infections [49]. In plants, genes of the PPIA family are rarely used as internal reference genes, but they are abundantly expressed in the T. ciliata transcriptome data and the expression level in each sample is very similar, which is the main reason for choosing them. As reference genes, they are also stably expressed in animals. For example, in different heart and disease conditions, PPIA is recognized by ReFinder as the best reference gene in different skeletal muscles of mice, and it is ranked first for human endometrial cancer [50, 51]. PPIB is believed as more optimal reference gene in analyzing the blood of Machado-Joseph disease (MJD) patients [52].
This study is the first screen and verification of expression stability analysis of a series of reference genes under different conditions in T. ciliata, showing that the optimal reference genes were TUB-α and PGK across all samples; PP2C59 and UBC5B in leaves and HIS1 and ACT7 in young stems under H. robusta treatment; TUB-α and PPIA95 in different tissues; 60S-18 and TUB-α under 4℃ treatment;UBC17 and PPIA95 under MeJA treatment; HIS1 and MUB1 under PEG6000 treatment. We believe this research is important for accurate quantification and expression analysis of genes under different conditions in T. ciliata. In future, it will play a vital role in the molecular breeding work of T. ciliata, such as the research on the H. robusta-resistant and drought-resistant varieties, as well as the research on the metabolic pathways of precious compounds in plants.
Plant materials
Five different experiments were conducted for data collection (Table 1). Experimental samples were all collected from one-year old T. ciliata, grown in a greenhouse in South China Agricultural University (SCAU). For samples from different tissues, mature leaves, young leaves, flowers, shoots and young stems were collected at 9:00 am on August 25, 2019. Before treatment with the H. robusta, PEG6000, 4℃, and MeJA, all the seedlings were pre-incubated in incubator for seven days with 16 hours of light at 28℃ and 8 hours of dark at 22℃ to mimic the wild environment. For H. robusta treatment, seedlings were exposed to herbivory, and leaves and young stems were harvested after 0, 6, 12, 24 and 36 h. After seedlings were treated with 0%, 5%, 10%, 20% and 30% (w/v), of PEG6000 for seven days, leaves were collected. For 4℃ treatment, seedlings were placed at 4℃, and samples (leaves) were taken at 0, 6, 12, 24and 36 h. Seedlings sprayed with MeJA(100 µM) were sealed in plastic bags and leaves were collected at 0, 6, 12, 24 and 36 h. Three biological replicates were taken for each sampling point, and all samples were immediately frozen in liquid nitrogen and stored at -80℃.
Total RNA was extracted from all samples using a HiPure HP Plant RNA Mini Kit (Magen) with DNase treatment to remove genomic DNA. The quality of the RNA was determined with NanoDrop ND1000 (Thermo Scientific). RNA samples with absorbance ratios A260/A280 and A260/A230 both around 2.0 were selected for further analysis. To synthesize cDNA, 0.5 µg of total RNA was used according to the manufacturer's instructions for the HiScript II Reverse Transcript kit (Vazyme). Five-fold diluted cDNA was used for subsequent RT-qPCR experiments.
Twenty candidate reference genes were selected from the T. ciliata leaf transcriptome database by reviewing previous literature: PGK, 60S-18, 60S-13, HIS1, HIS6, PP2C57, PP2C59, UBC5B, UBC17, ACT7, SAMDC, EF1, EF2, 18S, TUB-α, TUB-β, MUB1, PPIA26, PPIA95, TIP41. Since there is no genomic sequence data available for T. ciliata, we designed primers based on the sequences in the T. ciliata transcriptome database. Firstly, the Coding Sequence (CDS) and genomic DNA sequences (gDNA) of the candidate reference genes were amplified separately, and then the introns and exons of the candidate reference genes were obtained by sequence alignment (NCBI-Blast), and finally primers for RT-qPCR analysis of each reference gene were designed using the web-based Primer-Blast tool from NCBI (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). Apart from 60S-18, SAMDC, TUB-β and PPIA26, RT-qPCR primers were designed across introns. Details of these genes and primers are shown in Table 2.
The volume of each PCR amplification reaction mix was 20 µL, containing 10 µL of Phata Max Buffer, 2 µL of five-fold diluted cDNA, 2 µL of each primer (10 µM), 0.5 µL of dNTP, 0.5 µL of Phata Max Super-Fidelity DNA Polymerase, and 5 µL ddH2O. The PCR reaction procedure was as follows: 95℃ for 3 min, 35 cycles of 95℃ for 15 s, 55℃ for 15 s, 72℃ for 15 s, followed by 5 min extension at 72 °C. The RT-qPCR reaction mixture consisted of 10 µL ChamQ Universal SYBR qPCR Master Mix (Vazyme), 2 µL of cDNA, 0.4 µL of each primer (10 µM), and 7.2 µL ddH2O to a final volume of 20 µL. and it was performed on LightCycler480 (Roche Molecular Biochemicals, Mannheim, Germany) with optical 96-well plate. To test the specificity of the RT-qPCR primers, the products of PCR were analyzed by nucleic acid electrophoresis on a 2% (w/v) gel and the melting curve was included after amplification. All samples used for RT-qPCR analysis had three biological replicates, each containing three technical replicates. In order to calculate the gene specific PCR efficiency (E) and correlation coefficient (R2) of each gene, a standard curve was generated from the mixed complementary DNA (cDNA) using a fivefold dilution series.
CT values were obtained by RT-qPCR, and used to evaluate the expression levels of candidate genes in different experimental conditions and tissues. Three commonly employed algorithms, geNorm [24], NormFinder [25], and BestKeeper [26], were used to evaluate the stability of candidate reference genes in different experiments.
The package geNorm (Version3.5) screens stable reference genes by calculating the M value of the stability of each candidate gene, and the criterion is that the smaller the M value is, the better the stability of the candidate gene is. It also calculates pairwise variations of the normalized factor after introducing a new internal reference gene, and determines the number of optimal internal reference genes based on the ratio Vn / Vn+1. If the value of Vn / Vn+1 is less than 0.15, the number of optimal internal reference genes is n. If the value of Vn / Vn+1 is greater than 0.15, the number of optimal reference genes is n + 1. NormFinder selects the most suitable internal reference gene by calculating a stability value for gene expression of the candidates. The lower the stability value is, the more stable the gene is. Using BestKeeper, the SD and CV of expression of each candidate gene can be obtained. The CV ± SD values of different genes then compared to determine the relative stabilities of expression of the candidates. Finally, in order to generate an overall ranking of candidate genes from the data generated by geNorm, NormFinder and BestKeeper, we used the RankAggreg (version 0.6.5) software package in R as previously described [14, 53–55]. RankAggreg is an algorithmic package that can combine different ranking lists. Based on the size of the rankings list, we used the Cross-Entropy Monte Carlo algorithm [34]. The rankings list previously generated by the three packages were used as input with the following parameters: the distance was calculated using Spearman's Footrule function, rho was set at 0.1, the seed at 100, and the “convIn” argument at 50.
In order to verify the accuracy of the rankings and the stability of expression of the selected reference genes, the two most stable reference genes, alone and in combination, and the two most unstable reference genes, as recommended by RankAggreg were used to verify the relative expression of TcMYB3 under 4℃, MeJA, PEG6000, in different tissues and under H. robusta treatment (leaves and young stems). Finally, we used the 2−△△CT method to calculate the relative expression levels of the verified genes, where △CT = CT (target gene)-CT (reference gene), △△CT = △CT(treatment)-△ CT (control), 2−△△CT = relative expression. Three technical replicates were performed for each biological sample [56].
Quantitative real time PCR; T. ciliate:Toona ciliata Roem; H. robusta:Hypsipyla robusta; MeJA:methyl jasmonate; CT:cycle threshold; TcMYB3:transcription factor MYB3; CV:coefficient of variation; SD:standard deviation; CDS:Coding Sequence; gDNA:genomic DNA; cDNA:complementary DNA.
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Availability of data and materials
The data that support the results are included within the article and its additional files. Other relevant materials are available from the corresponding authors on reasonable request.
Competing interests
The authors declare that they have no competing interests.
Funding
This research is supported by grants from the National Key Research Projects, Forestry Resource Cultivation and Utilization Technology Innovation (Grant No.2016YFD0600606), the Natural Science Foundation of Guangdong Province of China (Grant No.2018A030313798) and Characteristic innovation projects of department of education of Guangdong province(Grant No.2019KTSCX017).
Authors’ contributions
PL designed the research; HS conducted the experiment; QQ, ZD and WM analyzed the data; PL and HS wrote the manuscript. All authors read and approved the manuscript.
Acknowledgments
We thank Guangdong Key Laboratory for Innovative Development and Utilization of Forest Plant Germplasm, for providing laboratory apparatus.