TBP, PPIA, YWHAZ and EF1A1 Are the Most Stably Expressed Genes during Osteogenic Differentiation

RT-qPCR is the gold standard and the most commonly used method for measuring gene expression. Selection of appropriate reference gene(s) for normalization is a crucial part of RT-qPCR experimental design, which allows accurate quantification and reliability of the results. Because there is no universal reference gene and even commonly used housekeeping genes’ expression can vary under certain conditions, careful selection of an appropriate internal control must be performed for each cell type or tissue and experimental design. The aim of this study was to identify the most stable reference genes during osteogenic differentiation of the human osteosarcoma cell lines MG-63, HOS, and SaOS-2 using the geNorm, NormFinder, and BestKeeper statistical algorithms. Our results show that TBP, PPIA, YWHAZ, and EF1A1 are the most stably expressed genes, while ACTB, and 18S rRNA expressions are most variable. These data provide a basis for future RT-qPCR normalizations when studying gene expression during osteogenic differentiation, for example, in studies of osteoporosis and other bone diseases.


Introduction
Osteoporosis is a common age-related bone disease characterized by low bone mineral density, increased bone fragility, and increased risk of bone fracture upon low-energy trauma. It represents a significant healthcare and economic burden, and can have a detrimental impact on the quality of life and life-expectancy of patients due to bone fractures [1]. With aging populations and new developments in research strategies, osteoporosis has gained more attention in the past decade. Several genome-wide association studies have proposed novel genes and proteins involved in bone biology, followed by in vitro and in vivo functional characterization of candidate genes and proposal of novel mechanisms of regulation of bone cell proliferation, differentiation, and function. This has significantly advanced our understanding of pathological changes in bone mineral density regulation. However, the role of many genes has not yet been explored [2,3].
Most of these studies have been performed on different in vitro cell models of osteoblast origin capable of differentiation and mineralization. This enables close observation of the mechanisms on a cellular level, which is frequently achieved through analysis of gene expression. RT-qPCR is thus the gold standard, a frequently used and widely accessible technique that is applied in almost all studies on osteoblast function and differentiation. However, not all studies account for certain drawbacks of the method or follow MIQE guidelines, which can result in poor inter-assay reproducibility or questionable reliability of results [4].

Results
The purpose of this study was to evaluate the stability of selected candidate reference genes in order to determine the most stably expressed genes in selected cell lines and possibly find a general set of genes suitable for other cells of osteogenic origin during differentiation. To obtain this, MG-63, HOS, and SaOS-2 cell lines were differentiated for up to 35 days and mRNA samples were collected every 7 days of differentiation. Gene expression of ten selected candidate reference genes was measured using RT-qPCR and corresponding quantification cycle (Cq) values were analysed by three types of commonly used software, geNorm [23], NormFinder [24], and BestKeeper [25], which calculate the most stably expressed gene(s) in each set based on specific statistical algorithms. The results were used to obtain comprehensive rankings of the most stable genes for each cell line and a suggested gene ranking for cells of bone origin in general.

Osteoblast Differentiation
To induce osteogenic differentiation and mineralization, MG-63, HOS, and SaOS-2 cells were exposed to osteogenic differentiation media for up to 35 days. SaOS-2 cells were only differentiated for up to 21 days, as the cells began to detach with longer culture times. To confirm mineralization, cells were stained with Alizarin Red S and observed under a microscope. As expected, the amount of stained and mineralized extracellular matrix increased in a time-dependent manner (Figures 1d and S1), which was further confirmed by spectrophotometric analysis of the stained samples (Figure 1a-c). The lowest degree of mineralization was observed in MG-63 cell line, while HOS and SaOS-2 showed considerable Alizarin Red S staining. In the non-differentiated controls, no mineralization of extracellular matrix occurred even after 35 days of culturing. This confirmed that a differentiated functional osteoblast cell type was obtained, allowing us to analyse gene stability during the differentiation process.

