Identification of valid reference genes during the differentiation of human myoblasts

Background Analysis of RNA expression using real-time PCR (qRT-PCR) traditionally includes reference genes (RG) as an internal control. This practice is being questioned as it becomes increasingly clear that RG may vary considerably under certain experimental conditions. Thus, the validity of a particular RG must be determined for each experimental setting. We used qRT-PCR to measure the levels of six RG, which have been reported in the literature to be invariant. The RG were analyzed in human myoblast cultures under differentiation conditions. We examined the expression by qRT-PCR of mRNA encoding Beta-actin (ACTB), Beta-2-microglobulin (B2M), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), peptidylprolyl isomerase A (PPIA), TATA box binding protein (TBP) and ribosomal protein, large, P0 (RPLPO). The mRNA expression of the following genes of interest (GOI) were analyzed: skeletal muscle alpha 1 actin (ACTA1), myogenin/myogenic factor 4 (MYOG), embryonic skeletal muscle myosin heavy chain 3 (MYH3) and the activity of creatine phosphokinase (CK). The geNorm, NormFinder and BestKeeper software programs were used to ascertain the most suitable RG to normalize the RNA input. Results Using the geNorm program, RPLPO and TBP were found to be the most stable genes, additionally a suitable normalization factor (NF) was calculated. The NormFinder software showed that RPLPO was the most stable, whereas TBP ranked second. BestKeeper program also revealed that RPLPO and TBP as stable genes, but PPIA was the most stable reference gene, whereas GAPDH and ACTB were the worst ranked. Conclusion RNA expression analyses including three independent softwares revealed that RPLPO, TBP as reference genes or NF calculated by geNorm software, are suitable to normalize the mRNA expression in myoblast after culture under differentiation conditions. Significant correlations can be identified between the differentiations markers ACTA1, MYOG, MYH3 and creatine phosphokinase (CK) activity, when the expression is normalized with the NF calculated with RPLPO and TBP.


Background
The real-time polymerase chain reaction (qRT-PCR) has revolutionized the field of gene expression analysis in living organisms. In comparison to classical semi-quantitative reverse transcription-PCR (sqRT-PCR), the main advantages of qRT-PCR are its higher sensitivity, specificity, and broad quantification range [1,2]. Despite being an extremely powerful technique, qRT-PCR suffers from certain pitfalls, the most important being the normalization with a reliable reference gene, habitually called 'housekeeping gene (HKG)'. The housekeeping term was initially given to genes that are necessary for cell function and being constitutively expressed in each cell type [3] but obviously it would be more precise to call it 'reference gene (RG)'. RG are often taken from the literature and used across a variety of experimental conditions, some of which may induce differences in the RG own expression under certain conditions [3]. Thus, experimental results are highly dependent on the RG gene chosen [4]. If unrecognized, unexpected changes in RG expression could result in erroneous conclusions about real biological effects such as responses to drugs [4][5][6]. Such errors could be introduced at a number of stages throughout the experimental protocol (input sample, RNA extraction, reverse transcription, etc.) [1,2,7]. Therefore, it is critical to correct any errors between samples when measuring RNA expression. During the selection of an adequate RG, several indications should be taken into consideration: 1) the expression of the RG must remain constant throughout the intervention; 2) the amplification efficiency of the RG should be similar to that of the genes of interest (GOI); and 3) the abundance of the RG should be similar to that of the genes of interest [8]. Therefore, appropriate RG validation of internal references is crucial in order to avoid misinterpretations of study findings [9]. If the selected RG fluctuate randomly among samples, then subtle differences between GOI will be lost. For example, variation in presumably stable RG was shown in studies examining their expression in serum-stimulated fibroblasts [9,10].
Skeletal muscle tissue engineering aims at the reconstruction of skeletal muscle loss. The field of skeletal muscle tissue engineering has involved to create efficient muscle tissue by using the regenerative potential of stem cells and their potential for proliferation and maturation [11]. The preferred sources of cells for skeletal muscle tissue engineering applications are primary satellite cells [12]. Human satellite cells can be successfully extracted and expanded in vitro and differentiated into myofibers using different stimuli. In this context, previous reports have used conventional RG [8], but there is no data available about the validity of these RG during serum-dependent differentiation of human myoblasts in vitro. Our main GOI were differentiation markers that have unknown RNA expression under experimental conditions. Conse-quently, we could anticipate that differences between the treatments might be small. Therefore, it was highly important to find an internal reference with minimal variability. Thus, we used qRT-PCR to investigate the levels of six RG expressed in cultured myoblasts. Here we report the most suitable RG to detect small variations during the differentiation of myoblasts under certain culture conditions.

