Expression profiling of 12 candidate reference genes in C. militaris
12 genes (namely ACT, TUB, UBC, EF-1α, GAPDH, PP2A, UBQ, PGK, RPS, FBOX, CYP and GTPB) were selected as candidate RGs to determine the most stable RGs under cold treatments. We first calculated the PCR efficiency of 12 pairs of primers according to the previously reported method [19]. As shown in Table 1, RT-qPCR efficiency of the 12 genes ranged from 95.5% (CYP) to 106.2% (EF-1α) which fell within the acceptable range (80–120%) [20]. The cycle threshold (Ct) values of the 12 genes exhibited a high variation ranging and shown in Fig. 1. The Ct values ranged from 16.6 to 28.4 and 16.9 to 29.0 of short time (Fig. 1A) and long time (Fig. 1B), respectively. RPS and PP2A was the most and least transcribed, respectively, across all the tested samples (Fig. 1A and 1B). Three different algorithms (geNorm, NormFinder and BestKeeper) were used to analyze the expression stability of the 12 genes in the next section.
Table 1
Details on Primers Used for Quantitative Reverse Transcriptase Polymerase Chain Reaction Analysis
Gene | Annotation | Accession No. | Primer Sequences (5’-3’) | Amplification Size (bp) | Efficiency (%) |
ACT | Actin cytoskeleton protein | XM_006669552.1 | F: CAACAACTTCCTGACGGGC R: TCCTTGGGCTTCTGCTGAC | 219 | 99.6 |
TUB | Tubulin | XM_006669203.1 | F: ATGTCGTTCGTCGTGAGG R: AGAGTGGCGTTGTAGGGT | 208 | 97.5 |
UBC | Ubiquitin-conjugating enzyme | XM_006672277.1 | F: ACCATTGACACGAGCCAGTT R: GCCCATGTAAGCCTCCTCA | 204 | 101.3 |
EF-1α | Elongation factor 1-alpha | XM_006665968.1 | F: TATCGGAACTGTGCCTGT R: CGTTACCACGACGGATTT | 200 | 106.2 |
GAPDH | Glyceraldehyde 3-phosphate dehydrogenase | XM_006669697.1 | F: CATCCACTCCTACACTGCTAC R: CTCAAGACGAACAGTCAGGT | 220 | 102.3 |
PP2A | Serine/threonine protein phosphatase 2A | XM_006673033.1 | F: CCTCCTACAGTCGTCATCAGC R: AGAAATGTCAAAGCGAGA | 198 | 105.7 |
UBQ | Polyubiquitin | XM_006672469.1 | F: TCAAAGAAGATAATGGTAACG R: GTATGGGTTCTCGGAAAGGT | 209 | 102.5 |
PGK | Phosphoglycerate kinase | XM_006673408.1 | F: GCTCAAGCCCGTCGTTTC R: CCCTCCTCCTCAATGTGG | 159 | 98.9 |
RPS | Ribosomal protein S25 | XM_006665399.1 | F: AAGTGGTCTAAGGGCAAGG R: TTCTCCTCCAGGTCGGTAA | 182 | 97.8 |
FBOX | F-box protein | XM_006672810.1 | F: CCGATGACAACGACAGCGAC R: GTAGTTGACCGTGGAGATGT | 224 | 101.3 |
CYP | cyclophilin | XM_006674778.1 | F: TTTTCCGCCTTATTCCACC R: TCCAGAGCATCAAATCCCT | 197 | 95.5 |
GTPB | GTP-binding protein | XM_006674648.1 | F: TAAGAAGCCCAAGAAGAAAA R: GTCCCACAGGTTCAGCGT | 181 | 103.4 |
F, forward; R, reverse |
geNorm Analysis
The geNorm algorithm calculates the value of measurement (M) to evaluate the stability of gene expression. The lowest M value and highest M value has the most stable and unstable expression, respectively. As determined by geNorm (Fig. 2), expression stability (M values) of the 12 genes both in short and long time were within the acceptable range (M < 1.5). The M value ranged from 0.91 (PGK) to 0.34 (UBC and PP2A) in short time (Fig. 2A) and 1.12 (PGK) to 0.25 (TUB and EF-1α) in long time (Fig. 2B). The genes were ranked from the highest M value (least stable) to the lowest M value (most stable): PGK, ACT, GAPDH, TUB, EF-1α, FBOX, RPS, UBQ, CYP, GTPB, UBC and PP2A in short time (Fig. 2A), PGK, UBC, ACT, FBOX, RPS, GTPB, PP2A, CYP, GAPDH, UBQ, TUB and EF-1α in long time (Fig. 2B).
NormFinder Analysis
Next, the expression stability of the 12 candidate RGs was analyzed using NormFinder to confirm the geNorm results. In short time, UBC and UBQ were the most stable; their stability values were 0.141 and 0.244, respectively (Table 2). In long time, CYP and UBQ were the most highly ranked; their stability values were 0.123 and 0.202, respectively (Table 3). Combining the results of geNorm and NormFinder analysis, UBC followed by PP2A, GTPB and UBQ were the most stable RGs in short time (Table 2 Mean Rank); UBQ followed by TUB, CYP and GAPDH were the most stable RGs in long time (Table 3 Mean Rank).
Table 2
Gene expression stability under short time cold stress ranked by geNorm, NormFinder and BestKeeper
Genes | geNorm | NormFinder | Mean Ranka | BestKeeper |
M Value | Rank Order | Stability Value | Rank Order | CV | SD |
TUB | 0.254 | 1.5 | 0.308 | 4 | 2 | 4.934 | 0.711 |
EF-1α | 0.254 | 1.5 | 0.360 | 6 | 5 | 6.359 | 1.930 |
UBQ | 0.470 | 3 | 0.202 | 2 | 1 | 4.066 | 0.701 |
GAPDH | 0.523 | 4 | 0.247 | 3 | 4 | 4.589 | 0.949 |
CYP | 0.548 | 5 | 0.123 | 1 | 3 | 3.560 | 0.693 |
PP2A | 0.571 | 6 | 0.404 | 7 | 7 | 5.381 | 1.442 |
GTPB | 0.599 | 7 | 0.319 | 5 | 6 | 5.670 | 1.426 |
RPS | 0.655 | 8 | 0.489 | 8 | 8 | 3.882 | 0.704 |
FBOX | 0.701 | 9 | 0.543 | 9 | 9 | 2.487 | 0.644 |
ACT | 0.764 | 10 | 0.641 | 10 | 10 | 2.935 | 0.743 |
UBC | 0.830 | 11 | 0.730 | 11 | 11 | 8.127 | 1.792 |
PGK | 1.121 | 12 | 1.733 | 12 | 12 | 7.325 | 1.786 |
Genes
|
geNorm
|
NormFinder
|
Mean Ranka
|
BestKeeper
|
M Value
|
Rank Order
|
Stability Value
|
Rank Order
|
CV
|
SD
|
TUB
|
0.254
|
1.5
|
0.308
|
4
|
2
|
4.934
|
0.711
|
EF-1α
|
0.254
|
1.5
|
0.360
|
6
|
5
|
6.359
|
1.930
|
UBQ
|
0.470
|
3
|
0.202
|
2
|
1
|
4.066
|
0.701
|
GAPDH
|
0.523
|
4
|
0.247
|
3
|
4
|
4.589
|
0.949
|
CYP
|
0.548
|
5
|
0.123
|
1
|
3
|
3.560
|
0.693
|
PP2A
|
0.571
|
6
|
0.404
|
7
|
7
|
5.381
|
1.442
|
GTPB
|
0.599
|
7
|
0.319
|
5
|
6
|
5.670
|
1.426
|
RPS
|
0.655
|
8
|
0.489
|
8
|
8
|
3.882
|
0.704
|
FBOX
|
0.701
|
9
|
0.543
|
9
|
9
|
2.487
|
0.644
|
ACT
|
0.764
|
10
|
0.641
|
10
|
10
|
2.935
|
0.743
|
UBC
|
0.830
|
11
|
0.730
|
11
|
11
|
8.127
|
1.792
|
PGK
|
1.121
|
12
|
1.733
|
12
|
12
|
7.325
|
1.786
|
Table 3 Gene expression stability under long time cold stress ranked by geNorm, NormFinder and BestKeeper
Optimal number of reference gene for normalization across the experimental sets
Next, we determine the minimal number of genes for RT-qPCR normalization by estimated pairwise variation (Vn/Vn + 1) in geNorm. Below 0.15 which was a proposed cut-off Vn/Vn + 1 value, adding an additional RG is not required [16]. According to this principle, the Vn/Vn + 1 value was calculated. In short time, the V2/3 value was 0.124 (< 0.15), therefore, UBC together with PP2A would be sufficient for purpose (Fig. 3A and Table 2). In long time, compared with V2/3 (0.188), V3/4 value was 0.125 (< 0.15), suggesting that three RGs, UBQ, TUB and CYP were identified as available for normalization (Fig. 3B and Table 3).
BestKeeper Analysis
We further used BestKeeper to analyze the stability of 12 candidate RGs by calculating the SD (standard deviation) and the CV (coefficient of variation) of their Ct values. It has been reported that any studied gene with the SD lower than 1 can be considered stable, and the most stably expressed gene, exhibiting the lowest CV [21]. As shown in Table 2, in short time, the top two stable genes UBC and PP2A with lower CV values 4.501 and 3.596, respectively. And the SD values of UBC and PP2A were 0.882 (< 1) and 0.872 (< 1), respectively. In long time, the top three stable genes UBQ, TUB and CYP with lower CV values 4.066, 4.934 and 3.560, respectively. And the SD values of UBQ, TUB and CYP were 0.701 (< 1), 0711 (< 1) and 0.693 (< 1), respectively (Table 3). These BestKeeper results also suggested that UBC/PP2A and UBQ/TUB/CYP were stable in short and long time, respectively.
Taken together, our results suggesting that two RGs UBC/PP2A and three RGs UBQ/TUB/CYP were identified as available for RT-qPCR normalization in cold treatments for short and long time, respectively, in C. militaris.
Quantitative effects of best and least ranked reference genes on target gene expression
To validate the utility of stable RGs on gene expression analysis, three different strategies, the least stable RG, the best ranked RG and multiple stable RGs were selected to normalize the expression of target gene in short and long time, respectively (PGK, UBC and UBC/PP2A for short time; PGK, UBQ and UBQ/TUB/CYP for long time). HK genes in cold stress response were used as target genes. Base on the annotation with C. militaris genome database (National Center for Biotechnology Information accession GSE28001 and AEVU00000000), 9 genes (CmHK1-9, Table S1) were identified as HK homologous genes. Normalization of the 9 HK genes using the three different strategies showed that a significant increase in transcription of HK2/HK3/HK6 and HK1/HK3/ HK5 were observed under cold stress conditions for short and long time, respectively (Fig. 4). The gene expression levels of other HK genes were not shown. When using the best ranked UBC or the multiple stable UBC/PP2A as the housekeeping genes in short time, the upregulated trends for HK2, HK3 and HK6 genes were consistent and the peak points were observed at 1 h (Fig. 4A-C). When using the best ranked UBQ or the multiple stable UBQ/TUB/CYP as the housekeeping genes in long time, the upregulated trends for HK1, HK3 and HK5 genes were consistent and the peak points were observed at 4 d (Fig. 4D-F), respectively. However, using the least ranked PGK as the housekeeping gene in short time, the peak points were observed at 2, 0.5, 2 h of HK2, HK3, HK6, respectively (Fig. 4A-C). In addition, when PGK was applied as the housekeeping gene in our experiment, the transcription patterns of HK genes were higher than those using the best ranked RG and multiple stable RGs as housekeeping genes (Fig. 4A-F), which indicated that normalization using the least stable PGK as housekeeping gene resulted in an overestimated relative expression level of the target genes. Thus, our experimental results suggest that it is feasible to use any one or more of UBC/PP2A and UBQ/TUB/CYP genes as the normalization gene in cold stress for short and long time, respectively, with C. militaris.