Osteoblast Differentiation
To induce osteogenic differentiation and mineralization, MG-63, HOS, and SaOS-2 cells were exposed to osteogenic differentiation media for up to 35 days. SaOS-2 cells were only differentiated for up to 21 days, as the cells began to detach with longer culture times. To confirm mineralization, cells were stained with Alizarin Red S and observed under a microscope. As expected, the amount of stained and mineralized extracellular matrix increased in a time-dependent manner (Figures 1d and S1), which was further confirmed by spectrophotometric analysis of the stained samples (Figure 1a-c). The lowest degree of mineralization was observed in MG-63 cell line, while HOS and SaOS-2 showed considerable Alizarin Red S staining. In the non-differentiated controls, no mineralization of extracellular matrix occurred even after 35 days of culturing. This confirmed that a differentiated functional osteoblast cell type was obtained, allowing us to analyse gene stability during the differentiation process.

Candidate Reference Gene Selection and Primer Validation
To determine the optimal internal control genes for study of differentiation processes, ten candidate reference genes were selected (see Table 4). The selected genes mostly belong to different functional classes, and were selected based on most stably expressed genes detected in other similar studies. The performance of all primers was

Candidate Reference Gene Selection and Primer Validation
To determine the optimal internal control genes for study of differentiation processes, ten candidate reference genes were selected (see Table 1). The selected genes mostly belong to different functional classes, and were selected based on most stably expressed genes detected in other similar studies. The performance of all primers was validated for each cell line. Amplification efficiencies for each gene were calculated from standard curves and ranged from 92.6% to 106.7%, with corresponding correlation coefficients (R 2 ) from 0.9956 to 1 ( Table 2). All primers had efficiencies within the recommended range (100 ± 10%) [26]. For each primer set only one melting curve was detected, confirming primer specificity ( Figure S2). These findings confirmed that the designed primers were suitable for use on all three selected osteoblast cell lines.

Analysis of Stability of Candidate Reference Genes
To determine the expression stability of the selected genes, their expression was determined at different time points during osteogenic differentiation in all three cell lines. Expression profiles were first evaluated by descriptive statistics. In all cell lines, each gene reached the detection threshold at a comparable Cqs, indicating comparable levels of expression in each cell line ( Figure 2). The lowest Cq values were from 10.16 ± 0.39 to 11.16 ± 0.59 for 18S rRNA, while the highest were from 24.94 ± 0.30 to 25.08 ± 0.31 for TBP. Figure 2 shows the distribution of Cq values for each gene at different differentiation time-points and provides an overview of the stability of the measured values during differentiation. Further analysis of stability was performed using three commonly applied statistical algorithms, geNorm, NormFinder, and BestKeeper.

Analysis of Stability of Candidate Reference Genes
To determine the expression stability of the selected genes, their expression was determined at different time points during osteogenic differentiation in all three cell lines. Expression profiles were first evaluated by descriptive statistics. In all cell lines, each gene reached the detection threshold at a comparable Cqs, indicating comparable levels of expression in each cell line ( Figure 2). The lowest Cq values were from 10.16 ± 0.39 to 11.16 ± 0.59 for 18S rRNA, while the highest were from 24.94 ± 0.30 to 25.08 ± 0.31 for TBP. Figure  2 shows the distribution of Cq values for each gene at different differentiation time-points and provides an overview of the stability of the measured values during differentiation. Further analysis of stability was performed using three commonly applied statistical algorithms, geNorm, NormFinder, and BestKeeper.

geNorm Analysis
geNorm's pair-wise approach ranks genes based on their M value, where a lower M value represents higher gene stability. As seen in Table 2, all calculated M values were < 1.5, therefore, all genes in every cell line are considered stable. However, the lowest M  (Table 3 and Figure 3a).
The pairwise variations analysis showed that the pairwise value V2/3 in every cell line is already below 0.15 (0.113 for MG-63, 0.079 for HOS, 0.069 for SaOS-2); therefore, the normalization to two most stable genes is sufficient and the addition of the third normalization gene does not contribute significantly to the reliability of the results (Figure 3b).