Differentiation of Myoblasts
Samples indicated various stages of differentiation. We initially verified that our samples truly represented differentiating cells by analyzing the creatine phosphokinase (CK) activity level, a well-known and accepted marker of differentiation [11].

RNA Quality
All RNA samples were examined as to their concentration, purity, and integrity. RNA purity was measured using the NanoDrop ® Spectrophotometer. Based on the absorbance ratio at 260 nm/280 nm (mean ± standard deviation [SD], 1.90 ± 0.05), all RNA samples were pure and protein free. The mean (± SEM) A260/230 ratio of 1.90 ± 0.20 (range from 1.65 -2.10) indicated that the RNA was free of phenol and ethanol. RNA integrity was assessed by the calculation of RIN values using the Agilent 2100 Bioanalyzer. Differentiated and undifferentiated myoblast samples revealed RIN values between 6.0 and 8.0.

Stability of RG within Sample Groups
First, the efficiencies of every run were calculated using standard curves. The mean of the efficiencies of genes of interest (GOI) and RG are similar in a range between 1.94 and 2.00 ( Table 1).
The stability of the six RG was assessed for each sample group. Gene expression levels were measured by qRT-PCR and the expression stabilities were evaluated by the three most commonly used software-based methods: 1-geNorm [13], 2-NormFinder [14] and 3-BestKeeper [15]. Figure 1 shows the average expression stability (M). M values indicate the stability of a given gene's expression. Based on the M value, the majority of RG are below the GeNorm's arbitrary cut-off level of 0.5 for stability ( Figure  1) automatically selected by the software, suggests that the use of any of these RG for normalization is applicable. Successive elimination of less stable genes according to lowest M values generated a ranking of genes and resulted in the identification of the RPLPO and TBP as the two most stable genes. In addition to the gene stability measure M, the geNorm program calculates a normalization factor (NF). The NF is determined from RPLPO and TBP genes, taking into account the variable as the pair wise var-iation between two sequential normalization factors (data not shown). Table 2 shows the ranking order of the six candidate reference genes mentioned above, using the NormFinder program to calculate their expression stability. Again, RPLPO and TBP were found to be the most stable genes. The NormFinder selected RPLPO as the most stable reference gene with a stability value of 0.093 (pooled data from myoblasts cultured in normal [N] or differentiation [D] medium) and 0.071 when only data of myoblasts cultured in medium N were analyzed (Table 1). RPLPO and TBP were identified as the best genes from myoblasts cultured in medium D with a stability value of 0.050 (Table  2).

