Transcriptome-based selection and validation of optimal house-keeping genes for skin research in goats (Capra hircus)

Background In quantitative real-time polymerase chain reaction (qRT-PCR) experiments, accurate and reliable target gene expression results are dependent on optimal amplification of house-keeping genes (HKGs). RNA-seq technology offers a novel approach to detect new HKGs with improved stability. Goat (Capra hircus) is an economically important livestock species and plays an indispensable role in the world animal fiber and meat industry. Unfortunately, uniform and reliable HKGs for skin research have not been identified in goat. Therefore, this study seeks to identify a set of stable HKGs for the skin tissue of C. hircus using high-throughput sequencing technology. Results Based on the transcriptome dataset of 39 goat skin tissue samples, 8 genes (SRP68, NCBP3, RRAGA, EIF4H, CTBP2, PTPRA, CNBP, and EEF2) with relatively stable expression levels were identified and selected as new candidate HKGs. Commonly used HKGs including SDHA and YWHAZ from a previous study, and 2 conventional genes (ACTB and GAPDH) were also examined. Four different experimental variables: (1) different development stages, (2) hair follicle cycle stages, (3) breeds, and (4) sampling sites were used for determination and validation. Four algorithms (geNorm, NormFinder, BestKeeper, and ΔCt method) and a comprehensive algorithm (ComprFinder, developed in-house) were used to assess the stability of each HKG. It was shown that NCBP3 + SDHA + PTPRA were more stably expressed than previously used genes in all conditions analysis, and that this combination was effective at normalizing target gene expression. Moreover, a new algorithm for comprehensive analysis, ComprFinder, was developed and released. Conclusion This study presents the first list of candidate HKGs for C. hircus skin tissues based on an RNA-seq dataset. We propose that the NCBP3 + SDHA + PTPRA combination could be regarded as a triplet set of HKGs in skin molecular biology experiments in C. hircus and other closely related species. In addition, we also encourage researchers who perform candidate HKG evaluations and who require comprehensive analysis to adopt our new algorithm, ComprFinder.


Background
In molecular biology research, determining the relative changes in target gene expression at the transcriptional level requires precise quantitative analysis. The emergence and development of quantitative real-time polymerase chain reaction (qRT-PCR) has enabled comprehensive mRNA quantification. Furthermore, qRT-PCR is a commonly used technique due to its accuracy, sensitivity, reproducibility, and cost-effectiveness in analyzing gene expression [1,2]. The copy number of nucleic acid was calculated through the changes in real-time fluorescence reaction. The changes is typically reported as a cycle threshold value (Ct) in the comparative Ct method [3]. The qRT-PCR assay relies on house-keeping genes (HKGs) to obtain relative gene expression data [4,5], thus choosing HKGs has become a major source of error and bottlenecks in qRT-PCR experiments.
In qRT-PCR experiments, inadequate HKG selection may lead to an inappropriate interpretation of target gene expression [6]. There are two common mistakes when selecting HKGs: (I) HKGs are selected based on experience without reviewing HKG research study, and (II) a single HKG with poor stability is used. In recent years, it has been reported with increasing frequency that the commonly used HKGs, such as ACTB, GAPDH, and 18sRNA, have critical limitations [7,8]. Ideal endogenous HKGs should exhibit consistent expression levels across all experimental conditions (e.g. cell types, physiological states, and growth conditions) [9,10]. Unfortunately, no HKGs are stable across all experimental conditions, which means that each experimental system may need to use unique HKG(s) to accurately explore the specific research question being investigated.
Goat (Capra hircus) is an economically important livestock species as a source of meat, hair, and dairy products [11]. Skin tissue, as the largest biological organ with important functions including physical protection from injury and infection, thermal insulation, and providing the substrate for growing hair. To reveal the molecular regulatory mechanism of hair follicle activity, it is necessary to clarify the pattern of target gene expression under different conditions, such as different stages of the hair follicle cycle. Unfortunately, most molecular studies examined goat skin have only included a single HKG such as ACTB [12][13][14] or GAPDH [15,16]. In 2014, Bai et al. [17] selected 10 commonly used HKGs based on a literature review to explore their stability in different hair follicle cycles of Liaoning cashmere goats. However, due to the limited number of animals used and testing only of commonly used HKGs, the previously published study [17] resulted in a limited impact. The development of high-throughput RNA-seq technology provides a method of determining spatiotemporal expression at the transcriptome level, and provides a novel approach for the identification of HKGs [18,19]. This strategy was successfully used to identify candidate HKGs for Artemisia sphaerocephala [7], Pyropia yezoensis [20], Euscaphis [21], Arabidopsis pumila [22], fish [23], tomato leaves [24], and holstein cows [25]. Therefore, we hypothesized that the novel, credible HKGs which serve goat skin research can be predicted and validated via transcriptome sequencing data.
In this study, the transcriptome dataset of 39 goat skin tissue samples was analyzed. Potential HKGs were predicted, of which 8 genes (SRP68, NCBP3, RRAGA, EIF4H, CTBP2, PTPRA, CNBP, and EEF2) were selected based on their relatively stable expression levels. Four commonly used HKGs (SDHA, YWHAZ, ACTB, and GAPDH) were selected for comparison. These 12 genes were amplified using qRT-PCR in four groups with different experimental treatments. Four different algorithms (geNorm [26], ΔCt method [27], NormFinder [28], and BestKeeper [29]) and a comprehensive method (ComprFinder, a newly developed method by our team) were used to evaluate the stability of each HKG. Finally, the reliability of the recommended optimal HKGs was validated and confirmed.