NormFinder Analysis
NormFinder uses a model-based approach by which it ranks genes according to their stability value. Again, the lowest stability value is attributed to the most stable gene. This algorithm identified PPIA as the most stable gene in both MG-63 and SaOS-2 (0.079 and 0.054, respectively). These two cell lines also share the least stable gene 18S rRNA (0.248 for MG-63 and 0.237 for SaOS-2). In HOS cell line, NormFinder identified YWHAZ as the most stable (0.114) and RPL13A as the least stable gene (0.478) ( Table 3). Table 3. Stability of selected candidate reference genes as ranked by geNorm, NormFinder, and BestKeeper in MG-63, HOS, and SaOS-2 cell lines. Candidates are listed from top to bottom by decreasing expression stability based on geNorm's M value, NormFinder's stability value, and BestKeeper's SD, CD, and r as well as by geometrical mean of comprehensive ranking.   The pairwise variations analysis showed that the pairwise value V2/3 in every cell line is already below 0.15 (0.113 for MG-63, 0.079 for HOS, 0.069 for SaOS-2); therefore, the normalization to two most stable genes is sufficient and the addition of the third normalization gene does not contribute significantly to the reliability of the results ( Figure  3b).

NormFinder Analysis
NormFinder uses a model-based approach by which it ranks genes according to their stability value. Again, the lowest stability value is attributed to the most stable gene. This algorithm identified PPIA as the most stable gene in both MG-63 and SaOS-2 (0.079 and 0.054, respectively). These two cell lines also share the least stable gene 18S rRNA (0.248 for MG-63 and 0.237 for SaOS-2). In HOS cell line, NormFinder identified YWHAZ as the most stable (0.114) and RPL13A as the least stable gene (0.478) ( Table 2).

BestKeeper Analysis
BestKeeper ranks genes' stability based on their standard deviation (SDCq value; ±Cq), coefficient of variation (CVCq value; % Cq), and correlation coefficient (r). MG-63 and HOS share EF1A1 as a gene with the lowest SD (0.259 for MG-63 and 0.124 for HOS) and ACTB as the gene with the highest SD (1.011 for MG-63 and 0.666 for HOS). In SaOS-2, YWHAZ exhibits the lowest SD (0.209) and RPLP0 the highest SD (0.514). As BestKeeper's requirement for stable gene is SD < 1, all genes, with the exception of ACTB in MG-63, were considered stable.

BestKeeper Analysis
BestKeeper ranks genes' stability based on their standard deviation (SD Cq value ; ±Cq), coefficient of variation (CV Cq value ; % Cq), and correlation coefficient (r). MG-63 and HOS share EF1A1 as a gene with the lowest SD (0.259 for MG-63 and 0.124 for HOS) and ACTB as the gene with the highest SD (1.011 for MG-63 and 0.666 for HOS). In SaOS-2, YWHAZ exhibits the lowest SD (0.209) and RPLP0 the highest SD (0.514). As BestKeeper's requirement for stable gene is SD < 1, all genes, with the exception of ACTB in MG-63, were considered stable.