3-BestKeeper (BK)
The results of reference gene evaluation by the BestKeeper tool are shown in Table 3. According to the variability observed, RG can be identified as the most stable genes, as they exhibited the lowest coefficient of variance (CV ± SD). In this context, we found RPLPO, TBP or PPIA, with CV ± SD of 4.5 ± 0.83, 4.0 ± 0.99 or 2.9 ± 0.7, respectively, to be the most stable reference genes (Table 3). B2M, ACTB and GAPDH exhibited the highest coefficient of variance with 5.8 ± 1.05, 6.0 ± 1.16 or 7.5 ± 1.12, indicating that these were the least stable NormFinder RG (Table 3, bold italic letters). However, genes that show a SD higher than 1 (= starting template variation by the factor 2) should be considered unacceptable [15]. A low SD of the crossing point (CP) values should be expected for useful RG. Corresponding to the estimation of the SD (± CP) of the CV [% CP], value was highest for ACTB, GAPDH and B2M (Table 3, bold italic letters). This constitutes a reason to exclude these genes from the BestKeeper index calculation, as they are not reliable RG in this setting [15]. The weighted index BestKeeper calculated for the six candidates showed a SD of CP = ± 0.98 cycles. After the exclusion of GAPDH, B2M, and ACTB from the index, its variation decreased (SD = ± 0.83 cycles) to a figure identical to RPLPO (Table 3). GOI expression data are statistically processed with the BestKeeper software in the same way as those of RG, e.g., their GM, AM, SD, CV, Min. and Max. values (Table 4). According to the variability observed, GOI can be identified as well as RG by choosing the most stable genes with the lowest variation, but none of the GOI tested (ACTA1, MYOG and MYH3) showed an acceptable CV ± SD variability (13.77 ± 3.06, 11.62 ± 2.88 and 9.11 ± 2.75 respectively, bold italic letters). A new version of the BestKeeper tool that also employs non-parametric methods is being prepared [15] by which genes with very different expression levels can be compared. On the other hand, all six RG correlated very well one with another (results not shown).

Evaluation of Selected Candidate RG and Normalization Approach
Following identification of the most stable RG and the normalization factor (NF) from the full panel of RG, a method was needed for their evaluation. Correlations among the status of the myoblast differentiation (identified by CK activity) and expression of RG and ACTA1, MYOG and MYH3 (normalized with NF [RPLPO+TBP]) are shown in Table 5. Therefore, we expected these observations to be mirrored at the mRNA level. More specifically, following treatment with the differentiation medium, a high and significant correlation between the CK activity and the expression of ACTA1, MYOG or MYH3 was found (Table 5), whereas RPLPO, TBP and GAPDH showed no significant correlation with CK (Table 5, bold italic letters). However, both positive (B2M) and negative (ACTB) correlation between gene expression and CK were also found ( Table 5). A significant positive correlation was found between the GOI (Table 5), whereas there was no correlation between RPLPO, TBP or GAPDH and ACTA1, MYOG or MYH3. However, B2M showed a good correlation (Table 5). ACTB has shown a negative correlation with ACTA1 and MYH3 (Table 5) but not with MYOG ( Table 5, bold italic letters). Figure 2 shows the percent of GOI variation normalized with NF and the best and worst RG candidates, after 12 days myoblasts differentiation. When the GOI were normalized with NF, RPLPO and TBP, we identify approximately 67.0% of variation; in contrast only 30.0% variation was observed, when MYH3 was normalized with ACTB or GAPDH.

