Validation of Housekeeping Genes to Study Human Gingival Stem Cells and Their In Vitro Osteogenic Differentiation Using Real-Time RT-qPCR

Gingival stem cells (GSCs) are recently isolated multipotent cells. Their osteogenic capacity has been validated in vitro and may be transferred to human cell therapy for maxillary large bone defects, as they share a neural crest cell origin with jaw bone cells. RT-qPCR is a widely used technique to study gene expression and may help us to follow osteoblast differentiation of GSCs. For accurate results, the choice of reliable housekeeping genes (HKGs) is crucial. The aim of this study was to select the most reliable HKGs for GSCs study and their osteogenic differentiation (dGSCs). The analysis was performed with ten selected HKGs using four algorithms: ΔCt comparative method, GeNorm, BestKeeper, and NormFinder. This study demonstrated that three HKGs, SDHA, ACTB, and B2M, were the most stable to study GSC, whereas TBP, SDHA, and ALAS1 were the most reliable to study dGSCs. The comparison to stem cells of mesenchymal origin (ASCs) showed that SDHA/HPRT1 were the most appropriate for ASCs study. The choice of suitable HKGs for GSCs is important as it gave access to an accurate analysis of osteogenic differentiation. It will allow further study of this interesting stem cells source for future human therapy.


Introduction
Gingival stem cells (GSCs) are emerging new tools that have been recently isolated from human gingival tissue [1], a tissue with high healing capacity [2]. These cells can differentiate into osteoblasts, chondrocytes, and adipocytes when appropriately induced [1,3]. An easy, noninvasive, and scar-free access to gingival tissue makes GSCs interesting cells and a potential source for future regenerative medicine applications, in order to replace or repair diseased tissues and organs [2].
The osteogenic differentiation of GSCs has been evaluated in vitro and confirmed using qualitative approaches such as Alizarin Red S staining [1]. This potential may be used to treat large maxillary bone defects. GSCs may be an alternative to stem cells derived from mesoderm, because they share their embryonic origin with facial bones, that is, cephalic neural crests [4,5]. In order to better characterize these cells and their in vitro osteogenic differentiation, molecular biology quantitative assays can be of high interest, using gene expression study. This approach is more accurate than histochemical qualitative methods to study GSCs and confirm their stemness and differentiation potencies. To achieve this purpose, reverse transcription quantitative Polymerase Chain Reaction (RT-qPCR) technique is considered the most accurate and reliable method, owing to its sensitivity, realtime detection of reaction progress, and the amplification of very small quantities of mRNA [6].
However, reproducibility problems and biases have been recently discussed [7]. Indeed, many factors may affect the RT-qPCR accuracy. These are related to operator's manipulation: pipetting, treating samples, cell culture, the amount 2 Stem Cells International of starting material, different conditions of harvesting cell samples, RNA extraction yield, mRNA quality, and homogeneity taken from tissues or cells. DNase treatment, reverse transcription enzyme specificity, the efficiency of reverse transcription reaction, the type of fluorescence, and the temperatures of qPCR cycles may also provoke differences [6].
In order to minimize errors related to these factors, data normalization is mandatory. Other solutions have been recently proposed, such as using robots, RNA extraction kits, or normalized kits from cell culture to cDNA synthesis, which showed high sensitivity, accuracy, and reproducibility [8]. However, these techniques remain expensive and not accessible to all research centers. The standardized normalization of qPCR results against reference genes thus remains the most common method to have reproducible and reliable results [6].
Reference genes or housekeeping genes (HKGs) are internal reaction control [9]. It is considered that these HKGs reflect basic metabolism implicated in the general processes and mechanisms of cell cycle and hence are supposed not to be affected or regulated by experimental conditions. Several criteria are required for a reliable reference gene: the expression level must be unchanged by experimental conditions, and its variability must be as minimal as possible in the different samples. It would be preferable that the cycle threshold (Ct) of such a HKG is close to the one of the gene of interest. It is often recommended to use at least two HKGs or more, because only one may lead to biases or may be regulated by experimental conditions without the possibility to control it [6]. Thus, HKG choice is a crucial step of the RT-qPCR protocols, and the investigator must carefully validate it. However, such information remains often lacking [10].
To avoid mischoices, several statistical approaches and algorithms have been proposed [11][12][13][14]. These methods have been employed on different mesenchymal stem cells, but to our knowledge not yet on oral stem cells [15]. The objective of this study is to choose the most reliable HKGs to study GSCs and their osteoblast differentiation by RT-qPCR. To validate these HKGs, we will assess the 10 most used HKGs with four published algorithms: the ΔCt comparative method, GeNorm, BestKeeper, and NormFinder.