Comprehensive Ranking
As the three applied algorithms use different logic to determine gene expression stability, each of them provided a different ranking based on the same Cq values for each cell line ( Table 3). As none of these approaches is considered more reliable than the others, we included all of them and performed a comprehensive ranking using two approaches: the geometric mean of the ranks and NormFinder's algorithm for estimating variability between subgroups. Based on the geometric mean comprehensive ranking, the most stably expressed gene in MG-63 is TBP, while YWHAZ is the most stable in both HOS and SaOS-2 ( Table 3). The least stably expressed are ACTB in MG-63, RPL13A in HOS, and 18S rRNA in SaOS-2. On the other hand, NormFinder's model-based algorithm identified GAPDH and EF1A1 as the best combination for normalization (Table S1).
As certain genes were consistently ranked as more stable in all three observed cell lines, we wanted to determine whether certain genes would show general stability across different cell lines, making them suitable for use as a starting point in experiments involving osteoblast differentiation. We therefore combined the comprehensive ranking obtained for each cell line in this study and calculated the geometric means of the ranks for each gene. Using this approach, YWHAZ and TBP were identified as most stable, closely followed by PPIA and EF1A1 (Table 4). To test our assumption and confirm the ranking obtained in this study, we collected the data from other similar studies where gene stability in cells of bone or osteoblast origin was evaluated [12][13][14][20][21][22]27]. Stability ranks obtained in these studies were combined with our data and a comprehensive ranking of each gene was determined using the geometric mean of the ranks. As in our study, PPIA, TBP, and YWHAZ ranked as the most stable and ACTB as the least stable gene (Table S2).

Normalization of Genes Involved in Osteogenic Differentiation against Selected Reference Genes vs. ACTB
To determine the influence of the reference gene, we determined the expression of three commonly used osteogenic differentiation markers, Alkaline phosphatase (ALP), Collagen Type I Alpha 1 (COL1A1), and Runt-related transcription factor 2 (RUNX2), in each cell line: the expression data was normalized using the geometric mean of the relative quantities of TBP, PPIA, YWHAZ, and EF1A1 genes or ACTB gene (Figure 4). Although the shapes of the obtained graph curves are similar and show the same trends, normalization against ACTB significantly increased the normalized expression of all three markers in MG-63 and HOS cells. Interestingly, there were no significant differences in expression obtained in SaOS-2 cells, and the expression of differentiation markers remained mostly stable during the differentiation process, suggesting that these cells were already finally differentiated and only underwent mineralization. This indicates that the choice of reference gene can strongly influence the fold increase of the observed expression changes after normalization. mostly stable during the differentiation process, suggesting that these cells were already finally differentiated and only underwent mineralization. This indicates that the choice of reference gene can strongly influence the fold increase of the observed expression changes after normalization. (Collagen type I alpha 1 chain) (b,e,h), and RUNX2 (Runt-related transcription factor 2) (c,f,i) differentiation markers was analysed at given time points using RT-qPCR. The relative quantity of each marker gene was normalized to ACTB (red) or to geometric mean of TBP, PPIA, YWHAZ, and EF1A1 (black) internal control genes. Results represent the mean and standard error of relative expression for three independent experiments. Statistical differences between the two normalizations were analysed using two-way ANOVA and are denoted as * p ≤ 0.05 and *** p ≤ 0.001. (c,f,i) differentiation markers was analysed at given time points using RT-qPCR. The relative quantity of each marker gene was normalized to ACTB (red) or to geometric mean of TBP, PPIA, YWHAZ, and EF1A1 (black) internal control genes. Results represent the mean and standard error of relative expression for three independent experiments. Statistical differences between the two normalizations were analysed using two-way ANOVA and are denoted as * p ≤ 0.05 and *** p ≤ 0.001.