Discussion
Myogenesis of skeletal muscle is a highly complex phenomenon controlled and regulated by multiple intracellular signalling mechanisms [16]. In order to investigate optimal culture conditions and characterization of human myoblasts, we analyzed cell behaviour under different culture media and showed that cultures treated with differentiation medium expressed several characteristic features of mature skeletal muscle [11]. Furthermore, in myoblasts cultured with growth medium, certain markers of differentiation were detected, demonstrating that, for an accurate characterization, multiple markers (e.g. genes, enzymes and structure proteins) must be analyzed [11]. During skeletal muscle differentiation, mononucleated proliferating myoblasts stop dividing, activate musclespecific gene expression, and fuse into multinucleated myotubes [11]; these events can produce drastic or small variation in the expression of genes and proteins. Markers of differentiation were also detected in cultures treated with non-differentiation medium, but there was no formation of myotubes. In the enzymatic assay of CK, a wellknown and accepted marker of differentiation, cultures treated with differentiation medium showed a higher activity, evidencing a higher degree of differentiation [11]. In this context, the strategy used to normalize the mRNA with adequate RG is critical. To identify the best RG as normalisers in gene expression studies, several strategies including computer programs have been recommended [5,6,9,[14][15][16][17][18][19][20]. The absence of the differential expression or variability of the candidate RG examined is the strongest proof of suitability and stability [20]. Therefore, we proposed to use several software programs to analyze six RG and select the RG or normalization factors to correct the RNA input. The major result of our study -evaluating six candidate RG in myoblasts under normal and differentiation conditions, during 8 and 12 days of treatmentwas that the genes RPLPO or TBP were identified by the three independent software programs as suitable RG that did not differ in their expression in differentiated and undifferentiated myoblasts. Consequently, only RPLPO, identified by geNorm, NormFinder and BestKeeper software programs fulfils the criterion of expression stability between samples and can be recommended as an accurate normalizer for relative gene quantification in cultured myoblasts. First, the expression of the candidate RG between the respective conditions was compared using an adequate test such as geNorm, NormFinder and Best-Keeper. It can be assumed that genes with significantly different expressions are not suited for target gene normalization, since they are affected by the study condition of interest. Thus, those genes should be excluded as normalizers. In the present study, the expression of all genes in the differentiated and undifferentiated myoblast samples showed different grades of variability but the most stable were RPLPO and TBP. We used the BestKeeper software to confirm the results obtained by geNorm and NormFinder. The raw CPs analyzed by BestKeeper seem to be good estimators of the expression levels as they are (in most cases) normally distributed and a parametric test can thus be performed [15]. This also gives the CP datasets the Gaussian distribution, justifying usage of parametric methods [15]. Low numbers of expressed genes where CPs were obtained somewhere surely show different variances compared to highly expressed genes with CPs. These higher levels of CV and SD invalidate the use of Best-Keeper as a reference to normalizing the expression of GOI [15]. The correlation analysis of the normalized gene with NF showed a bad correlation between the RG RPLPO, TBP or GAPDH and the different differentiation markers (CK, ACTA1, MYOG and MYH3). In this context, the absence of correlation confirmed that the RG selected  (14). The analysis of the expression stability was performed taken the expression of the RG from RNA samples of myoblasts. The analysis was performed using myoblasts cultured in normal medium (N), or differentiation medium (D). Additionally, the analysis of the combination of N and D together is shown. *High expression stability is indicated by a low stability value as an estimate of the combined intra-and intergroup variation of the individual gene. SE stands for standard error. Representative results from three experiments.  by geNorm, NormFinder and BestKeeper are suitable as reference for the normalization. In addition, we found a significant correlation between B2M or ACTB and the GOI which indicates a variation associated with the differentiation status of the myoblast. Thus, they cannot be used as reference genes.
Two or three genes represent a realistic calculation basis in a common laboratory and the minimum necessary number for a good analysis [15]. RPLPO STD dev CP [± SD] is identical to BestKeeper calculated with RPLPO, TBP and PPIA. In this context, the RG, RPLPO and TBP, detected by BestKeeper agree with the best RG detected by geNorm and NormFinder.
The normalization with NF, RPLPO or TBP allows measurements of differences in the expression of the differentiation gene-markers between myoblasts treated with differentiation medium or normal medium (control).

Conclusion
In conclusion, the RG RPLPO, TBP as well as a NF calculated by geNorm software are recommended as references for relative gene quantification in gene profiling studies of single genes in myoblasts under serum-differentiation conditions. The confirmation of RG by the three software programs is a suitable method to perform an adequate normalization of the mRNA input in quantitative PCR experiments.

Cell Culture
Human skeletal muscle biopsies were obtained from 15 patients during head and neck surgery. The median age of these patients was 58, ranging from 41 to 72. The study was approved by the Ethics Committee of the Mannheim Faculty of Medicine, University of Heidelberg, Germany, and the patients gave informed consent. Satellite cells were dissociated from the minced muscles by digestion with collagenase B (Roche, Mannheim, Germany) for 60 min and 0.05% trypsin-0.02% EDTA (Promo Cell, Heidelberg, Germany) for 45 min at 37°C, filtered through a sterile 70-μm cell strainer (Becton Dickinson, Franklin Lakes, NJ, USA) and purified with the pre-plating technique as recently described [12]. Purity of myoblast cultures (>80%) was evaluated by anti-desmin immunostainings. Cells were grown on 0.2% gelatinecoated culture flasks (Sigma, Deisenhofen, Germany) and in Ham's F10 growth medium, containing 1% penicillin/ streptomycin/fungizone-solution (PSF), 2 mM Lglutamine (all from PromoCell, Heidelberg, Germany) and 10% foetal bovine serum (FBS) (PAA Laboratories, Linz, Austria) (GM). Cells were maintained at 37°C in a humidified atmosphere of 5% CO2 and 95% air, and the medium was changed every 72 h.

Differentiation of skeletal muscle myoblasts
Skeletal myoblasts were cultured to ~60% confluence, washed with phosphate-buffered saline (PBS) and induced to differentiate by a change in medium (D) consisting of minimal essential medium (MEM) (PromoCell) supplemented with 2% horse serum (PAA Laboratories), 2 mM L-glutamine and PSF (both from PromoCell).

RNA Isolation
Total RNA was isolated with the RNA Mini Kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. RNA concentration, purity, and integrity was determined by A260 and A280 (A260/A280 = 1.7-2.0) measurements using a NanoDrop 8000 Spectrophotometer (Thermo Scientific, Schwerte, Germany) and Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany).  . Amplification was performed using the Mx3005P™ QPCR System (Stratagene). For relative quantification, a standard curve was generated in every individual run. Shortly, total RNA was pooled from muscle biopsies of healthy human volunteers, reverse transcription was performed and a serial dilution of the cDNA was used to perform the calibration curve. The data were analyzed using the relative standard curve method. For each unknown sample, the relative amount is calculated using linear regression analysis from their respective standard curves. Data were analyzed using the Mx3005P analysis software (Stratagene-Agilent Technologies, Waldbronn, Germany). The efficiencies of all GOI and RG were calculated in every individual run (Table 1).