Materials and Methods
The study has been carried out thanks to the collaboration between the Laboratory of Molecular Oral Physiopathology, INSERM UMRS 1138, Cordeliers Research Center, Paris, France, and the Hospital Complex Henri-Mondor Albert-Chenevier, CIC-BT-504, Creteil, France.
Written informed consents were obtained from patients in accordance with Helsinki Declaration Principles, taking into account the ethical, legal, and regulatory norms and standards for research in France (Loi Huriet number 91-73), as well as the applicable international norms and standards.
Our research protocol was submitted to the research ethics committee: Centre d'Investigations Cliniques de Creteil-Henri Mondor (Biothérapie) CIC-BT-504, directed by Pr. José COHEN, which specifically approved this study. (GSCs). GSC cultures were obtained from six healthy donors ( = 6), after their informed consent. Gingival tissues were collected from gingiva of tooth extraction sites from the donors during surgery. Stem cells were isolated either by explants technique, followed by fibroblast colony forming units (CFU-F) assay, or by enzymatic digestion with Collagenase II (Sigma-Aldrich) [16][17][18]. Eight GSCs samples were finally obtained (six samples from enzymatic digestion technique and two samples from explants technique). They were seeded in 35 mm Petri dishes, proliferated with 10 mL of Dulbecco's Modified Eagle Medium Low Glucose (DMEM-LG), GlutaMAX, and pyruvate Supplement (Gibco-Life Technologies, Carlsbad, CA, USA). The medium was supplemented with 10% of foetal calf serum (FCS; Gibco), 10 mL/L Penicillin-Streptomycin (5 UI/mL, Gibco), 1% of nonessential amino acids (NEAA; Gibco), and 2.5 mg/L of Amphotericin B (250 g/mL; Gibco). Culture medium was filtered in 0.22 m pore size filters before use. Cell cultures were maintained at 37 ∘ C in a humidified 5% CO 2 incubator. Culture media were changed twice a week, with supplementation of L-Ascorbic Acid 2-Phosphate (50 g/mL; Sigma-Aldrich), until they attempted 70 to 80% of confluence. Cells were then harvested after Trypsin-EDTA treatment and centrifugation. They were frozen at −80 ∘ C and stored until use for RNA extraction. GSCs were characterized using flow cytometry with stem cells markers: CD29, CD73, CD90, CD105, CD45, and HLA-DR [1].

Adipose Derived Stem Cells (ASCs; Purchased from ZenBio
Company, USA). Six samples of ASCs ( = 6), mesenchymal stem cells of adipose tissue origin, used as control, were also cultured in the same conditions, in early passages (1)(2)(3)(4)(5), in order to compare them to GSC. When reaching 70-80% of confluence, they were collected and frozen at −80 ∘ C for RNA extraction.

Histochemical Staining.
To confirm osteoinduction, fixed cells were washed with phosphate-buffered saline (PBS, Gibco), followed by two distilled water rinses and incubation in fresh Alizarin Red S solution (1 g/100 mL in distilled water, Sigma-Aldrich), pH around 4.10, for 30 min. The wells were then rinsed repeatedly with distilled water, and calcium deposits were observed (Figure 1). 2.5. Total RNA Extraction and cDNA Synthesis. Total RNA was isolated (GSC = 8, ASC = 6, and differentiated GSC = 8), using ReliaPrep RNA Cell Miniprep System kit (Promega), according to manufacturer's instructions. This kit fits cell culture and incorporates a recombinant DNase treatment step. RNA concentrations and purity were assessed by a NanoDrop spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The purity was determined by measuring the absorbance at 260 nm/280 nm. The ratio 260 / 280 is expected to be between 1.8 and 2.0. RNA quality was also confirmed by agarose gel electrophoresis, choosing 4 random samples, which confirmed the absence of ribosomal RNA degradation, with a 28S/18S ratio around 2 (Figure 1(f)).

