Selection of reference genes and verification of primer specificity
The nucleic acid sequences of the 10 reference genes were obtained from N. tangutorum transcriptome database of our research group, and then their nucleic acid sequences were quickly searched using BLASTX in the Arabidopsis database to determine their homology in N. tangutorum. Based on the unigene sequences, qRT-PCR primers were designed. The amplification efficiency of these 10 pairs of primers ranged from 93.12% to 109.07%, and the correlation coefficients (R2) ranged from 0.9762 to 0.9998 (Table 1). The melting curve only had a single melting peak, and there was no non-specific amplification (see Additional file 1: Figure S1). The results indicated that the primers of these 10 genes were designed reasonably and with good specificity, which met the qRT-PCR standards and allowed for use in subsequent experiments.
Ct value analyses of reference genes
The Ct value is inversely proportional to the transcription level of a gene, that is, a higher Ct value corresponds to lower expression abundance. A change in the Ct value of the same gene in different samples reflects the variation of the expression level. In this study, the Ct values of the 10 genes changed under different abiotic stresses (drought, salt, heat, cold, and ABA) and in different tissues (root, stem, and leaf), and the change trends were different (Fig. 1). Among the 10 reference genes, the lowest and highest Ct values were found in HSP and UBC, which were 19.96 and 32.65, respectively, indicating a wide range of expression abundances. The range of Ct values of ACT, GAPDH, TUA, TUB, CYP, UBC, His, PP2A, HSP, and EF1-α were 24.30–30.08, 21.92–27.53, 23.62–29.77, 23.50–28.80, 25.50–31.79, 25.72–32.65, 22.21–28.48, 25.54–30.54, 19.96–32.06 and 20.70–26.54, respectively. It can be seen that among the 10 genes, the expression abundances of EF1-α, GAPDH, His, ACT, TUB, and TUA were higher and the range of variation smaller. Followed by CYP, and UBC, the PP2A expression level had low variation, but its expression abundance was also low. The HSP expression level was the most variable. It can be preliminarily considered that the expression of EF1-α, GAPDH, His, ACT, TUB, and TUA was relatively stable, and the expression of HSP varied greatly. More accurate results need to be evaluated by reference gene analysis software.
GeNorm analysis
Using geNorm analysis, we obtained the stable expression M value of the reference gene, which is negatively correlated with the stability of the gene. That is, a larger M value represents a more unstable gene. When all the samples of different abiotic stresses and tissues were combined, the M value of His and EF1-α was the lowest and that of HSP was the highest, indicating that the expression of His and EF1-α was the most stable, and the expression of HSP was the least stable. For salt, drought, heat, cold, ABA stress, and different tissues, the top two genes were His and EF1-α, GAPDH and UBC, GAPDH and PP2A, TUA and PP2A, GAPDH and CYP, and GAPDH and His, respectively (Fig. 2).
GeNorm can also calculate the best number of genes when selecting multiple reference genes. It uses the relationship between the threshold V and 0.15 to determine whether an additional reference gene needs to be added. The principle is that when a new reference gene is introduced, the paired variation value (V value) will change. By calculating the ratio Vn/n+1, we can know whether a new reference gene needs to be introduced; 0.15 is the default Vn/n+1 threshold value of the software. If Vn/n+1 > 0.15, the n+1 reference gene needs to be introduced, otherwise it is not needed. In samples under salt, drought, heat, cold, and ABA stress, as well as in different tissues, the V2/V3 values were all less than 0.15, indicating that it was not necessary to introduce a third reference gene to calibrate the qRT-PCR data. However, when analyzing the samples of all experimental conditions, the V2/V3 value was 0.156, which exceeded the threshold of 0.15. The reason may be that the candidate reference genes of N. tangutorum changed significantly under different abiotic stresses (drought, salt, heat, cold, and ABA) and in different tissues (root, stem and leaf). At this time, the use of two reference genes cannot meet the needs of calibrating qRT-PCR data (Fig. 3).
NormFinder analysis
The NormFinder software directly evaluated the stability of reference genes according to the variance within and between groups. Generally, the smaller the stability value, the better the stability of the corresponding gene expression. As shown in Table 2, in all samples, EF1-α, His, and GAPDH had the lowest stability values, that is, 0.068, 0.080, and 0.080, respectively, This showed that among the 10 reference genes, EF1-α, His, and GAPDH had the best stability. HSP had the highest stability value, reaching 0.161, indicating that it had the highest variability. This coincides with the geNorm analysis. The highest ranked pairs of genes under different experimental conditions were not completely consistent. For example, PP2A and HSP, UBC and EF1-α, ACT and PP2A, EF1-α and CYP, GAPDH and His, and EF1-α and TUA were the top ranked genes under salt, drought, heat, cold, ABA stress, and in different tissues, respectively.
Table 2 Expression stability of the ten candidate reference genes calculated by NormFinder
Rank
|
Salt treatment
|
|
Drought treatment
|
|
Heat treatment
|
|
Cold treament
|
|
ABA treatment
|
|
Different tissues
|
|
Total samples
|
Gene
|
Stability
|
|
Gene
|
Stability
|
|
Gene
|
Stability
|
|
Gene
|
Stability
|
|
Gene
|
Stability
|
|
Gene
|
Stability
|
|
Gene
|
Stability
|
1
|
PP2A
|
0.06
|
|
UBC
|
0.033
|
|
ACT
|
0.008
|
|
EF1-α
|
0.018
|
|
GAPDH
|
0.017
|
|
EF1-α
|
0.045
|
|
EF1-α
|
0.068
|
2
|
HSP
|
0.066
|
|
EF1-α
|
0.047
|
|
PP2A
|
0.01
|
|
CYP
|
0.037
|
|
His
|
0.041
|
|
TUA
|
0.053
|
|
His
|
0.08
|
3
|
EF1-α
|
0.09
|
|
GAPDH
|
0.061
|
|
TUA
|
0.024
|
|
GAPDH
|
0.07
|
|
PP2A
|
0.049
|
|
GAPDH
|
0.067
|
|
GAPDH
|
0.08
|
4
|
CYP
|
0.096
|
|
PP2A
|
0.064
|
|
GAPDH
|
0.03
|
|
His
|
0.075
|
|
EF1-α
|
0.049
|
|
TUB
|
0.073
|
|
PP2A
|
0.1
|
5
|
GAPDH
|
0.099
|
|
His
|
0.067
|
|
CYP
|
0.036
|
|
ACT
|
0.077
|
|
CYP
|
0.056
|
|
His
|
0.084
|
|
TUB
|
0.116
|
6
|
TUB
|
0.1
|
|
TUB
|
0.079
|
|
TUB
|
0.037
|
|
UBC
|
0.089
|
|
TUB
|
0.061
|
|
UBC
|
0.114
|
|
CYP
|
0.123
|
7
|
TUA
|
0.102
|
|
ACT
|
0.095
|
|
His
|
0.038
|
|
PP2A
|
0.094
|
|
TUA
|
0.095
|
|
HSP
|
0.136
|
|
ACT
|
0.132
|
8
|
His
|
0.106
|
|
TUA
|
0.142
|
|
EF1-α
|
0.039
|
|
TUA
|
0.117
|
|
ACT
|
0.112
|
|
CYP
|
0.151
|
|
TUA
|
0.136
|
9
|
UBC
|
0.108
|
|
CYP
|
0.156
|
|
HSP
|
0.133
|
|
TUB
|
0.159
|
|
UBC
|
0.115
|
|
ACT
|
0.157
|
|
UBC
|
0.139
|
10
|
ACT
|
0.151
|
|
HSP
|
0.198
|
|
UBC
|
0.139
|
|
HSP
|
0.217
|
|
HSP
|
0.14
|
|
PP2A
|
0.206
|
|
HSP
|
0.161
|
BestKeeper analysis
The BestKeeper algorithm calculates SD and CV using Excel software based on the Ct value of the gene. The smaller the SD and the smaller the CV, the better the stability of the gene. If SD> 1, the gene is considered unstable. The analysis results are shown in Table 3. In all samples, only ACT and PP2A showed SD values <1, which met the requirements of reference genes. However, this is different from the analysis results of geNorm and NormFinder. There were differences in the number of genes that met the SD <1 condition under different experimental conditions. Under the salt, drought, heat, cold, ABA and in different tissues respectively, the highest ranked pair of genes were ACT and PP2A, TUA and ACT, UBC and CYP, TUB and TUA, PP2A and TUB, and CYP and TUB, respectively. Interestingly, three different analysis software programs showed the worst expression stability for HSP.
Table 3 Expression stability of the ten candidate reference genes calculated by BestKeeper analysis
Rank
|
Salt treatment
|
|
Drought treatment
|
|
Heat treatment
|
|
Cold teratment
|
|
ABA treatment
|
|
Different tissues
|
|
Total samples
|
Gene
|
CV±SD
|
|
Gene
|
CV±SD
|
|
Gene
|
CV±SD
|
|
Gene
|
CV±SD
|
|
Gene
|
CV±SD
|
|
Gene
|
CV±SD
|
|
Gene
|
CV±SD
|
1
|
ACT
|
1.46±0.38
|
|
TUA
|
2.10±0.56
|
|
UBC
|
2.6±0.74
|
|
TUB
|
3.18±0.84
|
|
PP2A
|
2.48±0.67
|
|
CYP
|
2.38±0.70
|
|
PP2A
|
3.52±0.98
|
2
|
PP2A
|
2.42±0.65
|
|
ACT
|
2.63±0.70
|
|
CYP
|
2.66±0.77
|
|
TUA
|
3.44±0.97
|
|
TUB
|
3.05±0.75
|
|
TUB
|
2.52±0.68
|
|
ACT
|
3.62±0.96
|
3
|
GAPDH
|
4.41±1.08
|
|
TUB
|
2.75±0.78
|
|
TUB
|
2.72±0.71
|
|
ACT
|
3.66±1.02
|
|
EF1-α
|
3.13±0.68
|
|
PP2A
|
2.53±0.73
|
|
TUB
|
3.84±1.00
|
4
|
TUA
|
4.23±1.11
|
|
PP2A
|
2.94±0.78
|
|
PP2A
|
2.75±0.78
|
|
PP2A
|
3.62±1.04
|
|
His
|
3.43±0.81
|
|
ACT
|
2.61±0.72
|
|
CYP
|
4.19±1.20
|
5
|
TUB
|
4.34±1.12
|
|
EF1-α
|
3.48±0.82
|
|
ACT
|
3.17±0.85
|
|
CYP
|
3.77±1.11
|
|
CYP
|
3.49±0.95
|
|
TUA
|
3.32±0.94
|
|
UBC
|
4.48±1.28
|
6
|
UBC
|
4.28±1.23
|
|
UBC
|
3.37±0.97
|
|
GAPDH
|
3.32±0.83
|
|
His
|
5.67±1.41
|
|
ACT
|
3.54±0.90
|
|
EF1-α
|
4.42±1.08
|
|
TUA
|
4.56±1.21
|
7
|
CYP
|
4.49±1.30
|
|
GAPDH
|
3.53±0.87
|
|
His
|
3.67±0.95
|
|
UBC
|
5.58±1.57
|
|
GAPDH
|
4.15±0.97
|
|
GAPDH
|
4.46±1.15
|
|
GAPDH
|
4.84±1.19
|
8
|
His
|
5.26±1.31
|
|
His
|
4.74±1.19
|
|
TUA
|
3.71±0.96
|
|
GAPDH
|
5.71±1.44
|
|
UBC
|
3.8±1.07
|
|
His
|
4.81±1.27
|
|
EF1-α
|
5.27±1.24
|
9
|
EF1-α
|
5.73±1.36
|
|
CYP
|
4.96±1.41
|
|
EF1-α
|
4.26±1.01
|
|
EF1-α
|
5.78±1.38
|
|
TUA
|
4.02±1.01
|
|
HSP
|
5.74±1.62
|
|
His
|
5.57±1.39
|
10
|
HSP
|
6.50±1.88
|
|
HSP
|
5.91±1.67
|
|
HSP
|
9.99±2.50
|
|
HSP
|
5.74±1.60
|
|
HSP
|
8.29±2.21
|
|
UBC
|
5.77±1.73
|
|
HSP
|
7.92±2.18
|
Reference gene validation
In order to verify the screening results of geNorm and NormFinder, we used the two genes (EF1-α and His) with the best ranking stability and the least stable gene (HSP) as reference genes and analyzed the expression pattern of NtCER7 related to the regulation of waxy synthesis of N. tangutorum under salt stress and in different tissues by qRT-PCR (Fig. 4). When EF1-α and His were used as reference genes under salt stress, the trend of NtCER7 expression showed an inverted V shape from 2 h to 24 h, with a peak at 4 h. When the data of NtCER7 expression were normalized with these two reference gene combinations, consistent results were obtained. When using HSP, the results showed that the lowest value appeared at 4 h, and the trend increased from 4 to 12 h. On the other hand, in different tissues, when EF1-α, His, and EF1-α+His were used to calibrate the expression data of NtCER7, they all showed consistent results. The expression level of NtCER7 showed a trend of stem> leaf> root. However, the expression pattern of NtCER7 normalized by HSP in different tissues showed the trend of being highest in roots, followed by stems, and lowest in leaves. It is clear that the expression of NtCER7 in the stem was underestimated. Therefore, based on the above experimental results, we believe that the results of geNorm and NormFinder were reliable.