Creatine Phosphokinase (CK) Analysis
Differentiation of myoblasts was identified by measuring the creatine phosphokinase (CK) activity of developing myofibers using the CK assay (Sigma-Aldrich Chemie, Taufkirchen, Germany). Cells were lysed with a buffer containing NP-40, protease inhibitors and were stored at -80°C. Samples were assayed according to the manufacturer's protocol. Enzymatic activity was normalized over the total protein content determined by the Lowry protein assay (Bio-Rad, Munich, Germany).

Analysis of Expression Stability
For stability comparisons of candidate RG within sample groups, the software geNorm, version 3.4 [13] (Visual Basic application tool for Microsoft Excel), NormFinder [14] and BestKeeper [15], were used according to developer's recommendations.
GeNorm uses a gene-stability measure M, which is defined as the average pair wise variation between a particular gene and all other control genes. It calculates the optimal number of genes necessary for normalization of a target gene and combines them into a normalization factor (NF). The NormFinder program uses a model-based approach for estimation of expression and enables the identification of the single best genes as well as giving a ranking order. The BestKeeper software determines the best suited reference genes and combines them into an index (BK). The index can be compared with further target genes to decide whether they are differentially expressed under certain conditions. In addition, the option in NormFinder to define groups was applied to compare the effect of the treatment of myoblasts with differentiation medium.

Statistical analysis
The software Sigma Plot was used to carry out statistical analysis, by paired Student's t-test or Mann-Whitney-U test, as well as Pearson's product moment correlation test of gene expression among experimental groups. Gene expression results were expressed as the mean. P < 0.05 was considered statistically significant.