A GSC A A A A A A A A A A A A A A A A
2 g of RNA of each sample was reverse-transcribed using SuperScript II Reverse Transcriptase (Invitrogen) initiated by random primers and oligo-dT primers, as described in the provided protocol. The final volume was 20 L. The latter was 20-fold diluted (v/v) to have the same final concentration of 5 ng/ L of cDNA. Reactions were also prepared without the reverse transcriptase enzyme, which gave minus RT products to confirm the absence of genomic DNA contamination in qPCR reactions.   (Table 1). They were selected with an amplicon size <250 bp, spanning an exon-exon junction when possible, with no polymorphism, and the ideal melting temperature was 60 ∘ C (with a maximum difference of 3 ∘ C between every two primers). For primer nucleotides, they should end with a C or G residue, and CG content was ranged between 50 and 60%. The cycle threshold (Ct) was ranged from 15 to 30. The same tool and settings were used to generate the primers for the analysis of stem cells osteogenic differentiation, such as RUNX2.
2.7. RT q-PCR Assay. The assay was carried out in triplicate in a 96-well format in the Bio-Rad CFX96 Real Time System (Bio-Rad Laboratories, USA). The reaction volume was 15 L. For each sample, the real-time PCR mixtures consisted of 3 L cDNA (≈15 ng of cDNA) and 12 L of mixture of 7.5 L of SYBR Green Supermix (Bio-Rad), 0.45 L of primer (250 nM final of each sense and antisense sequence of the primer), and 4.05 L of RNase/DNase-free sterile water. The assay included negative controls (nuclease-free water) and controls with minus RT to detect reagent contamination and the presence of genomic DNA, respectively. In order to define the efficiency ( ) of each primer, serial dilutions (five dilutions: from half to five times alternately) of one GSC sample which expresses the genes of interest and the housekeeping genes and relative standard curves were generated. Efficiency was calculated as follows: PCR efficiency = (10 [−1/slope] − 1) × 100.
The reactions were run for 35 cycles following this thermal cycling profile: (1) 95 ∘ C for 3 min, (2) denaturation at 95 ∘ C for 5 sec, (3) primer annealing step at 60 ∘ C for 20 sec (the most optimal temperature for all primer pairs after performing preliminary real-time RT q-PCR assays).
To confirm primer specificity, a melting curve analysis was performed after each amplification, ranging from 65 ∘ C to 95 ∘ C, with temperature increasing steps of 0.5 ∘ C every 5 sec. In all negative samples, no fluorescent signal was detected, or at very late cycle (more than 10 cycles after the Ct). This ensured of the quality of RNA extraction in all the samples.

Analysis of Gene Expression.
In order to assess the most stable reference genes to study GSC, ASC, and dGSC, Ct of all samples were calculated (Figures 2(c), 2(d), and 2(e)). Ct is defined as the number of cycles needed for the fluorescence signal to reach a specific threshold of 2.8.1. The ΔCt Comparative Method. This approach is based on the comparison of Ct differences, or ΔCt values of each "pairs of genes" within each sample. The objective is to determine if the reference gene has an increased or decreased level of deviation in ΔCt among the samples when compared to the other reference genes. If the ΔCt remains stable when analyzed in different samples, it means that these two genes are stably expressed all among the samples or that they are coregulated [13]. ΔCt mean is calculated for each pair of genes and between every two groups of HKGs. All genes are taken into account and all possible "gene combinations" are compared. Then, the standard deviation values of each ΔCt set of each pair of genes is determined. Reference genes are ranked according to the arithmetic mean of standard deviation values, which must be as low as possible to be the most stable. Box-and-whiskers charts can also be generated to show the distribution of ΔCt values in each pair of genes in the samples and allow us to compare ΔCt variability of each reference gene against all others. The bar is a line which represents the median (middle of dataset). The 75 and 25 percentiles represent the upper and lower limits of the boxes, and the whiskers refer to the highest and lowest ΔCt values among the samples excluding outsiders. The HKG with the lowest variability is the most stable [21].

