Evaluation of Angelica decursiva reference genes under various stimuli for RT-qPCR data normalization

Selection of candidate reference genes, evaluation of amplification specificity and PCR efficiency

Here, 11 candidate reference genes were selected based on the transcriptome data of A. decursiva (unpublished). Table 1 lists all the candidate reference gene names and abbreviations, homologous genes from Arabidopsis, primer sequences, amplification length, annealing temperature (°C), PCR efficiency (E), and correlation coefficient (R²). Next, conventional PCR and RT-qPCR were used to verify the primer-specific amplification of all candidate reference genes. As illustrated in Fig. S1, based on agarose gel electrophoresis and melting curve analysis which showed a peak, the 11 primer pairs were amplified by a single specific amplicon. Here, the amplification efficiency range of CYP2 and PP2A was 1.711 and 1.880 respectively whereas the correlation coefficient range of UBQ10 was 0.888, and that of PP2A, PTBP, TIP41, and ACT was 1.000.

Expression profile of candidate reference genes

Here, the cycle threshold values (Cp) showed the number of cycles when the generated fluorescent signal reached a level that could be detected. Therefore, in this study, the expression profile of the candidate reference genes were directly determined using the obtained Cp values. As shown in Fig. 1, the average Cp values of these 11 reference genes are distributed between 10 and 30, and a majority of them are distributed between 23 and 27. Among all the candidate reference genes, NCBP2 had the lowest average Cp values, whereas CYP2, PTBP, and TUBA had the highest average Cp values. Moreover, each reference gene had a different coefficient of variation under different conditions. Also, it was observed that SAND and PTBP had the lowest variability, whereas TUBA, which had the highest Cp value, had the highest variability. This value ranged between 25 and 32. Therefore, it may not be qualified as a reference gene. Through the Cp value combined with box-plot not only display the expression profile of the reference gene, but also give us a glimpse in their stability (Table S2). However, considering the complexity of their surroundings, we need to check the proper use of the references in every particular experimental condition. Thus, more statistical tools and further analyses are required.

The RT-qPCR Cp values for 11 candidate reference genes in all samples. The expression data is shown as the Cp value of each reference gene in the samples of A. decursiva. Boxes indicate the 25th/75th percentiles, the lines represent the median, squares represent the means and whiskers represent the maximum and minimum values.

Expression stability of candidate reference genes

Based on the relative expression levels of the 11 candidate reference genes, four algorithms such as BestKeeper, geNorm, NormFinder, and Delta Ct were used to examine their expression stability. Subsequently, the RefFinder tool was employed to sequence the expression stability of all these candidate reference genes and select the most suitable ones.

geNorm analysis: The geNorm analysis evaluated the stability of all the 11 candidate reference genes using the M value (reference expression stability measure). These M values were calculated from the average variation of the gene relative against other candidate reference genes, and the lower M values indicated a higher gene expression stability²¹. As illustrated in Fig. 2, all the candidate reference genes had different levels of stability under different treatments. Here, the M value of PP2A and SAND was the lowest in most treatments and was deliberated as the most stable reference genes. Besides, ACT showed good stability under H₂O₂, MeJA, NaCl, and UV treatments, whereas CYP2 was one of the most stable reference genes in the CuSO₄ and UV treatment groups. On the other hand, the stability of NCBP2 seemed to be less satisfactory than other candidate reference genes. Subsequently, the geNorm algorithm can also determine the optimal number of normalized reference genes²¹ by calculating pairwise mutations (V_n/n+1). Normally, the ideal paired variation (V) score is less than 0.15, which means that the addition of any other gene will not have a substantial impact on standardization. In our study, in the NaCl, Cold, UV, H₂O₂, CuSO₄, and WT subsets, their pairwise variation values (V _2/3) were all less than 0.15. This showed that the addition of a third reference gene lacked significant effects on the normalization of the target gene, thus, the number of the optimal reference genes that were determined was two. In contrast, as shown in Fig. 3 and Table S3, the pairwise variation values (V _3/4) of the PEG and MeJA subsets, was also less than 0.15. This indicated that the number of reference genes that were most suitable for these two treatments was three.