Selection of novel candidate HKGs based on RNA-seq data
From a complete transcriptome dataset, the fragments per kilobase of exon model per million mapped reads (FPKM) of all transcripts from each sample were obtained. We first removed some transcripts which did not have a credible function annotation, or exhibited low levels of expression (FPKM = 0). This resulted in 15,853 unigenes being found for further selection. Next, genes with a relatively high expression level (FPKM ≥10 or ≥ the 80th percentile) as determined by the mean FPKM value, and genes with low variability as determined by the coefficient of variation (CV, %), maximum fold change (MFC), and dispersion measure (DPM), were considered (see Methods section). As shown in Fig. 1, the probability density curve of all 15,853 unigenes was evaluated by these 4 indicators.
A Venn diagram was constructed for the 4-color blocks (green, red, yellow, and blue corresponding to those used in Fig. 1a-d, respectively). This showed that 1325 genes (Fig. 1e) met all 4 of the above requirements, and are significantly enriched in 11 signaling pathways (q < 0.05) as shown in Additional file 1: Figure S1. These   HKGs and of these,  8 genes (RRAGA, PTPRA, SRP68, EIF4H, NCBP3,  CTBP2, CNBP, and EEF2) that with lower CV value, higher FPKM value, and easier primers design were selected for further qualification. Besides, 4 genes outside of the initial 1325 were considered, including SDHA and YWHAZ as they had previously been proposed by other researchers [17], and ACTB and GAPDH genes were included as the most commonly used endogenous HKGs for exploring target gene expression in goats. In total, 12 candidate HKGs were analyzed in subsequent steps. Each gene was ranked based on its CV value with a lower CV receiving a higher-ranking order ( Table 1).
Amplification specificity and efficiency of the candidate HKGs and target genes A total of 15 primer pairs including 12 candidate HKGs and 3 target genes were designed for qRT-PCR experiments. Detailed information on gene symbol, primer sequence, and amplicon specifications are shown in Additional file 1: Table S1. Amplification efficiency for all 15 genes ranged from 96.4% for DKK1 to 103.9% for PTPRA, and the coefficient of determination (R 2 ) varied from 0.9986 to 0.9999. The specificity for each paired primer was validated by the melting curve analysis, which showed a single amplification peak (Additional file 1: Figure S2). Each pair of primers had good specificity and amplification efficiency around 100%.