Discussion
Differentiation of osteoblast-like cells is a highly complex multistep process of cell maturation during which the cells continuously deposit mineralized extracellular matrix [28]. In vitro differentiation is stimulated with dexamethasone to enhance expression of the transcription factors, L-ascorbic acid (Vitamin C) to increase collagen type I secretion, and β-glycerophosphate as a source of phosphate for hydroxyapatite formation and additional stimulation of signalling pathways [29]. The differentiation process can be followed by detecting the mineralization of the extracellular matrix using Alizarin Red S staining, as well as by following the expression of differentiation markers such as osteocalcin, ALP, RUNX2, COL1A1, etc. [22,28,30]. When measuring differentiation markers using RT-qPCR, their expression levels should be carefully normalized to reference genes that remain stably expressed during the differentiation process.
We exposed three commonly used human cell lines of osteogenic origin, MG-63, HOS, and SaOS-2, to osteogenic differentiation media in vitro. MG-63 and HOS cell lines were treated for 35 days, while SaOS-2 differentiation was terminated at day 21 due to the spontaneous detachment of cells from the plates. Time-dependent increase in mineralization was observed in each cell line (Figures 1 and S1). Despite being cultured for only 21 days, SaOS-2 exhibited extensive mineralization even at this time point. While HOS and SaOS-2 cell lines were extensively mineralized as judged by Alizarin Red S, MG-63 cells showed only modest mineralization at day 35 when examined spectrophotometrically. However, when observed under the microscope deposition of minerals was clearly seen and significantly different from the non-differentiated control. Therefore, by the last day of culture all of treated cells exhibited a differentiated/mineralized phenotype.
To determine the best reference genes for studying osteoblast differentiation, we first chose ten genes that are either commonly used as reference genes (GAPDH, 18S rRNA, ACTB) or have been found to be the most stable in other similar studies addressing gene stability in cells of bone/osteoblast origin (PPIA, EF1A1, RPL13A, RPLP0, TBP, HPRT1, YWHAZ). Expression profiling of the selected genes revealed that all of the genes were comparably expressed in all cell lines. While a majority of the genes exhibited similar Cqs (approximately [15][16][17][18][19][20], 18S rRNA had significantly lower Cq (~10) whereas HPRT1 and TBP Cqs exceeded 20. These differences arise from differences in basal expression of the genes, and could be avoided by modifying the amount of starting material. However, we chose to keep the amount of starting material constant in all experiments in order to avoid differences in sample preparation possibly leading to higher inter-experimental variability.
We applied geNorm, NormFinder, and BestKeeper software the selection of the appropriate reference gene. Each of these software types uses a different statistical algorithm; therefore, their integrated use might represent a way to avoid discrepancies in rankings. Such approach is a commonly used for determining the most stable genes [31][32][33]. As expected, we observed differences in the top-ranked genes between the different algorithms. However, in general, certain genes consistently ranked as more stable while others were consistently among the least stable, with the exceptions attributed mainly to BestKeeper's algorithm and ranking. For example, TBP in SaOS-2 drops from second place in Best-Keeper's SD and CV ranking to the tenth place in correlation coefficient (r) ranking. This parameter describes the correlation between a chosen candidate gene and the geometric mean of all housekeeping genes (BestKeeper index). It has been previously reported that correlation coefficient describes the stability of the gene better than SD, as the correlation coefficient correlates each gene with the BestKeeper index and SD can vary with the amount of input material in each sample [26,34]. Not surprisingly, the lowest correlation index was attributed to TBP, which exhibited the highest Cq values (~25), and 18S rRNA, which had the lowest Cq values (~10). On the other hand, the lowest-ranked genes by BestKeeper's correlation coefficient agree with the rankings by geNorm and NormFinder. Generally, the least stable genes correlate well between the algorithms.
The geNorm algorithm ranks genes based on their M value, which is inversely proportional to their stability. However, this algorithm cannot distinguish between the two most stable genes, and it is biased towards genes which could be co-regulated. Taking this issue into account, NormFinder and BestKeeper are considered to be more robust, as their algorithms are less subject to co-regulation [24,25]. When comparing the top-ranked genes by NormFinder and geNorm, similar rankings can be observed. At this point, it must be noted that based on the cut-off values of the chosen algorithms, namely, geNorm's M value < 1.5 and BestKeeper's SD < 1, all genes were considered stable with the exception of ACTB's SD in MG-63 (SD = 1.011). The differences in ranking values between consecutivelyranked genes are relatively small, showing similar degrees of stability. This outcome is not surprising, as the initial selection of most of the candidate genes was based on the findings of previous studies which identified them as stable in differentiation processes of cells of similar origin [13,[20][21][22][34][35][36][37].
Of the three algorithms, geNorm is the only one that proposes the number of genes required for normalization. Using pairwise variation analysis, it calculates the V value which describes the relevance of the additional gene to the normalization process. It is proposed that 0.15 is a cut-off value below which an additional gene does not contribute significantly to normalization [23]. In every cell line, geNorm calculated that two genes were sufficient for normalization. Using a comprehensive ranking by geometric mean which integrated data from the different algorithms, we chose two appropriate reference genes for normalization during osteoblastic differentiation in each cell line; these were TBP and PPIA in MG-63, YWHAZ and EF1A1 in HOS, and YWHAZ and PPIA in SaOS-2.
Even though MG-63, HOS, and SaOS-2 cells were of human osteosarcoma origin and were exposed to identical differentiation stimuli, different genes were identified as the most stable in these cell lines. These differences likely arise from their internal characteristics, e.g., differentiation stage and genetic instabilities, which offer different starting point for differentiation and mineralization. As observed in the comparison of the rankings of different algorithms for each cell line, certain genes were consistently ranked as more stable when comparing different cell lines. As we wanted to propose a group of stable reference genes for reliable use in studying osteoblast differentiation and mineralization, we compared NormFinder's inbuilt algorithm for comparison of different groups to the geometric mean of all compared ranks for each gene in this study.
NormFinder's model-based algorithm identified GAPDH and EF1A1 as the best combination for normalization (Table S1). However, GAPDH exhibits relatively high SDs and CVs (especially in HOS and SaOS-2 cell lines) and EF1A1 shows unfavourable correlation with the BestKeeper index. Additionally, selection of the best genes for normalization by NormFinder has previously been shown to be unreliable [38]. On the other hand, the geometric mean of the ranks for all three cell lines identified YWHAZ and TBP to be the most stable genes (Table 4).
Additionally, we combined our data with comprehensive rankings from other studies, selecting only cells of bone or osteoblast origin (Table S2) [12][13][14][20][21][22]27]. This ranking is comparable with the overall rankings obtained in this study, suggesting that YWHAZ, TBP, PPIA, and EF1A1 are stably expressed during osteoblast differentiation. Small differences in their geometrical means confirm that their expression stabilities are very similar.
To assess the selection of the most stable internal control genes, the same samples were analysed for expression of the differentiation markers ALP, RUNX2, and COL1A1 ( Figure 4). The obtained expression data (relative quantity) were normalized to the ACTB reference gene (a commonly used internal control) and to the geometric mean of the four internal control genes detected as most stable in this study (TBP, PPIA, YWHAZ, EF1A1). The results showed that for the MG-63 and HOS cell lines normalization to ACTB can considerably increase the obtained values, especially in later differentiation stages. This suggests that ACTB expression decreases during differentiation compared to the geometric means of expression of TBP, PPIA, YWHAZ, and EF1A1. Interestingly, in SaOS-2 cells both normalizations showed similar results. Furthermore, there were no changes in expression of observed differentiation markers despite significant mineralization, a phenomenon which has been observed previously [39,40]. This indicates that the SaOS-2 cell line consists of terminally differentiated osteoblast cells and that differentiation media only induces mineralization (Figures 1 and S1). The lack of changes in the differentiation state of these cells could explain the stability of ACTB expression during mineralization as well. These results indicate that the choice of internal control can have a significant impact on the obtained results. If the expression dynamics of certain genes are compared through several cell lines, the stability of the internal control genes used should be comparable in all cell lines.
Even though we recommend careful selection and validation of reference genes for each assay and cell line, any of YWHAZ, TBP, PPIA, and EF1A1 or their combination (we recommend at least two) could be used as starting reference genes in assays of osteoblastic differentiation when using other types of osteosarcoma cells or cell lines. We do not recommend the use of ACTB and 18S rRNA, as they are frequently ranked as the least stable genes.
For osteoblast differentiation, cells were seeded in 12-well plates. Following 24 h rest, differentiation was triggered by changing to osteogenic medium: full growth medium supplemented with 100 nM Dexamethasone (Sigma, Merck, Kenilworth, NJ, USA), 50 µg/mL 2-Phospho-L-ascorbic acid trisodium salt (Sigma), and 5 mM β-Glycerophosphate disodium salt hydrate (Sigma) for MG-63 and SaOS-2 cells and full growth medium supplemented with 100 nM Dexamethasone (Sigma), 50 µg/mL 2-Phospho-L-ascorbic acid trisodium salt (Sigma), and 10 mM β-Glycerophosphate disodium salt hydrate (Sigma) for HOS cells. Medium was changed every 3-4 days. Nondifferentiated controls were maintained in full growth medium. Samples were collected at day 0 (nondifferentiated control day 0) and then every 7 days of culture to obtain five time-points (days 0, 7, 21, and 35 for HOS and MG-63; days 0, 7, 14, and 21 for SaOS-2). Nondifferentiated controls were collected on the last day of the experiment. All experiments were performed three independent times in experimental duplicates.