Average expression stability values (M) of the 11 candidate reference genes using geNorm software. Expression stability was evaluated in samples from A. decursiva was submitted to cold stress, drought stress (20% PEG), Methyl jasmonate (MeJA) stress, salt stress (0.5 M NaCl), oxidative stress (H₂O₂), ultraviolet (UV) induction, Metal stress (CuSO₄), untreated (WT) and ‘total’ (all treated samples). The least stable genes are on the left with higher M-value and the most stable genes on the right with lower M-value.

Determination of the optimal numbers of reference genes for normalization by pairwise variation (V_n/n+1). Different treatments are marked as square frame with different colors. The cut-off value is 0.15 and used to determine the optimal number of candidate reference genes for qRT-PCR normalization.

NormFinder analysis: To normalize raw data and measure the expression stability of candidate reference genes through intra- and inter-group variation, the NormFinder analysis uses the 2^−ΔCt method. Like the geNorm analysis, its lower values show higher stability. Table 2 shows the ranking of all candidate genes as calculated using the NormFinder algorithm. Among the PEG, Cold, and ‘total’ treatment subsets, it was observed that SAND was the most stable candidate gene, and also ranked higher in other subsets. After MeJA treatment (0.078) on all sample subsets, TUBA was deliberated as the most stable candidate gene. Besides, among the other five groups (NaCl, UV, H₂O₂, CuSO₄, and WT), PTBP (0.670), TIP41 (0.307), GAPDH (0.495), CYP2 (0.483), and UBQ10 (0.451) were observed to be the most stable genes, respectively. Similar to the geNorm analysis, our study found that the NCBP2 was the most unstable candidate gene.

BestKeeper analysis: Here, the raw data of the Ct value was directly calculated and the more stable candidate reference genes was evaluated by calculating the standard deviation (SD) and the Ct value²². Notably, lower SD and Ct values indicated higher gene expression stability, particularly when the SD was greater than 1 which indicated that the reference gene was unstable²³ and could not be utilized for normalization. Table 3 shows that in the NaCl, MeJA, UV, WT, and ‘total’ treatment subsets, the PTBP gene was considered the most stable. However, the MeJA and ‘total’ treatment subsets were postulated to be expression-insensitive since their SD values were exceed 1 and could not be used for normalization. For this reason, such values should be excluded. Besides, under H₂O₂, CuSO₄, Cold, and PEG treatments, the most stable genes were observed to be CYP2, SAND, TUBA, and PP2A. Generally, as illustrated in Table 3, it was observed that the NCBP2 gene was still the most unstable under all treatments.

Delta Ct analysis: This method evaluated gene expression stability by calculating the mean standard deviation (SD) of each gene. Here, the smaller the value, the higher the stability²⁴. As shown in Table S4, the results of this analysis are consistent with the geNorm analysis. The only difference is that in the Cold and H₂O₂ subsets, TUBA and PP2A candidate genes are the most stable respectively. According to the geNorm analysis, the two most stable candidate genes are SAND and ACT. Hence, SAND and PP2A are the most qualified reference genes.

RefFinder analysis: As shown in Fig. 4, we further calculated the geometric mean of the ranking of each candidate gene using the RefFinder algorithm (http://150.216.56.64/referencegene.php#). This was based on the results obtained from the three statistical algorithms such as geNorm, NormFinder, and BestKeeper. Tables S5 and S6, show the comprehensive index ranking, whereby, the smaller the index, the more stable the gene expression¹⁹. This study showed that SAND and PP2A ranked the highest in most subsets, whereas NCBP2 and TUBA ranked the lowest, making them the most unstable reference genes. In contrast, in the MeJA and Cold subsets, TUBA seemed to be a relatively stable reference gene. Despite the different assessment methods, this resulting difference is reasonable and acceptable. In summary, the stability of these 11 candidate reference genes from the highest to the lowest is: SAND, PP2A, PTBP, ACT, CYP2, EXP-1, GAPDH, UBQ10, TIP41, TUBA, and NCBP2. These results were similar to those obtained from the geNorm and NormFinder analysis, but slightly different from those of the BestKeeper analysis.

Comprehensive ranking of BestKeeper, NormFinder, and geNorm. Using BestKeeper, NormFinder, and geNorm algorithms, the 11 candidate reference genes were ranked comprehensively, and the geometric mean of each reference gene was calculated.