Expression profiles of the candidate HKGs
The mean Ct (the average of 3 technical replicates in the same sample) values were used to calculate gene expression levels among samples with distinct experimental factors. As shown in Fig. 2 and Additional file 1: . The narrower standard deviation range of a gene means it has higher expression stability in different samples. Although some genes had a lower standard deviation than others, experimental errors are always possible. Therefore, to obtain a reliable evaluation of these candidate HKGs, further analysis with more scientific algorithms is needed.

Analysis of HKG expression stability
In this study, 4 publically available algorithms were used to evaluate HKGs for higher-accuracy stability rankings: geNorm, NormFinder, BestKeeper, and the ΔCt method.

geNorm analysis
Gene expression stability was determined by the Mvalue 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 (V n /V n + 1 ). This method recognizes V n /V n + 1 < 0.15 as a threshold value, and "n" as an appropriate number of HKG needed. The V 2 /V 3 values for all 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, 2, 3, 4, and all samples, respectively), which indicate that using double HKGs (first two genes in each group) is sufficiently accurate for use in normalizing qRT-PCR derived gene expression data (Additional file 1: Figure S3). The triplet or more gene combinations can also be used as V n /V n + 1 < 0.15 (n ≥ 3).

NormFinder analysis
Expression stability values, as determined by NormFinder, are shown in Table 2. For group 1, SDHA and EIF4H were the most stable HKGs, and ACTB was the least stable gene, which was the same as was determined by geNorm. In group 2, SDHA and NCBP3 were the most stable HKGs while ACTB was the least stable gene. In group 3, SDHA and YWHAZ got the top rank, while ACTB ranked at the lowest. In group 4, PTPRA and NCBP3 were the most stable, while GAPDH ranked at the lowest. In all samples, SDHA and NCBP3 were the most stable, while ACTB was the least.

BestKeeper analysis
The BestKeeper algorithm used std-values to assess HKG stability with the lower the std-value, the more stable HKG expression was. As shown in Table 3, in group 1, SDHA and PTPRA were the most stable HKGs, whereas RRAGA was the least stable. The same was observed with the geNorm analysis. In group 2, SDHA and EEF2 were the most stable HKGs, while ACTB was the least stable. In group 3, SDHA and YWHAZ got the top rank, while SRP68 ranked at the lowest. In group 4, EEF2 and NCBP3 were most stable, while GAPDH was the least. In all samples, SDHA and EEF2 were most stable, while ACTB was the least.