GeNorm
Analysis. This method determines, among set of candidate genes, the two most stable ones that share a similar expression profile throughout all studied samples 6 Stem Cells International [22]. For each two candidate genes (gene and gene ), using normalized values ( , here ), the algorithm calculates the array of log 2 transformed expression ratios and the standard deviation of the pairwise variation of this gene toward all other genes, as shown in the following formula: where is array of pairwise variation of normalized values of gene j toward gene k for all samples; is standard deviation of values for gene j toward gene . Then the stability value ( ), defined as the average or arithmetic mean of the standard deviations of this gene, is calculated [14]: The algorithm ranks HKGs toward their values, from the most stable with the lowest to the least stable. If value is more than 1.5, the gene is considered not stable and is removed from the analysis. GeNorm recalculates, hence, values for the remaining stable genes and ranks reference genes from the most stable to the least stable one. Stepwise exclusion of the least stable gene with the highest value will ultimately result in the two most stable genes that cannot be further ranked and selected a pair of genes as the most stable. The other genes are hence ranked based on their highest compatibility with of the first pair [14].
In order to define accurately the number of HKGs needed for the normalization, GeNorm proposes a second step of calculation of the normalization factor, which is the variation between each pair of genes consecutively ranked from the most to the least stable. If the values are smaller than 0.15, there is no need to add a new HKG. An Excel Add-In method is available to make GeNorm analysis and graphs excel, http://medgen.ugent.be/∼jvdesomp/genorm/.

BestKeeper Analysis. Unlike ΔCt comparative method
and GeNorm software, BestKeeper considers not only "intergene" relation, but takes also into account the "intragene" variation, which assesses samples variation in the same gene. In this approach, the ideal HKGs are expected to have a stable expression in the same tissue or sample [12]. For this, standard deviation (SD: ±Ct) and covariance (CV: %Ct) are the two important values to assess HKG stability. For a studied reference gene, SD of Ct values must be as low as possible and must not be >1; otherwise it will be excluded. A second step consists of studying the correlation with BestKeeper index and selecting the most correlated reference gene to this index. The tool simulates an ideal HKG, or the "best" HKG, for which it calculates, for each sample, a BestKeeper index defined as the geometric mean of Ct values of all HKGs for this sample, as shown in the following formula: see [23].
Then it compares this index to each reference gene by a pairwise correlation and regression analyses and calculates a coefficient of correlation ( ) and the probability . " " should be as close to 1 as possible and " " as low as possible. The BestKeeper index can also be compared to a further ten genes, for the same samples, to decide whether they are differentially expressed. Moreover, this algorithm calculates invariances, InVar (±Ct), of each sample to validate its stability; these values must be <3-fold the least one. A freely Excel based spreadsheet software exists and allows us to calculate automatically these values to make the choice of the most stable HKGs easier. It has been established that sample size has a minimal effect on this tool [23], http://www .gene-quantification.de/bestkeeper.html#download.