Alizarin Red S Staining
To confirm mineralization in the differentiated samples, Alizarin staining was performed. Cells were washed 3× with PBS and fixed with 4% formalin for 10 min at room temperature. Formalin was removed and the cells were washed 3× with distilled water. Cells were then stained with 2% Alizarin red S solution (pH 4.2) (Sigma) for 30 min and washed 5× with distilled water. Plates were analysed under the microscope. For spectrophotometric analysis, Alizarin was dissolved in 10% acetic acid and measured at 405 nm. Alizarin concentration was calculated from a standard curve.

RNA Isolation and Reverse Transcription
Total RNA was isolated using a peqGOLD Total RNA Kit (VWR, Radnor, PA, USA) following the manufacturer's instructions. Yield and purity of RNA were determined with a NanoDrop One C spectrophotometer (Thermo Scientific, Waltham, MA, USA); 3.0 µg of total RNA was reverse transcribed in a 60 µL reaction using a High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Thermo Scientific, Waltham, MA, USA) following the manufacturer's instructions.

Gene Selection, Primer Design, and Validation
The reference genes used in this study were selected based on genes frequently used for normalization in differentiation experiments on these cell lines (ACTB, GAPDH, 18S rRNA) and genes and that were stable in studies analysing the stability of internal controls in bone-derived cells or cell lines [13,[20][21][22][34][35][36][37] (Table 1). Genes and their protein products cover different functions in cells, from structural proteins to enzymes, ribosomal proteins, and different transcription factors. Primer BLAST was used for primer design. If possible, primers were selected to span an exon-exon junction, to have an optimal melting temperature around 60 • C, and to result in a PCR product length of 70-170 bp. Primers were ordered from Macrogen Europe B.V (Amsterdam, the Netherlands, EU). Primers were validated with a standard curve to determine their efficiency. The specificity of the primers was checked using melting curves.

RT-qPCR
RT-qPCR was performed using 5× Hot FirePol EvaGreen qPCR Supermix (Solis, BioDyne, Tartu, Estonia) following the manufacturer's instructions. For each reaction, 2.5 ng of cDNA was used in a 15 µL reaction. The amplification programme consisted of initial denaturation at 95 • C for 15 min and 40 cycles of 95 • C for 5 s (denaturation), 60 • C for 20 s (annealing), and 20 s at 72 • C (elongation). The amplification programme was terminated with a thermal denaturing cycle to obtain the melting (thermal dissociation) curve. The RT-qPCR was performed on a LightCycler 480 (Roche Diagnostics, Rotkreuz, Switzerland). All samples were quantified in duplicates.

Primer Efficiency
A seven-point serial dilution of pooled cDNA samples was prepared to obtain a relative standard curve for each gene in each cell line. The standard curve was obtained by plotting the Cq (y axis) against the logarithm of the total cDNA concentration (x axis). The efficiency of each primer was calculated using Equation (1): where E is primer efficiency and slope is a log-linear relationship between Cq values and the cDNA concentrations of the standard curve [4].
geNorm is a Microsoft Excel-based tool which calculates gene expression stability based on the geometric mean of candidate genes' 2 −∆Cq , calculated using Equation (2): where min Cq is the minimal Cq value of the sample set and sample Cq is the Cq of the selected sample in the sample set. geNorm ranks candidate genes based on an internal control gene stability measure called the M value, which describes how one gene varies with respect to other candidate genes. The M value is defined as the average pairwise variation of a particular gene with all other control genes. Lower M values indicate more stably expressed genes, and genes with an M value < 1.5 are considered stable. Using stepwise exclusion of the least stable gene followed by the recalculation of the M values of the remaining genes, this algorithm can identify the two most stably expressed reference genes. In addition, geNorm calculates the pairwise variation of the selected reference genes with the others and determines how many genes are needed for accurate normalization. For this, a V value is determined between two subsequent normalization steps (V = Vn/(Vn + 1)). When V is less than 0.15, additional genes do not contribute significantly to normalization [23].
NormFinder is a model-based approach that uses 2 −∆Cq values [Equation (2)] to estimate intra-and intergroup variations. It then combines the two to calculate a stability value, with the lowest value indicating the most stably expressed gene. Because of the model-based approach, it is less subject to gene co-regulation than the pair-wise variation approach of geNorm [24].
BestKeeper estimates inter-gene relations of candidate reference genes using pair-wise correlation analysis. It uses raw Cq and PCR efficiencies to calculate standard deviation (SD Cq value ) and coefficient of variation (CV Cq value ), and calculates a BestKeeper index (geometric mean) from the candidate genes' Cqs. Then, it calculates the pair-wise correlation between each candidate gene and the BestKeeper index as the correlation coefficient. The most stably expressed genes have the lowest SD Cq value and CV Cq value and the highest correlation coefficient. Generally, genes are considered as stable when SD Cq value < 1 [25].
To assess the overall stability of reference genes, an appropriate weight based on the ranking by each algorithm was attributed to each gene. The comprehensive ranking was performed using the geometric mean of the weights, and was calculated separately for each cell line as well as for overall stability in all cell lines combined.
Expression of differentiation markers (ALP, RUNX2, COL1A1) ( Table 5) was evaluated by calculating the relative quantity (Rq) for each sample based on the PCR efficiency of each gene (E) using Equation (3): The obtained relative quantities of marker gene expression were normalized to the relative quantity of the ACTB gene or to the geometric mean of the TBP, PPIA, EF1A1, and YWHAZ internal control genes. Statistical differences between normalization strategies were determined at each time point using two-way Anova (GraphPad Prism 8.0.1). Table 5. Osteogenic differentiation-involved genes, their function, and the primer sequences used in this study.

Symbol Accession Number
Gene

Conclusions
In this study, we identified the most stable reference genes for accurate normalization of RT-qPCR results of osteogenic differentiation of the human osteosarcoma cell lines MG-63, HOS, and SaOS-2. Using the geNorm, NormFinder, and BestKeeper programs, we identified that each cell line requires two reference genes for normalization, which are TBP and PPIA in MG-63, YWHAZ and EF1A1 in HOS, and YWHAZ and PPIA in SaOS-2. Additionally, we compared our results with the outcomes of previous studies focused on identification of reference genes for other bone-derived cell lines. Our findings strongly correlate with them, confirming that TBP, PPIA, YWHAZ, and EF1A1 are stably expressed during osteogenic processes in different cell types and cell lines of bone origin. Certain commonly used reference genes such as ACTB and 18S rRNA are not recommended for normalization during osteogenic differentiation, as they were proven to be the least stably expressed in this study as well as in others. The internal control genes identified in this study will help researchers design better RT-qPCR experiments for the study of different processes during osteoblast differentiation of the MG-63, HOS, and SaOS-2 cell lines. Moreover, we propose a group of stably expressed genes (TBP, PPIA, YWHAZ, and EF1A1) that can be used as a starting point for further optimization of RT-qPCR experiments involving differentiation of osteoblasts in general.