Δ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 stdvalue, 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.
A comprehensive ranking of the four methods examined Next, the ComprFinder algorithm was employed to obtain a comprehensive score that was used to rank the potential HKGs (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. In group 4, NCBP3, PTPRA, and EEF2 were the most stable genes. The overall rankings, from the highest to the lowest sta- It is interesting to note that the top 3 genes in 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. NCBP3, SDHA, PTPRA were the 3 most stable HKGs across all samples with scores within a tight range, calculated final score (FS) of 0.096, 0.099, and 0.108, respectively. They were also preferably ranked in groups 1-4 relative to other genes and were therefore considered to be the 3 most promising candidate HKGs, and were advanced for further validation.
Validation of the recommended HKGs by DKK1, SHH, and FGF5 genes Based on the above analyses, 3 target genes (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 using 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 their ComprFinder FS values. 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 added to this 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 normalization factors.
As is shown in Fig. 4a, the expression profiles of DKK1 were similarly obtained using the 8 stable singlegene and multi-gene combinations. Furthermore, it was observed that DKK1 was more highly expressed in T2 compared to T1, and it was most highly expressed during the T3. Among the unstable single-and multi-gene combinations, only ACTB and ACTB + GAPDH performed similarly to the stable genes. However, the gene expression profile as normalized by GAPDH was different from the other conditions, and no significant     difference has been identified among T1, T2, and T3. Expression of the SHH gene was even during the T1 and T2, but there was a significant decrease in T3 (Fig. 4b).
The 5 multi-gene combinations and NCBP3, SDHA identified this trend, but PTPRA did not. Though the GAPDH-normalized gene expression profile had similar trends to stable multi-gene combinations, ACTB was different. 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 expression levels were observed in T2, but no statistical significance 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, significant differences in ACTB were also identified in T2 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 multiplegene combinations. 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 shown in Fig. 5, the normalized results using 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 doublegene 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

Discussion
Standard criteria for HKG screening for skin tissue research in goats Which candidate HKGs should we choose?
Four original algorithms were used to identify the expression stability values of 12 candidate HKGs and their FS values were determined using a comprehensive algorithm. However, even for the final ComprFinder value, the results varied between different groups. If the top 3 genes were considered, groups 1-4 should theoretically be EIF4H + PTPRA + SDHA, SDHA + NCBP3 + EEF2, SDHA + PTPRA + EIF4H, and NCBP3 + PTPRA + EEF2, and a total of 5 HKGs (EIF4H, PTPRA, SDHA, NCBP3, EEF2) would be needed in goat skin research. In theory, it is preferable to use multiple high-performing HKGs as a normalization factor. However, in practice, the additional cost and excessive number of HKGs, limit the number of samples that can be tested. Therefore, the minimum number of HKGs should be used to meet the relevant statistical needs, in addition to reducing experimental costs [10,30]. In this study, NCBP3, SDHA, and PTPRA were the top 3 most stable HKGs for all samples, and 2 of those ranked in the top 3 of each group. They were therefore considered as the most three stably expressed HKGs, and for further validation.
Considering that the 3-gene combination of NCBP3 + SDHA + PTPRA exhibited better stability, it can be applied to various factors in goat dermatologic research, and 3 HKGs is still an acceptable number for qRT-PCR experiments. Therefore, it is recommended that NCBP3 + SDHA + PTPRA be adopted as the HKG combination for skin research in goats.

The HKGs of skin tissue in goats and other species
As noted previously, common HKGs used in traditional skin research of goats were either ACTB [12][13][14] or GAPDH [15,16], so erroneous data might be obtained. Therefore, the advantage of a set of appropriate HKGs is very valuable. Regarding target genes that have undergone significant changes, these can be identified by less stable HKGs. But for target genes that show slight changes, these can only be identified by optimal HKGs [30]. As far as we know, only one previous HKG study reported on goat skin tissue, Bai et al. [17] selected ten commonly used HKGs by consulting the literature. The selected HKGs were tested on 3 stages of hair follicle cycle in Liaoning cashmere goats (referred to here as T1, T2, and T3 of IMCG), and authors finally recommended the SDHA + YWHAZ + UBC as the HKGs. But their geNorm values (V 2 /V 3 = 0.159 and V 3 /V 4 = 0.144) imply that the combination of 3 genes was not ideal. In the present study, large number of biological samples were provided for determination and validation, and multiple algorithms were used for evaluation, with the RNA-seq dataset was used for prediction and selection. Therefore, in terms of both the number and quality of HKGs, this study is a significant stepforward from previous studies.
When studying the target gene expression level in skin tissue from other species, such as Angora rabbits [31], mink [32], mice [33], and humans [34], ACTB or GAPD H are generally used as the HKGs. The selection of HKGs from the same type of tissue within neighboring species has been widely recognized and accepted. Thus, the data presented here could prove HKGs are suited to skin tissue research, not only for goats, but also for other species.

Selection and validation of HKGs based on RNA-seq data
Selection and validation of HKGs using RNA-seq produced more reproducible results, had greater sensitivity, yielded better correlation with protein expression levels, and had more accurate detection and higher coverage [35]. To the best of our knowledge, this study is the first one to report on the selection and validation of novel HKGs for qRT-PCR analysis in goats. Two novel HKGs (NCBP3 and PTPRA) and a known HKG (SDHA) belonging to the NCBP3 + SDHA + PTPRA combination were recommended. Using a similar approach as other studies [19][20][21], new and improved HKGs were identified through analyzing an RNAseq dataset. While this study demonstrated the advantages of using RNA-seq datasets in the discovery of new HKGs, it is also possible that the prediction of HKGs by RNA-seq datasets may be lacking in some respects. For example, the ranking order of these candidate HKGs (Table 1) and the determined FS (Table 5) did not match (compared in Additional file 1: Table S3). Specifically, the CV values of RRAGA and SRP68 were in the top 3, although in the final ranking they did not appear in the top 3 positions of any group (groups 1-4 and all samples). This might be the reason why RNA-seq samples in the selection stage and the qRT-PCR samples used in the determination stage did not completely overlap. This phenomenon is consistent with those reported by Gao et al. [20] and implies that the HKGs predicted by RNA-seq screening were not fully reliable, and need further validation by qRT-PCR experiments.
As mentioned before, ACTB and GAPDH are currently the most popular HKG in the literature, but their limited normalizing capacity was verified here. This suggests that scientists must be cautious when selecting traditional HKGs, especially when identifying target genes that have slight changes in expression. Therefore, it is recommended that common HKGs be included as a comparison, to provide direct evidence. Of course, it must also be acknowledged that mining reliable HKGs require scientific experimental design, complete experimental materials, more algorithmic tools, and a certain amount of scientific research time [30], but these are not available to every research laboratory. Therefore, those experimental systems that do not meet the above conditions were recommended to search for HKGs in close species using the ICG platform [36]. Meantime, use as many HKGs as possible and calculate their arithmetic mean as the normalization factor to increase experimental stability, instead of simply using a single HKG such as ACTB or GAPDH.

Comprehensive analysis
The requirement of comprehensive analysis, and the shortcomings of the previous algorithms After evaluating candidate HKGs with the abovementioned algorithms (geNorm [26], NormFinder [28], BestKeeper [29], and the ΔCt method [27]), it is not surprising that the rankings of candidate genes may vary depending on the algorithm used [30]. Thus, another algorithm is needed for comprehensive ranking. After reviewing the literature on HKGs, it was determined that 2 types of comprehensive algorithms were mainly used: (1) the primary ranking order was used to calculate the arithmetic average, and then get the final ranking [24,37]; (2) using the primary ranking order, the geometric average is then calculated to get the final ranking. RefFinder [38] is a typical representative (https://www.heartcure.com.au/ reffinder/?type=reference) and many studies [7,19,20,22,[39][40][41] have used it (Times Cited: 352, on Web of Science, 2020/2/3), which illustrates its vast impact. Both the above comprehensive algorithm types depend on the ranking number of the original algorithms. Due to this, it is possible that the use of ranking numbers may cause some errors.
These ranking numbers reflect the true size of stability values of these candidate HKGs, but they should not be used as the input numerical value for the next calculation. Doing so would excessively reduce or enlarge the real differences among them. Referring to our experimental data as an example, these candidate HKGs were evaluated and unevenly distributed (Fig. 6a) on the axis of the four algorithms. The RefFinder algorithm provides them with the uniform rank of 1-12 (Fig. 6b), and then calculates the geometric average for each candidate HKG. The RefFinder algorithm increases the gap between PTPRA and EIF4H, and reduces the gap between CNBP and ACTB (the axis of the geNorm algorithm in Fig. 6a and b). We consider this was unreasonable for the comprehensive evaluation and therefore, in this study, a new algorithm for comprehensive analysis was developed.

The ComprFinder algorithm
ComprFinder was intentionally developed to replace RefFinder and the comparison of the two algorithm's results are shown in Additional file 1: Table S4. It can be seen that the FS and ranking order of the 12 candidate HKGs calculated using the two algorithms are different. For example, NCBP3 (0.096), SDHA (0.099), and PTPRA (0.108) were the top 3 genes for all samples calculated by ComprFinder, whereas PTPRA (2.14), SDHA (2.45), and EEF2 (3.13) were the top three calculated by RefFinder. The ComprFinder algorithm directly standardizes the results of the original algorithms, unlike RefFinder, which uses the ranking numbers of original algorithms. Theoretically, if the series of HKG expression stability values calculated by each original algorithm is uniformly distributed (likely does not exist), the results of Compr-Finder and RefFinder algorithms will be the same. Using the ranking numbers, the real differences among the candidate HKGs will excessively reduce or increase. Therefore, the new comprehensive algorithm can overcome the intrinsic errors caused by artificial assignment. Standardized processing allows the results of different algorithms to have the same dimension and makes them essentially comparable. Finally, the standardized results can be integrated to get a series of scores and a final evaluation. Therefore, ComprFinder would be a more reasonable algorithm than RefFinder for the comprehensive evaluation of HKGs.
We provide a ComprFinder algorithm tool (Additional file 2) for researchers who have comprehensive evaluation needs for candidate HKGs. This tool is based in Microsoft Excel and can be downloaded from the supplementary materials of this article. Briefly, after inputting the original algorithm results into the input area, the ComprFinder algorithm automatically processes the data and all candidate HKGs will be scored and presented in the output area. Although the use of ComprFinder in this study was based on geNorm, BestKeeper, NormFinder, and the ΔCt method, analysis is not limited to these 4 algorithms.

Conclusion
In this study, we present the first list of candidate HKG selection for goat skin tissue based on transcriptome data. The NCBP3 + SDHA + PTPRA combination was identified and recommended as the triplet HKGs for skin molecular biology studies in goats and other closely related species. In addition, a comprehensive algorithm tool based in Microsoft Excel, ComprFinder, was developed for the comprehensive evaluation of candidate HKGs.

Animals and skin tissue samples
All animals and sampling procedures in this study were supervised and approved by the Institutional Animal Care and Use Committee of Southwest University. Each 1 cm 2 skin tissue was sampled from the Inner Mongolia cashmere goat (IMCG), Dazu black goat (DBG), Hechuan white goat (HCWG), or the first filial generation (F 1 , DBG♂ × IMCG♀). Detailed information regarding the animal source, method of the animal sacrificed, anesthesia procedure, sampling procedure, and sample preservation is found in Additional file 1: Tables S5 and S6. All samples were stored at − 80°C until further usage.
A total of 48 skin tissue samples were collected to determine (determination stage) and to validate (validation stage) the potential HKGs. The sample sizes (n) for every group (factor) and every level were described in the legend of Fig. 7. In the determination stage of this study (Fig. 7c, Additional file 1: Table S5), all samples were collected from does. Four groups were used, including age (4 development stages, group 1), sampling time (3-stages of hair follicle cycle, group 2), breed (4 different breeds, group 3), and sampling site (5 different sampling sites on the body of the goat, group 4). As shown in Fig. 7c, group 1 included F 1 _P0, F 1 _P60, F 1 _P240, and F 1 _Adult, which were sampled at 0-day, 2 months, 8 months, and 2 years after birth from F 1 . Group 2 included IMCG_T1, IMCG_T2, and IMCG_T3, which were sampled during the anagen (September), catagen (December), and telogen (March) from IMCG. Group 3 consisted of animals from 4 breeds and included IMCG, DBG, F 1 , and HCWG. Skin samples from each of the breeds were sampled in the anagen phase of the hair follicle cycle. Skin biopsies collected from group 4 (#4, #5, #6, #12, and #14 from IMCG, as described in our previous publication [42]), were taken from the forearm, dorsal chest, lateral chest, thigh, and the inner side of the forearm as described in our previous study. Except for samples #4, #5, #12, and #14, all samples analyzed here were collected from the lateral chest of the goat body. Except for samples belonging to F 1 _P0, F 1 _P60, and F 1 _ P240, all samples were collected from adult goats. In the validation stage of this study (Fig. 7e, Additional file 1: Table S6), 3 bucks of IMCG were added to group 2, to enhance validation accuracy.

RNA isolation and cDNA synthesis
Total RNA was extracted using the RNAiso Plus kit (#9109, TaKaRa, China) according to the manufacturer's Fig. 6 Gene stability values and rank order. The same color presents the same gene. a Twelve candidate HKGs unevenly distributed on the axis, this is their true distribution; b Twelve candidate HKGs evenly distributed on the axis ordered by [1][2][3][4][5][6][7][8][9][10][11][12] instructions. The concentration and purity were determined using the Nanodrop2000 (Thermo, USA) with the 260/280 ratios being between 1.8 and 2.0, and the 260/230 ratios were greater than 1.6 in all analyzed RNA samples. First-strand cDNA was synthesized using the 5X All-In-One RT MasterMix (with AccuRT Genomic DNA Removal Kit) (#G492, ABM, Canada) and 1:4 volumes of DEPC water was added to dilute the samples.

Selection of candidate HKGs
Transcriptome sequencing data of 39 goat skin tissues (unpublished data, Fig. 7a) was performed using the paired-end sequencing technology on an Illumina Nova-Seq 6000 platform. After assembly and annotation, the gene expression profiles and read counts of unique transcripts were converted into FPKM values on the platform BMKCloud (www.biocloud.net), according to the formula FPKM = cDNA fragments/[mapped fragments (millions) × transcript length (kb)]. Based on the FPKM value of every gene in each transcriptome, the CV and MFC were calculated using Microsoft Excel and DPM was calculated using a jar package [43]. The CV is defined as the ratio of the SD FPKM to the mean of the FPKM of all samples for one gene. The DPM parameter was introduced for the identification of the HKGs on pattern gene finder (PaGeFinder) [43], and their jar package was downloaded and run from the PaGeFinder website. The MFC, which is defined as the fold change between the largest and smallest FPKM values within 39 RNA-seq were calculated. The standard criteria for HKGs are relatively high expression level, and low expression variation, therefore, a candidate HKG should have a relatively high FPKM value, and low CV, DPM and MFC values.
Genes with RPKM, CV, MFC, and DPM fulfilling the criteria of HKGs were retained for further analysis (Fig. 7b). Moreover, two HKGs (SDHA and YWHAZ) from a previous study by Bai [17], and two commonly used HKGs (ACTB and GAPDH) were also considered as HKGs. All candidate HKGs were amplified using qRT-PCR for subsequent determination and validation. The probability density curve was drawn by an in house script (Additional file 3) using the Matlab software (https://ww2.mathworks.cn/ products/matlab.html). Venn diagram analysis was performed using the OmicShare online platform tools (http:// www.omicshare.com/tools).

Primer design and amplification efficiency analysis
Specific primers were designed using the Primer-BLAST [44] web tool (https://www.ncbi.nlm.nih.gov/tools/ Fig. 7 The workflow of this study. a The sample information of transcriptome sequence data of these IMCG (n = 3*3), DBG (n = 3*3) and F1_Adult (n = 7*3); b Candidate housekeeping genes were preliminarily selected by four indicators which including FPKM, CV, DPM, and MFC, and were further selected by Venn diagram analysis; c The sample information of the qRT-PCR experiments on the 4 experimental groups, with 3 biological replicates in every level of each group. Group 1, different development stages (n = 4*3); Group 2, hair follicle cycle stages (n = 3*3); Group 3, breeds (n = 3*3); Group 4, sampling sites (n = 5*3). d Candidate housekeeping genes were determined using 4 algorithms, including geNorm, NormFinder, BestKeeper and the ΔCt method. An additional comprehensive analysis was conducted using ComprFinder, a new algorithm developed by the authors. e The selected HKGs were validated by 3 target genes, and were performed on 18 skin samples (n = 6*3) primer-blast/) based on the sequences of the unigenes. The criteria for primer design were as follows: primer lengths of 17-24 bp, GC content of 50-66%, theoretical anneal at around 60°C, and amplicon lengths of 100-200 bp. All primers were synthesized by the Beijing Genomics Institute (Beijing, China).

qRT-PCR analysis
Sample reactions were performed in a 10 μL reaction volume with 5 μL of 2× qPCR MasterMix (#MasterMix-S, ABM, Canada), 1 μL cDNA template, 0.3 μL each primer, and 3.4 μL DNase/RNase-free water and run on the Bio-Rad CFX96 Real-Time PCR Detection System. The thermal cycling conditions were conducted according to the reagent kit instructions as follows: enzyme activation at 95°C for 10 mins, followed by 40 cycles of denaturation at 95°C for 15 s, annealing/extension at 60°C for 60 s. The specificity of the SYBR green PCR signal was confirmed by melting curve analysis. All samples were analyzed in 3 replicates. Serial tenfold dilutions (dilution ratio 1:10 3 -1:10 10 ) of cDNA template (note: in primer amplification test, cDNA template is PCR product with the same primer) were used to generate a slope of the standard curve to calculate the amplification efficiency and R 2 of each paired primer. All qRT-PCR experiments and data analyses in the present study were performed following the MIQE guidelines [45,46].
A comprehensive analysis was performed using a newly developed algorithm After using the above-mentioned traditional evaluation algorithms, another algorithm was developed for comprehensive ranking. The new algorithm, ComprFinder, standardizes the output values from the above 4 algorithms then arithmetically averages them to get their FS and final ranking order (Fig. 7d). Fig. 8 Schematic diagram of the ComprFinder algorithm. The same color presents the same gene. From step 1 to step 2, the uneven distribution of these 12 HKGs will not change, but they are proportionally enlarged or reduced to the range of 0-1. Calculated the arithmetic mean value for each gene to determine the final ranking order The schematic diagram of the ComprFinder algorithm is presented in Fig. 8. Specifically, STEP 1, according to the original algorithm results ordered the values from small to large to find the minimum and the maximum. From this, the range was calculated (1) where A is one of the original algorithms, and is also used to calculate Δ B , Δ C , and Δ D . In STEP 2, the values from the original algorithm were standardized in the interval [0,1], where the minimum value = 0 and the maximum value = 1. The normalized value (A' i ) was calculated (2) where i is one of these candidate HKGs, and is also used to calculate B' i , C' i , and, D' i . All the other data were assigned normalized values between 0 and 1. This step makes these stability values that belong to the same HKG but from different original algorithms abide by the additive property. In STEP 3, arithmetic averages of the standardized values for each candidate HKG were calculated, and then these FS were sorted to obtain their final rankings (Fig. 8). The FS i (3) was determined with n being the number of original algorithms we used.
Experimental validation of the HKGs To verify the results, the 3 best candidate HKGs were selected (NCBP3, SDHA, and PTPRA), as well as the most unstable and most commonly used HKGs (ACTB and GAPDH). Next, the HKGs were verified and evaluated using 3 target genes (DKK1, SHH, and FGF5), which are the hot genes in hair follicle research (Fig. 7e).
Considering the need to accurately evaluate target gene expression profile, 3 new bucks were added in the original sample size of 3 does, (6 adult IMCG, 3♂ and 3♀) and were sampled over 3 time-points (T1, T2, and T3) (Additional file 1: Table S6). The qRT-PCR was conducted as described in the determination stage. The paired sample t-test was performed using Microsoft Excel, and the graph was plotted using GraphPad Prism 6. The results are presented as Mean ± SEM, * P < 0.05, ** P < 0.01. For multiple gene combinations, the geometric average of their Ct value was calculated [26]. The relative gene expression level was calculated as 2 -ΔCt , ΔCt = Δ (Ct target gene -Ct HKGs ).
To further evaluate the internal relationship of these candidate HKGs, a correlation analysis was performed. First, the target genes were normalized by different HKGs or HKG combinations, to calculate their ΔCt values. Then their copy numbers were converted to relative expression levels with 2 -ΔCt . Finally, their normalized-based expression levels were examined by the correlation analysis.