NormFinder Analysis.
This model takes into account the influence of variation of gene expression in the same sample and coregulation between different individuals [11]. For each candidate gene, the algorithm calculates a stability value , taking into account the intra-and intergroup variances. The reference gene which has the smallest value is ranked as the most stable, as it has the smallest variation over all samples. By creating subgroups, the algorithm can also select the most appropriate candidate genes to study two groups of samples, choosing the ones with the minimal intraand intergroup variations. Intragroup variation must be as small as possible, and intergroup variation ≈ 0. The results are affected by the number of samples and are more accurate when increases ( ≥ 8). It also requires at least 5 to 10 HKGs.
In order to compare the 4 algorithms and select the most stable HKGs, we applied RefFinder (http://www.leonxie .com/referencegene.php), a web-based comprehensive tool, which uses the currently available algorithms, ΔCt comparative methods [13], GeNorm [14], BestKeeper [12], and NormFinder [11], and assigns an appropriate weight to an individual gene and calculates the geometric mean in order to rank all the studied reference genes.

GSC Osteogenic Induction and RNA Quality Control.
Eight GSC samples were obtained from healthy patients during surgical treatment. Six samples of ASCs were purchased and cultured in the same conditions. Both GSCs and ASCs were induced into osteogenic differentiation in duplicate for 21 days. Cellular matrix in proliferation medium was negative to the Alizarin Red S staining (Figures 1(a) and 1(b)). Mineral nodules and matrix were present in osteogenic medium and confirmed with this histochemical staining for both stem cells (Figures 1(c) and 1(d)). This stresses that GSCs are capable of osteogenic differentiation in vitro like ASCs, which Stem Cells International 7 are already known to form minerals in vitro under these conditions of culture. GSCs were CD29+, CD73+, CD90+, CD105+, CD45−, and HLA-DR− by flow cytometry analysis (data not shown) and form CFU-F (Figure 1(e)). RNA of both GSC and ASC was collected, with a yield of more than 1 g. The analyses of purity with the NanoDrop spectrophotometer showed suitable ratios 260 / 280 with values around 2.00. RNA quality was confirmed randomly for 4 samples of GSCs and ASCs by agarose gel electrophoresis, insuring the integrity of ribosomal RNA, with two bands of 18S and 28S (Figure 1(f)).

The ΔCt Comparative Method.
This method is based on comparing the variability of expression levels of each "pair of genes" within all the samples, in order to identify the reference gene with the lowest variability to be ranked as the most stable. The method was applied in each group of samples separately: GSCs, dGSCs, and ASCs. For each pair of genes, ΔCt values were calculated in all the samples, as well as the Ct mean and the standard deviations (StdDevs) of ΔCt values (for GSCs, Table 3; for ASCs and dGSCs, Supplementary Material 2). The arithmetic mean of StdDevs was then identified for each reference gene and HKGs were ranked according to these values. This is shown in boxand-whiskers charts with all possibilities for all the 10 HKGs (Figures 3(a), 3(b), and 3(c)). In GSC samples, SDHA, ACTB, and B2M were the most stable genes (StdDev mean = 0.97, 0.98, and 0.98, resp.). HPRT1 and RPII were the least stable and had a higher variability within the samples (StdDev mean = 1.38 and 1.64, resp.). In dGSCs, TBP, SDHA, and ALAS1 showed the lowest variability in the samples as compared to the other genes, with a StdDev mean of 0.69, 0.72, and 0.78, respectively, whereas GAPDH and RPS18 showed the highest variability. In ASCs, the most stable genes were SDHA, HPRT1, and TBP, with the lowest StdDev means (0.66, 0.71, and 0.77, resp.), while ALAS1 and ACTB were the least stable ones.

GeNorm Analysis.
Normalized expression values of each group: GSCs, dGSCs, and ASCs were entered in the Excel spreadsheet under the format recommended by the authors, respectively [14].
Firstly, the stability value was calculated for all HKGs, and reference genes were ranked from the most stable, with the lowest values, to the least stable ones, as shown in Table 4. Graphs were generated showing the average expression stability values, selecting the most stable pair of genes, according to which it ranks the other genes (Figure 4(a)). For GSCs ( = 8), the most stable genes were SDHA, B2M, and ACTB, ( = 0.892, 0.919, and 0.924, resp.). RPII was the least stable ( = 1.588) and excluded from analysis ( > 1.5). Hence, RPS 18 and ALAS1 were the least stable ( = 1.162, 1.163, resp.). SDHA and ACTB were selected as the most stable pair of genes (Figure 4(a1)). Regarding dGSCs ( = 8), TBP, SDHA, ALAS1, and RPII were the most stable with values, 0.656, 0.677, 0.729, and 0.771, respectively. UBC and RPS18 had the highest values ( = 0.858; 0.925) and were considered the least stable genes. TBP and ALAS1 were thus selected as the best pair of genes (Figure 4(a2)). As for ASCs ( = 6), SDHA, HPRT1, and TBP were the most stable ( = 0.619, 0.656, and 0.717, resp.). ALAS1 and ACTB were the least stable ( = 1,081; 1,491). GAPDH and B2M were the optimal pair of genes according to the software (Figure 4(a3)). Pairwise variation analysis by calculating two sequential normalization factors (NF and NF +1 ) suggested that the optimal number of reference genes to study teach group of stem cells was three HKGs for GSCs and dGSCs ( 3/4 = 0.165 for both of them, which is around 0.15 as suggested by the software) (Figures 4(b1) and 4(b2)) and two reference genes for ASCs ( 2/3 = 0.134 below 0.15) (Figure 4(b3)).

BestKeeper Analysis.
Raw Ct values of each group (GSCs ( = 8), dGSCs ( = 8), and ASCs ( = 6)) were uploaded to BestKeeper software in the form of an Excel spreadsheet. Firstly, BestKeeper excluded candidate HKGs with highest standard deviation (SD) and the covariance (CV). For GSCs, RPII and ALAS1 had the highest SD (1.66, 1.31) and CV (6.99, 5.18), respectively, and were therefore excluded from analysis. Using pairwise correlation analysis and regression analysis, BestKeeper then compared the intergene relations of remaining candidate genes based on its index. UBC, HPRT1, TBP, and RPS18 had a low correlation with BestKeeper index ( = 0.685, −0.050, 0.622, and 0.684, resp.) and thus were also excluded from analysis, even though they may have low SD and CV like UBC, HPRT1, and TBP ( Figure 5 Table 5(c).

NormFinder Analysis.
Normalized values of the three groups were introduced into NormFinder software as recommended by the authors [11]. The algorithm is presented as a free complement in Excel program. The analysis calculates the value of stability of each reference gene and ranks them from the most stable with the lowest to the least stable.
For intragroup analysis, NormFinder ranked SDHA, B2M, and ACTB as the most stable ( = 0.248, 0.300, and 0.317, resp.) for GSCs. TBP, SDHA, and ALAS1 were the most stable for dGSC ( = 0.222, 0.234, and 0.335, resp.)  TBP, and RPS18 for ASCs ( = 0.063, 0.224, and 0.246, resp.) ( Table 6). Further analysis was realized in order to select the most appropriate reference genes to compare gene expression of GSCs to dGSCs and to ASCs. For this purpose, two pairs of subgroups were created: GSCs versus dGSCs and GSCs versus ASCs. NormFinder calculated intra-and intergroup variances, as shown in charts in which error bars presented intragroup variances and bars referred to intergroup variances (Figures 6(a) and 6(b)).
As for GSCs versus dGSCs, although ALAS1 and UBC displayed low intergroup variances, ALAS1 had a high intragroup variance in GSCs and UBC in dGSC, hence not suitable for the study (Figure 6(a)). RPII had the highest intergroup variance in GSCs and dGSCs and also the highest intragroup variance in GSCs; this gene was then ranked as the least stable one. SDHA and ACTB were selected by the algorithm as the most suitable to compare GSCs to dGSCs, with the lowest intergroup variances and intragroup variances close to zero.
Concerning GSCs versus ASCs, GAPDH and RPS18 had low intervariances, but their intragroup variances were high: GAPDH in both GSCs and ASCs and RPS18 in GSCs. Hence, they were not the most stable genes to study these subgroups. RPII had the highest intergroup variance and was the least stable gene. SDHA and B2M had the lowest intergroup variances and were also stable inside each group, with the lowest intragroup variances. Consequently, they were selected to compare gene expression between GSCs and ASCs.
Thus, NormFinder allowed us to confirm gene ranking of GeNorm for the three groups of samples, GSCs, dGSCs, and ASCs and, moreover, enabled us to select the most   Table 7.

Effect of the Choice of Stable HKGs.
In order to validate the importance of choosing the optimal reference genes to normalize RT-qPCR data, two analyses of gene expression of GSCs were performed based on the former results. The first analysis compared four random samples of GSC, GSC1, GSC3, GSC4, and GSC6, and their expression of a target gene, Collagen 1 alpha 1 (COLL1A1  Figure 6: Determination of the most stable reference genes with NormFinder For GSC versus dGSC (a) and GSC versus ASC (b). Bars represent intergroup variances, while error bars represented the average of intragroup variances. Ideal reference gene had intergroup variation as close to zero as possible and error bars as small as possible. fold expression with no normalization of this gene showed differences between the four samples: GSC6 had the highest relative fold expression (1.00 ± 0.028), followed by GSC3 (0.66 ± 0.001), GSC4 (0.38 ± 0.001), and GSC1 (0.11 ± 0.003). GSC3 had a 6-fold higher expression than GSC1 (Figure 7(a)). Data were then normalized with three groups of genes: SDHA, ACTB, and B2M were the most stable genes as suggested by our present study, GAPDH was the most commonly used reference gene, and ALAS1 and RPS18 were the least stable ones according to our results. The normalization with SDHA, ACTB, and B2M showed a different ranking: GSC3 (2.21 ± 0.021) had the highest normalized fold expression, followed by GSC4 (1.31 ± 0.015), GSC6 (1.09 ± 0.141), and GSC1 (0.62 ± 0.023). In this ranking, GSC3 had a 3.5-fold higher expression than GSC1 (Figure 7(b)). When normalized with GAPDH, the ranking was the same as the normalization with the stable HKGs, but GSC3 (1.83 ± 0.020) had an 8.5-fold higher expression than GSC1 (0.21 ± 0.014) (Figure 7(c)). At last, the normalization with ALAS1 and RPII showed a different ranking, unlike the former analyses, GSC1 (38.98 ± 0.295) had the highest normalized expression, followed by GSC3 (5.02 ± 0.097), GSC4 (1.37 ± 0.070), and GSC6 (1.01 ± 0.300). GSC1 had a 0.15-fold expression as compared to GSC3. The second analysis was performed on samples of dGSCs ( = 6) and compared to ASCs ( = 6) after osteogenic differentiation. Expression level of RUNX2 was assessed on different time points: day 0 (D0), day 7 (D7) and day 14 (D14) after osteogenic differentiation, confirmed by the microscopic observation of mineral nodules formation and Alizarin Red S staining (Supplementary Material 3). The normalization was performed against the two most stable reference genes and compared to the least stable one as shown in Figure 8. For dGSCs, the normalization to SDHA and TBP showed a significant increase of RUNX2 expression at D7 and D14 when compared to the relative gene expression ( Figure 8(a)). When data were normalized with RPS18, the least stable reference gene for dGSCs, there was a decrease in the expression of RUNX2, with a statistically significant difference when compared to SDHA and TBP. Concerning ASCs, the normalization with SDHA and HPRT1 showed also an increasing expression level of RUNX2 at D7 and D14 and an insignificant increase at D7 of this target gene when normalized with ACTB, the least stable reference gene (Figure 8(b)). This confirms the importance and the dramatic effect of the choice of suitable HKGs for GSCs before and after osteogenic differentiation.

Discussion
This study was performed to confirm previously wellrecognized properties of GSCs [1] as well as their osteogenic potential [2]. GSCs might be an interesting tool in human cell therapy and especially in bone regeneration of craniofacial bones that share the same embryonic origin. Our study compared human GSCs to adipose-derived mesenchymal stem cells (ASCs), as they have been well studied in their differentiation capacities, thoroughly highlighted in vitro and in vivo [25]. The osteogenic differentiation capacity of GSCs needs to be supported by q-PCR, which measures the relative expression of genes implicated in osteogenic differentiation. q-PCR is a very accurate method for gene expression analysis; however, data analyses are often variable. This is due to many factors linked to each step of the investigations, from harvesting cells to cDNA synthesis. That is why q-PCR results have to be normalized. The validation of stable and adequate reference genes, to which data are normalized, is the most commonly used method and may allow bypassing variability factors. The choice of a reference gene must be accurate, because any improper HKG may give false interpretation of q-PCR results [26]. GAPDH is the most used reference gene, but recent studies showed that this gene is not suitable for all tissues or cell types or experimental settings [6]. An ideal HKG is expected to have a constant level of expression in all cell types, at all-time points, and under different experimental conditions. Studies have been carried out to reach a "universal HKG" in all species but, to our knowledge, not one was found [13]. To solve this problem of choosing suitable HKGs, many algorithms were generated by mathematicians and statisticians in order to select accurately the most stable reference genes.
In this work, we aimed to select the most appropriate HKGs to study human GSCs and their osteogenic differentiation in vitro in order to use them for clinical application. 10 reference genes were studied on GSCs, dGSCs, and ASCs (mesenchymal stem cells, with known osteogenic capacity) and analysed by four algorithms: ΔCt comparative method, GeNorm, BestKeeper, and NormFinder.
Experiments were rigorously conducted to obtain comparable and reproducible results; a special attention was given to prepare accurately with the same conditions growing and differentiation media while culturing different stem cells. The used RNA extraction kit was appropriate to cell culture and the products were carefully used at the adequate amount corresponding to cells number. RT q-PCR reactions were performed according to the MIQE Guidelines (Supplementary Material 4) [7]. Consequently, in our study, despite a slight dissimilarity in the ranking of HKGs among the four algorithms, results showed a highly homogenous reproducibility; indeed, a similar tendency was found for the three most stable genes and an overall comparable order of genes.
ΔCt comparative method and GeNorm were used for intergene study; ΔCt comparative method identifies the most reliable reference gene as the one that varies the less when compared to all the other genes among all the sample. This method needs no high-level mathematical methodology and is ideal for the nonspecialist to determine the best reference genes [13]. GeNorm software not only defines the optimal combination of genes or the ideal "pair of genes" with low variability, but also proposes the minimal number of stable HKGs. These two methods can be used if the amount of the starting material of studied samples is not enough to study many HKGs, because they are minimally affected by expression intensity [23]. However, they do not take into account coregulated genes and if used alone they may lead to errors. These algorithms showed the same ranking of the ten reference genes in the three groups. SDHA, ACTB, and B2M were the most stable to study GSCs; TBP, SDHA, and ALAS1 for dGSCs; and SDHA, HPRT1, and TBP for ASCs.
For more accuracy, we used BestKeeper and NormFinder, which study also intragene variability. BestKeeper is an interesting tool, because it considers not only intra-and intergroup variation, but also sample integrity. Our samples seemed to have conserved their integrity when applying this analysis throughout all the groups (22 samples of GSCs, dGSCs, and ASCs). NormFinder is considered to be an extra control for GeNorm [22] because it uses a different algorithm and relies on calculating the stability value for each HKG. These two methods seem to be very sensitive to the size of samples and to expression intensity. The ranking of stable reference genes was the same as ΔCt comparative method and GeNorm for GSCs but slightly different for dGSCs and ASCs. BestKeeper ranked B2M, SDHA, and HPRT1 as the most stable for dGSCs and HPRT1, SDHA, and RPS18 for ASCs. The least stable genes were ranked differently than ΔCt comparative method and GeNorm for the three groups of samples. NormFinder ranking was less different from the two precedent algorithms, and the only difference was noticed on ASCs ranking. HPRT1 was ranked 8th rather than with the three most stable reference genes. This can be due to the decreased sample size of ASC ( = 6) which may affect BestKeeper and NormFinder analyses.
Ideally, the ranking of stable HKGs must rely on the results found by the four methods, because the reproducibility and accuracy increase, as the mathematical and statistical bases are different. RefFinder by calculating the geometric mean of overall ordering may help us to identify more easily the stable reference genes; however, it does not take into account q-PCR efficiencies of different primers and only uses Ct values.
Finally, when we analysed each group of samples separately, we found that SDHA, B2M, and ACTB were ranked as the most stable for GSCs; TBP, SDHA, and ALAS1 for dGSCs; and SDHA and HPRT1 for ASCs.
The analysis of gene expression for GSCs either differentiated or not showed that the normalization with random reference genes leads to errors. Our results regarding osteogenic differentiation of GSCs were in agreement with previous publications; TBP was found as a suitable reference gene to study osteogenic differentiation in stem cells [24]. Regardless of the origin of stem cells and the stage of differentiation, SDHA seems to be a good reference gene, as it was stable in all conditions but is not enough alone to study GSC. Ribosomal genes were the least stable, but they might be useful in other conditions [21].
This present study provides methods to determine suitable reference genes. These tools are crucial to studying further the GSC properties and compare the stem cell lineages. Indeed, the osteogenic differentiation can be thoroughly investigated under different conditions (conditioned media, biomaterials, etc.) in order to use GSCs in human bone regeneration.

Conclusion
Regardless of the algorithm used in this study, all of the software used has ranked the same set of reference genes as the most stable. Finally, based on our results, we recommend the use of SDHA, B2M, and ACTB for GSCs; TBP, SDHA, and ALAS1 for dGSCs; and SDHA and HPRT1 for ASCs study. We stress on the importance to select the most suitable HKGs before conducting a study of q-PCR of any stem cell type.