In this study, 4 publically available algorithms were used to evaluate HKGs for higher-accuracy stability rankings: geNorm, NormFinder, BestKeeper, and ΔCt method.
geNorm analysis
Gene expression stability was determined by M-value in geNorm analysis; the lower M value suggests a higher gene expression stability. For group 1, the two most stable genes were EIF4H and EEF2 with the lowest M value, and GAPDH was the most unstable gene (Fig. 3A). For group 2, the two most stable genes were EIF4H and PTPRA, and ACTB was the most unstable gene (Fig. 3B). For group 3, the two most stable genes were EIF4H and PTPRA, whereas ACTB was the most unstable gene (Fig. 3C). For group 4, the two most stable genes were NCBP3 and PTPRA, and GAPDH was the most unstable gene (Fig. 3D). For all samples, geNorm analysis was conducted on 39 samples and 12 HKGs. It was determined that the 3 most stable genes were PTPRA, EIF4H, and NCBP3. Conversely, ACTB, CNBP, and GAPDH were the most unstable genes (Fig. 3E).
geNorm can be used to determine the minimum optimal number of HKGs needed for accurate normalization under different experimental treatments by analyzing pairwise variation (Vn/Vn+1). This method recognizes Vn/Vn+1 <0.15 as a threshold value, and “n” as an appropriate number of the HKG needed. The V2/3 values for all of the experimental variables were below the cut-off value of 0.15 (0.067, 0.078, 0.099, 0.091 and 0.081 for group 1, group 2, group 3, group 4 and all samples, respectively), which indicate that using the double HKGs (first two genes in each group) are sufficiently accurate for use in normalizing qRT-PCR derived gene expression data (Additional file 1: Figure S2). The triplet or more gene combinations can also be used as Vn/Vn+1 <0.15.
ΔCt analysis
The 12 candidate HKGs were analyzed using the Delta Ct method, the data of which is presented in Table 4. The stability of the gene is inversely related to the std-value, thus a lower value indicates greater stability. In group 1, the two most stably expressed genes were PTPRA and SDHA, and the lowest were GAPDH and ACTB. In group 2, the two most stable genes were EEF2 and SDHA, and the least were ACTB and CNBP. In group 3, SDHA and PTPRA were the most stably expressed, whereas ACTB and CNBP were the least. In group 4, the top two stably expressed genes were NCBP3 and EEF2, whereas CNBP and GAPDH were the least. In all samples, the 3 most stable genes were PTPRA, EEF2, and SDHA, while GAPDH, CNBP, and ACTB were the least stable genes.
Table 4
Gene expression stability calculated by the ΔCt method
Gene name | Group 1 | Group 2 | Group 3 | Group 4 | All samples |
PTPRA | 0.391(1) | 0.439(4) | 0.573(2) | 0.443(3) | 0.499(1) |
EEF2 | 0.423(4) | 0.412(1) | 0.587(4) | 0.428(2) | 0.500(2) |
SDHA | 0.391(2) | 0.417(2) | 0.535(1) | 0.475(4) | 0.503(3) |
NCBP3 | 0.457(6) | 0.424(3) | 0.643(7) | 0.422(1) | 0.512(4) |
EIF4H | 0.392(3) | 0.441(5) | 0.604(5) | 0.483(6) | 0.520(5) |
CTBP2 | 0.486(8) | 0.516(7) | 0.619(6) | 0.493(8) | 0.553(6) |
YWHAZ | 0.425(5) | 0.549(8) | 0.583(3) | 0.486(7) | 0.557(7) |
RRAGA | 0.573(10) | 0.442(6) | 0.672(8) | 0.479(5) | 0.568(8) |
SRP68 | 0.460(7) | 0.583(9) | 0.702(10) | 0.499(9) | 0.590(9) |
GAPDH | 0.740(12) | 0.594(10) | 0.699(9) | 0.731(12) | 0.741(10) |
CNBP | 0.555(9) | 0.690(11) | 0.957(11) | 0.644(11) | 0.757(11) |
ACTB | 0.582(11) | 0.826(12) | 1.198(12) | 0.597(10) | 0.973(12) |
A comprehensive ranking of the four methods examined
Next, the ComprFinder algorithm was employed to obtain a comprehensive score by which to rank the potential HKGs. The results are presented in Table 5. In group 1, the 3 most stable HKGs were EIF4H, PTPRA, and SDHA. In group 2, SDHA, NCBP3, and EEF2 were the most stable HKGs analyzed. In group 3, SDHA, PTPRA, and EIF4H were the three most stable HKGs analyzed, whereas NCBP3, PTPRA, and EEF2 were the most stable genes from group 4. The overall rankings, from the highest to the lowest stability, were NCBP3 > SDHA > PTPRA > EEF2 > EIF4H > SRP68 > CTBP2 > YWHAZ > RRAGA > GAPDH > CNBP > ACTB. It is interesting to note that the top 3 of different group rankings have at least 2 of NCBP3, SDHA and PTPRA. In contrast, the commonly used HKGs, ACTB and GAPDH, were relegated to the bottom 2 and 4 positions, respectively.
Table 5
Comprehensive rankings calculated using the ComprFinder method
Ranking No | Group 1 | Group 2 | Group 3 | Group 4 | All samples |
Gene | Score | Gene | Score | Gene | Score | Gene | Score | Gene | Score |
1 | EIF4H | 0.063 | SDHA | 0.059 | SDHA | 0.129 | NCBP3 | 0.105 | NCBP3 | 0.096 |
2 | PTPRA | 0.090 | NCBP3 | 0.082 | PTPRA | 0.170 | PTPRA | 0.105 | SDHA | 0.099 |
3 | SDHA | 0.093 | EEF2 | 0.090 | EIF4H | 0.180 | EEF2 | 0.193 | PTPRA | 0.108 |
4 | EEF2 | 0.171 | EIF4H | 0.210 | SRP68 | 0.230 | SDHA | 0.211 | EEF2 | 0.129 |
5 | SRP68 | 0.174 | RRAGA | 0.227 | EEF2 | 0.247 | SRP68 | 0.245 | EIF4H | 0.143 |
6 | YWHAZ | 0.256 | PTPRA | 0.236 | NCBP3 | 0.252 | CTBP2 | 0.263 | SRP68 | 0.192 |
7 | NCBP3 | 0.282 | CTBP2 | 0.322 | YWHAZ | 0.277 | EIF4H | 0.309 | CTBP2 | 0.248 |
8 | CTBP2 | 0.358 | SRP68 | 0.430 | CTBP2 | 0.293 | RRAGA | 0.327 | YWHAZ | 0.311 |
9 | RRAGA | 0.605 | GAPDH | 0.609 | GAPDH | 0.399 | YWHAZ | 0.361 | RRAGA | 0.320 |
10 | CNBP | 0.637 | YWHAZ | 0.630 | RRAGA | 0.404 | ACTB | 0.795 | GAPDH | 0.603 |
11 | ACTB | 0.677 | CNBP | 0.697 | CNBP | 0.730 | CNBP | 0.820 | CNBP | 0.680 |
12 | GAPDH | 0.971 | ACTB | 0.943 | ACTB | 1.000 | GAPDH | 0.994 | ACTB | 1.000 |
NCBP3, SDHA, PTPRA were the 3 most stable HKGs across all samples. Furthermore, their scores were within a tight range, calculated at 0.096, 0.099, and 0.108, respectively. They were also preferably ranked in groups 1–4 relative to other genes. Therefore, they were considered to be the 3 most promising candidate HKGs, and were advanced for further validation.
Validation of HKGs by DKK1, SHH, and FGF5 genes
Based on the above analyses, 3 target genes, including DKK1, SHH, and FGF5 were further characterized based on their changes in expression levels during the secondary hair follicle cycle (T1, T2, T3) with normalizations of different single HKG and multi-gene combinations. It was observed that NCBP3, SDHA, and EEF2 were the top 3 HKGs in group 2 (factor: hair follicle cycle) based on expression stability. Therefore, it can be concluded that the combination of NCBP3 + SDHA + EEF2 was the best-normalized gene set for Group 2. Since these 3 genes (NCBP3, SDHA, and PTPRA) are possibly the most important candidate HKGs, they were further characterized to determine optimal combinations for normalization of gene expression studies. Four multi-gene combinations, including NCBP3 + SDHA + PTPRA, NCBP3 + SDHA, NCBP3 + PTPRA, and SDHA + PTPRA, in addition to 3 single-genes (NCBP3, SDHA, and PTPRA) were used for the analysis. Conversely, ACTB and GAPDH were used for comparison, and were also examined as the multi-gene combination ACTB + GAPDH. In total, 11 multi-gene combinations or single genes were used as the normalization factors. For multiple gene combinations, the geometric average of their Ct value was calculated. The relative gene expression level was calculated as 2−ΔCt, ΔCt = Δ (Cttarget gene-CtHKGs).
As is shown in Fig. 4A, the expression profiles of DKK1 were similarly obtained using the 7 stable single-gene and multi-gene combinations. Furthermore, it was observed that DKK1 was more highly expressed in the T2 period compared to T1, and it was the most highly expressed during the T3 period. Among the unstable single- and multi-gene combinations, only the ACTB and ACTB + GAPDH performed similarly to the stable genes. However, the gene expression profile as normalized by GAPDH different from the other conditions, and no significant difference has been identified among T1, T2 and T3 periods. Expression of the SHH gene was even during the T1 and T2 periods, but a significant decrease in the T3 period compared to T1 and T2 was observed (Fig. 4B). Similarly, the seven stable internal HKG combinations identified this trend, but ACTB did not. The combination of ACTB + GAPDH identified this expression change as a trend, but was not able to detect significant changes in expression. The expression profile of the FGF5 gene, when normalized by the most stable candidate HKGs used individually or in combination here, were very similar. High levels of expression was observed in the T2 period, but no statistical significance expression was identified relative to T1 and T3 (Fig. 4C). The combination of ACTB + GAPDH showed a similar pattern to the stable HKGs, but when ACTB and GAPDH were used individually, the expression patterns were completely different. Furthermore, ACTB also identified significant differences in the T2 period relative to T1.
The above-mentioned results derived from Fig. 4 reflect the differences of expression profiles of a single target gene normalized by 11 types of single or multiple-gene combinations. In order to further understand the relationship of those single or multi-HKG combinations, a correlation analysis on these relative expression data (2−ΔCt) of 3 target genes was performed. As is shown in Fig. 5, the normalized results by NCBP3 + SDHA + EEF2 and NCBP3 + SDHA + PTPRA had a high correlation coefficient (R = 0.990, P < 0.001), suggesting that they have extremely similar normalization capabilities. Other double-gene combinations including NCBP3 + SDHA, NCBP3 + PTPRA, and SDHA + PTPRA had high correlation coefficients, ranging from 0.969–0.997 with NCBP3 + SDHA + EEF2. Also, these double-gene combinations had high correlation coefficients of 0.989–0.994 with NCBP3 + SDHA + PTPRA. This indicated that these 3 types of double-gene combinations exhibited similar normalization capabilities like NCBP3 + SDHA + EEF2 and NCBP3 + SDHA + PTPRA. For single stable HKGs NCBP3, SDHA, and PTPRA also exhibited high correlation coefficients with NCBP3 + SDHA + EEF2 (0.942–0.973) and with NCBP3 + SDHA + PTPRA (0.952–0.977). The ACTB, GAPDH, and ACTB + GAPDH combinations had relatively low correlation coefficients with anyone of stable single- (0.513–0.780) and multi-gene combinations (0.548–0.738).