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.

1. 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, real-time 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 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 [1114]. 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.

2. 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.

2.1. Isolation and Growth of Stem Cells: Gingival Stem Cells (GSCs)

GSC cultures were obtained from six healthy donors (), 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) [1618]. 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% CO2 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].

2.2. Adipose Derived Stem Cells (ASCs; Purchased from ZenBio Company, USA)

Six samples of ASCs (), mesenchymal stem cells of adipose tissue origin, used as control, were also cultured in the same conditions, in early passages (1–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.

2.3. Osteoinduction

GSC () and ASC () were cultured in 6-well plates. When attempting 70 to 80% of confluence, osteogenic differentiation was induced by DMEM-LG supplemented with 20% of FCS, 10 mL/L of Penicillin-Streptomycin (5 UI/mL), 1% of NEAA, 2.5 mg/L of Amphotericin B (250 μg/mL), adding to 50 μg/mL of L-Ascorbic Acid 2-Phosphate, 100 nM of dexamethasone, and 10 mM of β-glycerophosphate [19]. Culture media were changed every 72 hours, by adding L-Ascorbic Acid 2-Phosphate. Dexamethasone was supplemented every 7 days. After 21 days, differentiated GSC (dGSC) and ASC were either harvested after Trypsin-EDTA treatment and centrifugation, frozen at −80°C for RNA extraction, or fixed in PFA 4%/Sucrose 5% and conserved at 4°C for staining.

2.4. 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 , ASC , and differentiated GSC ), 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 / 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)).

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.

2.6. Selection of Reference Genes

Primer sequences for the 10 candidate housekeeping genes in the study of human stem cells, TBP, SDHA, HPRT1, GAPDH, ALAS1, ACTIN, RPS18, RPII, UBC, and B2M, were obtained from NCBI primer-blast tool, USA (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 = () × 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.

2.8. 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 detection and is therefore reversely correlated to the input amount of total RNA [20]. These values were used in Ct and BestKeeper algorithms. Values of relative gene expression or normalized values () were also calculated for use in GeNorm and NormFinder analyses (see Supplementary Material 1 in Supplementary Material available online at http://dx.doi.org/10.1155/2016/6261490). For each HKG, Ct values were logarithmically transformed for all the samples; the lowest Ct was used as a calibrator against which change is given. For each sample, Ct (the subtraction of cycle threshold between the studied sample and the calibrator) is calculated. The investigators use the efficiency of each primer, and the formula is the following:

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].

2.8.2. 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 [22]. For each two candidate genes (gene and gene ), using normalized values (, here ), the algorithm calculates the array of log2 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 toward gene for all samples; is standard deviation of values for gene 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/.

2.8.3. 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.

2.8.4. 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 intra- and 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 (). 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.

3. Results

3.1. 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 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 / 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)).

3.2. The Choice of the Most Reliable HKGs for GSCs, dGSCs, and ASCs

In order to select the most appropriate and stable HKGs to study GSC and their osteogenic differentiation, expression of ten reference genes chosen according to previous published studies [15, 23, 24] was assessed. These reference genes were TBP (TATA-binding protein), SDHA (succinate dehydrogenase complex, subunit A, flavoprotein), HPRT1 (hypoxanthine guanine phosphoribosyltransferase I), GAPDH (glyceraldehyde-3-phosphate dehydrogenase), RPS18 (40S ribosomal protein S18), ALAS1 (5-Aminolaevulinate synthase), ACTB (Beta-actin (β-actin)), B2M (Beta-2-microglobulin), UBC (Ubiquitin C), and RPII (50S ribosomal protein L9). Gene information are available and primer’s efficiencies were between 82.5% and 103% (Table 1).

RT-qPCR were performed for 35 cycles in GSC samples (), GSCs after osteogenic induction (referred to as differentiated GSCs: dGSC) (), and compared to ASCs (). A melting curve analysis showed a single product at the expected melt temperature for each reference gene (Figure 2(a)). Results were analysed with a logarithmic scale and a threshold level was manually selected for each reference gene to detect the expression levels as cycle threshold (Ct) values for all samples (Figure 2(b)). Ct values were collected for each group of samples (Figures 2(c), 2(d), and 2(e)). In all samples, GAPDH, RPS18, B2M, and ACTB had Ct values below 20 cycles and hence had an early expression, unlike SDHA, TBP, HPRT1, ALAS1, UBC, and RPII that had lowest expression values. The normalized values (): were calculated for each sample using Ct values (Supplementary Material 1). Data were treated with four algorithms: the Ct comparative method [13] and GeNorm [14] for inter-gene relation, BestKeeper [12] and NormFinder [11] for intra- and intergene correlation (Table 2).

3.3. 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 box-and-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.

Finally, the intergene analysis by Ct comparative method ended with a recommendation to normalize GSCs data using SDHA, ACTB, and B2M, dGSCs with TBP, SDHA, and ALAS1, and ASCs with SDHA, HPRT1, and TBP.

