Expression profile of the candidate reference genes
The mean Ct values of the 10 candidate reference genes ranged from 9.95-33.16, and most were distributed between 18 and 24. PP2A presented the highest expression level, with the lowest Ct values. TIP41 had the least variability, with the lowest SD value of 0.72, followed by YLS8, CYP2, PP2A, and ACT, suggesting that these genes possessed stable expression levels under different stress conditions. In contrast, PP2A2 showed the highest mean Ct value (33.16 ± 2.44), which indicated that this reference was not worthy of consideration (Fig. 1).
Expression stability of the reference genes
BestKeeper analysis
The results of BestKeeper analysis are shown in Table 2. As shown in the table, the CV ± SD value increases form top to bottom (the lowest value is shown in the first row of the table). In all groups, TIP41 showed the highest stability with the lowest CV ± SD value of 2.45 ± 0.56, followed by CYP2 and YLS8. Accordingly, under H2O2, NaCl, CuSO4, high temperature, UV, and wild type stress, TIP41 showed the most stable expression level, with CV ± SD values as low as 2.35 ± 0.54, 1.54 ± 0.36, 0.95 ± 0.21, 1.82 ± 0.40, and 1.87 ± 0.45, respectively. With the lowest CV ± SD value of 1.64 ± 0.31 and 1.55 ± 0.29, respectively, YLS8 and GAPDH showed the highest expression stability under low-temperature and PEG stress. It is worth noting that the high CV ± SD value of PP2A2under each condition indicated that its gene expression was not stable, which corresponded with the result of the Ct value. However, PP2A showed the lowest stability with the lowest CV ± SD value under most stress treatments, which was contrary to the outcomeof the Ct value. Therefore, further evaluation is of great significance to obtain accurate results in reference gene selection.
Table 2
Expression stability rank of 10 candidate reference genes using BestKeeper in P. corylifolia (coefficient of variation (CV) ± standard deviation (SD)).
Rank
|
H2O2
|
NaCl
|
CuSO4
|
high-
temperature
|
low-
temperature
|
PEG
|
UV
|
WT
|
Total
|
1
|
TIP41
|
TIP41
|
TIP41
|
TIP41
|
YLS8
|
GAPDH
|
TIP41
|
TIP41
|
TIP41
|
2.35 ± 0.54
|
1.54 ± 0.36
|
0.95 ± 0.21
|
1.82 ± 0.40
|
1.64 ± 0.31
|
1.55 ± 0.29
|
1.87 ± 0.45
|
1.86 ± 0.42
|
2.45 ± 0.56
|
2
|
CYP2
|
PTBP
|
EF-1α
|
EF-1α
|
PTBP
|
TIP41
|
GAPDH
|
PTBP
|
CYP2
|
3.03 ± 0.68
|
1.93 ± 0.41
|
1.58 ± 0.34
|
1.88 ± 0.39
|
1.74 ± 0.37
|
1.69 ± 0.39
|
2.23 ± 0.48
|
2.24 ± 0.48
|
3.23 ± 0.73
|
3
|
EF-1α
|
CYP2
|
PTBP
|
GAPDH
|
EF-1α
|
EF-1α
|
CYP2
|
EF-1α
|
YLS8
|
3.38 ± 0.74
|
1.97 ± 0.44
|
1.99 ± 0.42
|
2.03 ± 0.40
|
2.26 ± 0.49
|
1.99 ± 0.43
|
2.23 ± 0.54
|
2.66 ± 0.55
|
3.64 ± 0.69
|
4
|
YLS8
|
YLS8
|
YLS8
|
ACT
|
TIP41
|
CYP2
|
ACT
|
CYP2
|
PTBP
|
3.69 ± 0.70
|
2.32 ± 0.43
|
2.07 ± 0.39
|
2.25 ± 0.45
|
2.30 ± 0.52
|
2.02 ± 0.46
|
3.12 ± 0.69
|
2.73 ± 0.60
|
3.71 ± 0.80
|
5
|
ACT
|
ACT
|
GAPDH
|
PTBP
|
CYP2
|
YLS8
|
PTBP
|
YLS8
|
EF-1α
|
3.75 ± 0.71
|
3.23 ± 0.61
|
2.23 ± 0.36
|
2.32 ± 0.52
|
2.43 ± 0.55
|
2.02 ± 0.38
|
3.35 ± 0.77
|
3.09 ± 0.56
|
3.96 ± 0.85
|
6
|
PTBP
|
GAPDH
|
CYP2
|
CYP2
|
ACT
|
ACT
|
EF-1α
|
ACT
|
ACT
|
3.88 ± 0.83
|
3.67 ± 0.62
|
2.31 ± 0.51
|
2.98 ± 0.65
|
2.85 ± 0.55
|
2.15 ± 0.43
|
3.89 ± 0.92
|
3.15 ± 0.61
|
4.94 ± 0.97
|
7
|
GAPDH
|
EF-1α
|
ACT
|
YLS8
|
GAPDH
|
PP2A2
|
YLS8
|
PP2A2
|
PP2A2
|
4.24 ± 0.70
|
4.28 ± 0.89
|
2.65 ± 0.49
|
3.49 ± 0.64
|
3.70 ± 0.67
|
6.29 ± 2.05
|
4.61 ± 0.94
|
4.00 ± 1.33
|
5.66 ± 1.88
|
8
|
PP2A2
|
PP2A2
|
NCBP2
|
NCBP2
|
PP2A2
|
NCBP2
|
PP2A
|
NCBP2
|
NCBP2
|
4.80 ± 1.62
|
4.35 ± 1.44
|
5.88 ± 1.35
|
5.09 ± 1.20
|
5.52 ± 1.82
|
9.16 ± 2.26
|
5.08 ± 0.49
|
4.80 ± 1.10
|
7.14 ± 1.70
|
9
|
NCBP2
|
NCBP2
|
PP2A2
|
PP2A
|
PP2A
|
PTBP
|
NCBP2
|
GAPDH
|
GAPDH
|
6.02 ± 1.43
|
7.67 ± 1.77
|
6.20 ± 2.06
|
5.91 ± 0.59
|
5.67 ± 0.56
|
9.42 ± 1.97
|
5.84 ± 1.52
|
5.16 ± 0.89
|
8.04 ± 1.45
|
10
|
PP2A
|
PP2A
|
PP2A
|
PP2A2
|
NCBP2
|
PP2A
|
PP2A2
|
PP2A
|
PP2A
|
7.97 ± 0.81
|
8.82 ± 0.90
|
10.67 ± 1.03
|
6.70 ± 2.22
|
6.25 ± 1.47
|
10.35 ± 1.03
|
6.22 ± 2.06
|
10.71 ± 1.07
|
8.19 ± 0.82
|
NormFinder analysis
Similar to that in BestKeeperanalysis, a lower Ct value indicates a higher stability of gene expression. As shown in Table 3, EF-1α and GAPDH displayed the highest expression stability, showing the lowest Ct values under H2O2 and PEG stress. ACT showed the highest expression stability in the NaCl and UV stress subsets. In the CuSO4 and Cold subsets, YLS8 was the most stable reference gene. TIP41 had the highest stability in the high temperature and wild type subsets. NCPB2 and PP2A2 showed the lowest expression stability under the different stress conditions. Among all sample subsets, YLS8 was the most stable reference gene, followed by TIP41 and CYP2. Interestingly, in the overall ranking, the top three genes from the BestKeeper and NormFinder analyses were consistent.
Table 3
Expression stability rank of 10 candidate reference genes using NormFinder in P. corylifolia.
Rank
|
H2O2
|
NaCl
|
CuSO4
|
high-
temperature
|
low-
temperature
|
PEG
|
UV
|
WT
|
Total
|
1
|
EF-1α
|
ACT
|
YLS8
|
TIP41
|
YLS8
|
GAPDH
|
ACT
|
TIP41
|
YLS8
|
|
0.090
|
0.072
|
0.075
|
0.083
|
0.080
|
0.045
|
0.039
|
0.034
|
0.242
|
2
|
YLS8
|
GAPDH
|
PTBP
|
ACT
|
PTBP
|
YLS8
|
PTBP
|
YLS8
|
TIP41
|
|
0.102
|
0.118
|
0.075
|
0.117
|
0.080
|
0.069
|
0.257
|
0.099
|
0.348
|
3
|
GAPDH
|
YLS8
|
GAPDH
|
GAPDH
|
TIP41
|
ACT
|
GAPDH
|
CYP2
|
CYP2
|
|
0.232
|
0.167
|
0.115
|
0.137
|
0.122
|
0.080
|
0.316
|
0.278
|
0.365
|
4
|
TIP41
|
TIP41
|
TIP41
|
EF-1α
|
EF-1α
|
EF-1α
|
YLS8
|
PTBP
|
ACT
|
|
0.289
|
0.225
|
0.115
|
0.270
|
0.227
|
0.170
|
0.415
|
0.210
|
0.420
|
5
|
ACT
|
PTBP
|
EF-1α
|
YLS8
|
ACT
|
CYP2
|
EF-1α
|
CYP2
|
EF-1α
|
|
0.342
|
0.333
|
0.135
|
0.289
|
0.334
|
0.274
|
0.427
|
0.278
|
0.420
|
6
|
CYP2
|
CYP2
|
CYP2
|
PTBP
|
CYP2
|
TIP41
|
TIP41
|
EF-1α
|
PTBP
|
|
0.365
|
0.358
|
0.194
|
0.320
|
0.360
|
0.408
|
0.466
|
0.443
|
0.778
|
7
|
PTBP
|
EF-1α
|
ACT
|
CYP2
|
GAPDH
|
PP2A
|
CYP2
|
GAPDH
|
GAPDH
|
|
0.382
|
0.367
|
0.329
|
0.524
|
0.413
|
1.187
|
0.480
|
0.599
|
0.928
|
8
|
PP2A
|
PP2A
|
PP2A
|
PP2A
|
PP2A
|
PP2A2
|
PP2A
|
NCBP2
|
PP2A
|
|
0.747
|
0.543
|
0.901
|
0.790
|
0.599
|
1.624
|
0.680
|
1.042
|
1.006
|
9
|
NCBP2
|
PP2A2
|
NCBP2
|
NCBP2
|
NCBP2
|
PTBP
|
NCBP2
|
PP2A
|
NCBP2
|
|
1.390
|
1.282
|
1.561
|
1.135
|
1.233
|
2.121
|
1.367
|
1.335
|
1.436
|
10
|
PP2A2
|
NCBP2
|
PP2A2
|
PP2A2
|
PP2A2
|
NCBP2
|
PP2A2
|
PP2A2
|
PP2A2
|
|
1.668
|
1.420
|
2.302
|
2.161
|
1.679
|
2.330
|
1.431
|
1.614
|
1.745
|
geNorm analysis
As shown in Fig. 2, a lower average M valueindicates a higher gene expression stability. Moreover, a gene with M value > 1.5, was considered unsuitable as a reference gene. In all treatment subsets, the top three reference genes all showed excellent stability with low M values. YLS8 was the most stable reference gene with the lowest M values in all treatment subsets, except for the high temperature subset. PTBP also presented highest expression stability under all stress coditions, except for PEG stress. TIP41 showed the least fluctuation of M value in all stress subsets, and maintained the highest gene expression stability in the total and control subsets. On the contrary, PP2A, PP2A2, and NCBP2 had the lowest gene expression stability under all stress conditions, indicating that they are the least suitable as reference genes. This result is consistent with that of BestKeeper and NormFinder analysis.
In addition to comparing and ranking the expression stability of the reference genes, geNorm can also optimize the number of reference genes by using pairwise variation (Vn/Vn+1) as a parameter. A threshold value of 0.15 was used to determine whether one more reference gene should be added for normalization. A value of < 0.15 indicates that an additional reference gene has no notable significance for normalization. In this study, we calculated the Vn/n+1 value by using geNorm (Fig. 3 and Table S3). In all sterss condition except low temperature, two reference genes were sufficient for normalization, with V2/3 < 0.15. In the low-temperature subset, the V3/4 value was below 0.15, which suggested that three reference genes should be applied for normalization.
Comprehensive analysis of reference gene validation
RefFinder was used to systematically select the optimal reference genes under different treatment conditions. By comparing the results of the above three algorithms with the ΔCt method, RefFinder can assess and screen out the most stable reference genes [40]. The results are listed in Table S4. In addition, all gene pairs were compared with each other and ranked according to the value of ΔCt, from low to high, as shown in Table S5. According to these results, TIP41 showed high stability in most treatment groups, which was consistent with the results of geNorm analysis.
Next, the geometric means of the three algorithms were calculated according to the ranks of the reference genes. TIP41, CYP2, and YLS8, ranked top three in all algorithms. In contrast, PP2A, NCBP2, and PP2A2, which suggested that these genes were not suitable for use as reference genes under different stress conditions (Fig. 4). Table S6 shows the ranking of the reference genes in all four algorithm tools revealing the genes that had stable expression under different stress conditions.