3.4. 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 (), the most stable genes were SDHA, B2M, and ACTB, (, 0.919, and 0.924, resp.). RPII was the least stable () and excluded from analysis (). 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 (), 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.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 (), SDHA, HPRT1, and TBP were the most stable (, 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 ( and ) suggested that the optimal number of reference genes to study teach group of stem cells was three HKGs for GSCs and dGSCs ( 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 ( below 0.15) (Figure 4(b3)).

Finally, the intergene analysis by GeNorm recommended SDHA/ACTB and UBC to study GSCs, TBP/ALAS1 and SDHA for dGSCs, and GAPDH/B2M or SDHA/HPRT1 to study ASCs (Table 4).

3.5. BestKeeper Analysis

Raw Ct values of each group (GSCs (), dGSCs (), and ASCs ()) 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 and CV , 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.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(a)). Among the remaining genes, SDHA had the lowest variation (SD = 0.79, CV = 3.55) and the best correlation with BestKeeper index (, ) followed by B2M (, ), ACTB (, ), and GAPDH (, ) as summarized in Table 5(a).

For dGSCs (), all HKGs had low variances, with the least variability for TBP (SD = 0.37, CV = 1.40) and ACTB (SD = 0.46, CV = 2.48) and the highest variability for RPII (SD = 0.79, CV = 3.27) and UBC (SD = 0.80, CV = 2.92). Then the correlation with BestKeeper index excluded GAPDH (), RPS18 (), ACTB (), and ALAS1 (). Gene ranking from the most correlated was B2M (, ), SDHA (, ), HPRT1 (, ), and TBP (, ). RPII and UBC had a high correlation with BestKeeper index (, and , , resp.); however, they were the least stable, as shown formerly with the highest SD and CV values (Figure 5(b)). Consequently, the final ranking was B2M, SDHA, HPRT1, and TBP (Table 5(b)).

For ASCs (), ALAS 1 was first excluded (SD = 1.01, CV = 4.23). Then, the correlation with BestKeeper index showed a low for ACTB (), GAPDH (), B2M (), and UBC (), which were also excluded (Figure 5(c)). The most correlated HKGs with BestKeeper index were HPRT1 (; ), SDHA (; ), RPS18 (; ), and RPII (; ) as shown in Table 5(c).

Based on these results, the most reliable HKGs with BestKeeper software were SDHA, B2M, and ACTB for GSCs, B2M and SDHA for dGSCs, and HPRT1, SDHA, and RPS18 for ASCs.

3.6. 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.300, and 0.317, resp.) for GSCs. TBP, SDHA, and ALAS1 were the most stable for dGSC (, 0.234, and 0.335, resp.), and SDHA, TBP, and RPS18 for ASCs (, 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 appropriate genes to compare GSCs versus dGSCs: SDHA and ACTB, and GSCs versus ASCs: SDHA and B2M.

3.7. Final Ranking of HKGs for GSCs, dGSCs, and ASCs

Based on the four algorithms, an additional control was performed with RefFinder, a tool that uses only Ct values and calculates the geometric mean of overall ordering of all reference genes and suggested a final ranking of the most stable housekeeping genes. Results are summarized in Table 7.

3.8. 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). The relative fold expression with no normalization of this gene showed differences between the four samples: GSC6 had the highest relative fold expression (), followed by GSC3 (), GSC4 (), and GSC1 (). 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 () had the highest normalized fold expression, followed by GSC4 (), GSC6 (), and GSC1 (). 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 () had an 8.5-fold higher expression than GSC1 () (Figure 7(c)). At last, the normalization with ALAS1 and RPII showed a different ranking, unlike the former analyses, GSC1 () had the highest normalized expression, followed by GSC3 (), GSC4 (), and GSC6 (). GSC1 had a 0.15-fold expression as compared to GSC3.

The second analysis was performed on samples of dGSCs () and compared to ASCs () 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.

4. Discussion

This study was performed to confirm previously well-recognized 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 () 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.

5. 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.


ASCs:Adipose-derived stem cells
dGSCs:Differentiated gingival stem cells
GSCs:Gingival stem cells
MIQE:Minimum Information for Publication of Quantitative Real-Time PCR Experiments
qPCR:Quantitative Polymerase Chain Reaction
:Relative expression.

Conflict of Interests

The authors declare that they have no competing interests.

Authors’ Contribution

Ihsène Taïhi, a Ph.D. student, helped to design the study, performed cell culture of GSC, dGSCs, and ASC, RNA extraction, and cDNA synthesis, carried out the qPCR studies, executed the CT, GeNorm, BestKeeper, and NormFinder analyses, and drafted the paper. Benjamin Philippe Fournier designed or extracted primers from previous publications for qPCR, conceived the study, selected algorithms, provided daily supervision for Ihsène Taïhi, and contributed in drafting the paper. Ali Nassif participated in the analyses of data with different algorithms and drafting the paper. Tsouria Berbar participated in RNA extraction, DNA synthesis, and the confirmation of qPCR results. Juliane Isaac, Ariane Berdal, and Bruno Gogly conceived the study in its design, coordinated all efforts, and finalized the paper. All authors read and approved the final paper. Ali Nassif and Tsouria Berbar equally contributed to this work.


This work has been performed at the Cordeliers Research Center, in the Laboratory of Molecular Oral Physiopathology, INSERM UMRS 1138, in collaboration with the Odontology Department of the Hospital Complex Henri-Mondor Albert-Chenevier, Créteil, France. The authors thank the “Gueules Cassées” Foundation for their support of this work.

Supplementary Materials

Supplementary Material 1: Three tables showing normalized values (R) of the three groups of samples GSCs, dGSCs and ASCs with the ten studied HKGs. R were calculated using formula and depend on E the efficiency of each primer (1).

Supplementary Material 2: Two tables describing the ΔCt comprative method for dGSCs and ASCs and the ranking of the ten HKGs according to StdDev mean of each gene.

Supplementary Material 3: Mineral nodules were observed for dGSCs at different time points D7, D14 and D21. This confirmed the osteogenic differentiation in accordance with the increase of RUNX2 gene expression. Alizarin Red S staining at D21 confirmed the calcium content of these nodules.

Supplementary Material 4: MIQE checklist was almost respected by the different steps of our study, from collecting samples to gene expression analyses, which is in accordance with the accuracy of our results.

  1. Supplementary Materials