Reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) remains one of the best-established techniques to assess gene expression patterns. However, appropriate reference gene(s) selection remains a critical and challenging subject in which inappropriate reference gene selction can distort results leading to false interpretations. To date, mixed opinions still exist in how to choose the most optimal reference gene sets in accodrance to the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guideline. Therefore, the purpose of this study was to investigate which schemes were the most feasible for the identification of reference genes in a bone and cartilage bioengineering experimental setting. In this study, rat bone mesenchymal stem cells (rBMSCs), skeletal muscle tissue and adipose tissue were utilized, undergoing either chondrogenic or osteogenic induction, to investigate the optimal reference gene set identification scheme that would subsequently ensure stable and accurate interpretation of gene expression in bone and cartilage bioengineering.
The stability and pairwise variance of eight candidate reference genes were analyzed using geNorm. The V0.15- vs. Vmin-based normalization scheme in rBMSCs had no significant effect on the eventual normalization of target genes. In terms of the muscle tissue, the results of the correlation of NF values between the V0.15 and Vmin schemes and the variance of target genes expression levels generated by these two schemes showed that different schemes do indeed have a significant effect on the eventual normalization of target genes. Three selection schemes were adopted in terms of the adipose tissue, including the three optimal reference genes (Opt3), V0.20 and Vmin schemes, and the analysis of NF values with eventual normalization of target genes showed that the different selection schemes also have a significant effect on the eventual normalization of target genes.
Based on these results, the proposed cut-off value of Vn/n + 1 under 0.15, according to the geNorm algorithm, should be considered with caution. For cell only experiments, at least rBMSCs, a Vn/n + 1 under 0.15 is sufficient in RT-qPCR studies. However, when using certain tissue types such as skeletal muscle and adipose tissue the minimum Vn/n + 1 should be used instead as this provides a far superior mode of generating accurate gene expression results. We thus recommended that when the stability and variation of a candidate reference genes in a specific study is unclear the minimum Vn/n + 1 should always be used as this ensures the best and most accurate gene expression value is achieved during RT-qPCR assays.
The successful re-formation of bone and cartilage remains an unsolved riddle to achieve clinically. Whilst many new bioengineering concepts have shown great potential to possibly someday replace the golden standard utilizing autografts, bone or cartilage [1,2,3], their translatability clinically remain problematic. To interpret and manipulate the nature’s biological defined processes successful for clinical applications, detailed and accurate assays in genes expression patterns and modulatory mechanism are a critical fundamental knowledge that if not properly defined will continue to generate more questions than answers. Reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) has emerged as one of the fundamental techniques to properly solve the questions sought for bioengineering principles . If utilized properly, it has shown to generate reliable, comparable and reproducible gene expression information on how tissues respond during bioengineering  making RT-qPCR a benchmark for gene analysis [4, 6] but also a critical validation tool to support NextGen sequencing and microarray assay results [7, 8]. Nevertheless, improper optimization and standardization have shown to significantly affect the variability of gene expression results generated by RT-qPCR thereby impairing the reproducibility that subsequently compromises the translation efficiency of present bioengineering techniques [4, 9,10,11,12].
To prevent the misinterpretation of results caused by non-standard experimental procedures and details , the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines was established, which targets the promotion of reliability of results including the consistency and transparency between laboratories [9, 14,15,16]. Particularly, the selection of reference genes has drawn considerable attention, because the expressions of these so-called “ideal” non-variant genes can become unstable under certain conditions [15, 17, 18]. Hence, the proper selection and combination of multiple reference genes was established to minimize the instability and variation . Additionally, various methodology articles have made further efforts to ensure the stability and optimal quantity of reference genes to obtain accurately reproducible data, providing considerable impetus towards perfecting RT-qPCR [19,20,21,22]. Subsequent bone and cartilage bioengineering studies have suggested that the stability and normalization quantity of reference genes should not only be determined by cell or tissue type but should also be re-optimized for experiments under different processing conditions [23, 24], which has provided an extra step in perfecting the understanding of the MIQE guidelines.
GeNorm, a bioinformatics tool, is commonly used to find the most stable reference genes and determine the proper quantity by calculating the M value, normalization factor (NF) and the pairwise variation (Vn/n + 1-score) . The value of Vn/n + 1 with 0.15 (V0.15) was generally accepted as the cut-off for choosing optimal number of reference genes, below which the participation of more reference genes was thought redundant . Besides the Vn/n + 1-score-based scheme, the geNorm algorithm also provides an alternative to select the best three reference genes based on geNorm results. Moreover, other schemes were raised due to the limitation of the V0.15 method [25,26,27]. In the study performed by Ayers et al.  and Hosseini et al. , 0.20 was set as a trade-off where in some tissue types (e.g. adipose tissue) the minimum Vn/n + 1-score (Vmin) was higher than 0.15 [18, 19]. To select stable reference genes for a cell-based study, Lee et al.  assessed the stability of twelve candidate genes across experimental conditions by geNorm analysis, finally abandoned the strategy of selecting the three most stably expressed reference genes and instead adopted the Vmin strategy. Hence, debates persist on the schemes for choosing reference genes to improve the accuracy of RT-qPCR assays [25, 27,28,29,30].
Therefore, this study sought to investigate which schemes (V0.15-, V0.20-, VOpt3- and Vmin-based normalization schemes) were the most feasible for the identification of reference genes in bone and cartilage bioengineering experiments, as this is our primary research direction. Our hypothesis was that identification of reference genes based on Vmin was the optimal scheme for the normalization of RT-qPCR, thereby gaining accurate and reliable gene expression data for bone and cartilage bioengineering.
V0.15- vs. Vmin-based normalization scheme for analyzing gene expression data in rBMSCs
The correlation of NF values between the V0.15 and Vmin schemes was analyzed and the variance of target genes expression levels generated by these two schemes was compared in rat bone mesenchymal stem cells (rBMSCs).
In the osteogenic sub-study, the V0.15 was 0.078, while the Vmin was 0.055 (Fig. 1a). Hence, combining the sequencing of eight candidate reference genes based on M-value, the V0.15-based reference gene set contained ribosomal protein L13α (Rpl13α) and actin beta (Actb), while the Vmin-based reference gene set contained Rpl13α, Actb and RNA polymerase II subunit e (Polr2e) (Fig. 1b). The r- value in the Spearman rank correlation analysis between NFV0.15 and NFVmin was 0.9762 (Fig. 1c), which showed a very strong correlation between these two schemes in terms of rBMSCs. The variance of target genes expression levels using these two schemes was then compared. The osteogenic-related target genes, including bone morphogenetic protein-2 (Bmp-2), Bmp-6, osteocalcin (Ocn) and runt-related transcription factor 2 (Runx2) were normalized to the two reference gene sets, and the calibrated normalized relative quantity (CNRQ) values were obtained. The CNRQ values of all the investigated target genes, either using the V0.15 or Vmin scheme, did not show a significant difference (P > 0.05) (Fig. 1d), revealing that the two reference genes identification schemes had no significant effect on the eventual normalization of target genes.
In the chondrogenic sub-study, V0.15 was 0.068, while Vmin was 0.058 (Fig. 2a). Hence, combining the sequences of eight candidate reference genes based on the M-value, the V0.15-based reference gene set contained Actb and Rpl13α, and The Vmin-based reference gene set contained Rpl13α, Actb and Polr2e(Fig. 2b). The r-value in the Spearman rank correlation analysis between NFV0.15 and NFVmin was 0.9524 (Fig. 2c). These results showed that there was a strong correlation between these two selection schemes. The chondrogenic-related target genes, including aggrecan (Acan), Sex determining region Y-box 9 (Sox9), transforming growth factor beta 1 (Tgf-β1) and Tgf-β3, were normalized to the two reference gene sets. CNRQ values showed that the relative expression levels of Acan, Sox9, Tgf-β1 and Tgf-β3 did not show a significant deference between the V0.15- and Vmin- based CNRQ values (P > 0.05), which revealed that different selection schemes had no significant effect on the eventual normalization of target genes (Fig. 2d).
V0.15- vs. Vmin-based normalization scheme for analyzing gene expression data in skeletal muscle tissue
The correlation of NF values between the V0.15 and Vmin schemes was analyzed and the variance of target genes expression levels generated by these two schemes was compared in terms of the muscle tissue.
In the osteogenic sub-study, the V0.15 was 0.108, while Vmin was 0.067 (Fig. 3a). Hence, combining the sequencing of eight candidate reference genes based on M-value, the V0.15-based reference gene set contained Polr2e and Rpl13α, while the Vmin-based reference gene set contained Rpl13α, Polr2e, glyceraldehyde 3-phosphate dehydrogenase (Gapdh), TATA-binding protein (Tbp), Actb and ribosomal protein lateral stalk subunit P0 (Rplp0) (Fig. 3b). The r-value in the Spearman rank correlation analysis between NFV0.15 and NFVmin was 0.6429 (Fig. 3c). This result showed that there was no strong correlation between these two selection schemes. The osteogenic-related target genes, including Bmp-2, Bmp-6, Ocn and Runx2, were normalized to the two reference gene sets. The relative expression levels of Ocn and Runx2 showed no significant deference between the V0.15- and Vmin- based CNRQ values (P > 0.05), while the results of Bmp-2 and Bmp-6 showed significant deference (P < 0.05), which indicated that different schemes may have a significant effect on the eventual normalization of specific target genes (Fig. 3d).
In the chondrogenic sub-study, the V0.15 was 0.079, while the Vmin was 0.057 (Fig. 4a). Hence, combining the sequencing of eight candidate reference genes based on the M-value, the V0.15-based reference gene set contained Gapdh and Rplp0, and the Vmin-based reference gene set contained Rplp0, Gapdh, Actb, Polr2e and Tbp(Fig. 4b). The correlation coefficient in the Spearman rank correlation analysis between NFV0.15 and NFVmin was 0.833. This result showed that there was a rare correlation between these two selection schemes (Fig. 4c). The chondrogenic-related target genes, including Acan, Sox9, Tgf-β1 and Tgf-β3, were normalized to the two reference gene sets. The relative expression levels of Acan, Sox9, Tgf-β1 and Tgf-β3 showed a significant deference between the V0.15- and Vmin- based CNRQ values (P < 0.05), revealing that here the different selection schemes have a significant effect on the eventual normalization of target genes (Fig. 4d).
V0.20/Opt3 vs. Vmin-based normalization scheme for analyzing gene expression data in adipose tissue
Three selection schemes were adopted in terms of the adipose tissue, either in the osteogenic or the chondrogenic sub-study including optimally three reference genes (Opt3), V0.20 and Vmin schemes. The correlation of NF values between the V0.20 and Vmin schemes with Opt3 and Vmin schemes were analyzed including the variance of target genes expression levels generated by these three schemes was compared.
In the osteogenic adipose sub-study, V0.20 was 0.194, while the Vmin was 0.170 (Fig. 5a). Hence, combining the sequencing of eight candidate reference genes based on the M-value, the Opt3-based reference gene set contained RNA 28S ribosomal 4 (Rna28s4), Gapdh and Actb, the V0.20-based reference gene set contained Rna28s4, Gapdh, Actb and Rpl13α, while the Vmin-based reference gene set contained Rna28s4, Gapdh, Actb, Rpl13α, Polr2e and Tbp(Fig. 5b). The r-value in the Spearman rank correlation analysis between NFopt3 and NFVmin was 0.833, with that of NFV0.20 and NFVmin was 0.762 (Fig. 5c, e). These results showed that there was neither a correlation between Vmin and V0.20 schemes nor between Vmin and Opt3 schemes. The osteogenic-related target genes, including Bmp-2, Bmp-6, Ocn and Runx2, were normalized by three reference gene sets. The relative expression levels of Bmp-2, Bmp-6, Ocn and Runx2 all showed significant deference between the Opt3- and Vmin- based CNRQ values including between the V0.20- and Vmin- based CNRQ values (P < 0.05), which revealed that different selection schemes have significant effects on the eventual normalization of target genes (Fig. 5d, f).
In the chondrogenic adipose sub-study, V0.20 was 0.181, while the Vmin was 0.168 (Fig. 6a). Hence, combining the sequencing of eight candidate reference genes based on M value, the Opt3-based reference gene set contained Rplp0, Gapdh and Polr2e, the V0.20-based reference gene set contained Rplp0, Gapdh, Polr2e and Rpl13α, while the Vmin-based reference gene set contained Rplp0, Gapdh, Polr2e, Rpl13α, Rna28s4, Tbp and succinate dehydrogenase complex flavoprotein subunit A (Sdha)(Fig. 6b). The r-value in the Spearman rank correlation analysis between NFopt3 and NFVmin was 0.4286, whilst between NFV0.20 and NFVmin was 0.6429 (Fig. 6c, e). These results showed that there was neither a correlation between Vmin and V0.20 schemes nor between Vmin and Opt3 schemes. The chondrogenic-related target genes, including Acan, Sox9, Tgf-β1 and Tgf-β3, were normalized by three reference gene sets, and the CNRQ values were obtained. The relative expression levels of Acan, Sox9, Tgf-β1 and Tgf-β3 showed significant deference between the Vopt3- and Vmin- based CNRQ values or between the V0.20- and Vmin- based CNRQ values (P < 0.05), which revealed that different selection schemes have significant effect on the eventual normalization of target genes (Fig. 6d, f).
Accurate gene analysis remains essential for objective evaluation of the efficacy of bioengineering [31, 32]. Although the establishment of MIQE Guidelines and the advocacy of multiple reference genes has increased the reliability of the results [9, 15], there are still no unified solutions for the choice of the quantity and priority of reference genes, making much of the results across the bioengineering spectrum inconsistent [33, 34]. Several mathematical algorithms have been developed to solve this dilemma intended at detecting the stability of reference genes [14, 35, 36], such as geNorm , NormFinder , and BestKeeper . In the present study, we as such sought to optimize the choice in geNorm schemes for reference genes determination where normally 0.15, commonly regarded as the cut-off point in the pairwise variation analysis, is a dangerous generalization as mono-cell cultures versus tissue-based experiments show clear deviations when using Vmin or V0.15.
By integrating the results of geNorm, NormFinder, and BestKeeper, Tong et al.  and Zhang et al.  suggested that they had identified suitable reference gene(s) for normalization in gene expression with regards to their own specific experiments. It seems that the combination of different algorithms is a popular method [12, 14, 37, 38]. However, in our previous study, we indicated that the algorithms utilized were irrelevant with regard to the stability of reference genes evaluation; in other words, the different programs may not affect which reference genes should be selected and the optimum number needed to generate accurate results . GeNorm remains one of the most accessible platforms in which a user-friendly interface, intuitive data presentation and subsequent infinite sample inclusion including candidate reference genes, make highly popular compared to the other algorithms . Whilst one of the limitations of geNorm considers that co-related genes could occupy adjacent positions in the ranking, possibly causing a selection bias in the pairwise measurements [23, 38], our results keep validating the opposite in which the consistent stability of reference genes was maintained when these were removed successively. This indicates there is no co-regulation in the genes selected (data not shown). It was of critical importance to apply the mathematical algorithm accurately, which was the same purpose of this study and geNorm was considered a reliable and convenient method for analyzing the stability and determining the optimal number of the reference genes in a specific research endeavor [15, 23].
In terms of geNorm, the M-value reflects the stability of each reference gene, where the lowest value embodies the most stable expression, while the determination of the optimal number of reference genes for accurate gene expression analysis relies on the pattern of pairwise variation, the V-score . The usual V-curve showed that adding reference genes is a double-edged sword in the normalization process. Whilst the nonspecific variation was eliminated by adding the stable reference genes and proven by the decrease in V-score, the increase after the minimum V-score indicated that unstable reference genes could interfere with the normalization process. Hence, it was reasonable that the participation of more reference genes after the minimum V-score do not contribute to minimizing the instability and variation of the reference gene set. In the initial design, Spearman rank correlation analysis was used to analyze the correlation of NF and then showed the pairwise variation (Vn/n + 1) to assist determining the number of reference genes required for accurate normalization, in which the high correlation coefficients corresponded to the low V-scores . Vandesompele et al.  decided to take 0.15 as the cut-off, below which the high correlation suggested that it was not necessary to include more reference genes. Indeed, our results of rBMSCs in the present study showed that the high correlation of the NFs obtained from two different gene sets, determined by the V0.15 and Vmin (r = 0.95 for chondrogenesis; 0.97 for osteogenesis), respectively, confirmed the previous suggestion. However, the analysis of osteogenic induced muscle fragments presented a relatively low correlation coefficient when comparing the NFV0.15 with NFVmin (r = 0.83 for chondrogenesis; 0.64 for osteogenesis), which suggested that the variation caused by the additional reference genes only contained in the scheme based on Vmin was significant. Erkens et al.  and Pérez et al.  all attempted to evaluate the Vmin-based scheme and they acknowledged that normalization became more accurate, but the improvement was minimal. Unfortunately, no corresponding verifications were performed in their reports, such as Spearman rank correlation analysis. Ragni et al.  also noticed the differences caused by diverse reference gene sets, but the analysis method used was to observe whether significant changes in the target genes occurred when different combinations of reference genes were applied for normalization, which supplied our experiments with a new evaluation method. The relative expression levels of chondrogenic- and part of osteogenic-related genes in our muscle tissue study showed significant differences when normalized by two reference gene sets (Vmin- or V0.15-based), but not in the milieu of all target gene when using rBMSCs. The results revealed that the variation erased by Vmin based sets cannot be neglected in muscle tissue, which was consistent with the outcomes in Spearman rank correlation analysis. On the basis of the current data, regarding rBMSCs, the selection of reference genes based on V0.15 is fully viable. However, when using muscle tissue, the conditions need be stricter, where we recommend using the Vmin scheme to obtain a more optimized set for normalization.
Interestingly, the pairwise variation analysis in the adipose tissue showed that none of the V-score was lower than 0.15, causing the V0.15 based scheme inapplicable. Facing a similar dilemma, Ayers et al.  flexibly elevated the threshold to 0.2, coincidentally, where the V-score corresponding to the decided number of reference genes was the minimum one in their analysis results. In a bovine adipose tissue explants study, Hosseini et al.  also drew on this new cut-off (V0.2) when determining the optimal reference genes due to the lack of values lower than 0.15. In the subsequent Spearman rank correlation analysis, although involving more reference genes provided a more stable normalization (r = 0.93), compared with using less reference genes based on the cut-off of 0.2 (r = 0.85), Hosseini et al. assumed it as “marginal” and hence ignored this difference . In our study, the correlation between two NFs which corresponded to the Vmin and V0.2 schemes were also weak, suggesting that Vmin, rather than a fixed threshold (e.g. 0.2), was still a reliable choice when all V-scores were higher than 0.15. Additionally, only choosing the three best reference genes was a compromise recommended by Vandesompele et al.  and Lu et al. , which we consider improper as according to the correlation analysis or the comparison of the relative expression levels of target genes in our adipose-related study the accuracy of normalization can be significantly improved when using Vmin based reference genes, highlighting the advantages of the minimum scheme once again.
It is worth considering why different reference gene sets have significantly diverse performance in tissue fragments while keeping stable in rBMSCs, which may rely on the inherent difference between cells and tissue. Research conducted by Vandesompele et al.  also showed that the participation of more reference genes was required to remove non-specific variations in the normalization of genes in a tissue study. Compared with specific cell culture systems tissue models are more complex in which the heterogeneity of various cell types generates a multifaceted response when exposed to a single stimulus [42, 43]. Studies have shown that for different cell types, the expression of reference genes was not as stable [23, 44], not to mention when multiple reference genes and numerous cell types required to be weighed simultaneously, the change of a single factor, such as the number of reference genes, probably led to significant difference. Ren et al.  demonstrated that the response to osteogenic stimulation in muscle fragment with or without fascia was significantly different, suggesting that endogenous gene expression induced by the integration of external stimuli would change with the increase of cell types or tissue structures. Furthermore, concerning the RNA extraction, it was difficult for cells in tissues to be isolated independently from the extracellular matrix, while the influence of mRNA and protein in the matrix during the purification resulted in the enhancement of differences . Homogenization was another potential threat, which is essential for RNA extraction from tissue models. However, the reduction of total RNA due to inadequate grinding, the limitation of the number of samples per homogenization and local temperature changes due to the high shearing forces could lead to loss, inactivation and degeneration of unstable cellular components .
From the present study, we recommend that the proposed 0.15 cut-off value, according to the geNorm algorithm, should be carefully considered as cell and tissue experiments show clear variations. Whilst our results presently only reflect this for skeletal muscle and adipose tissue, some of the favorite tissues used in bioengineering experiments, whether this is applicable to all tissues needs to be further verified. However, from the results we can with certainty state that if the stability and variation of candidate reference genes, in a specific study are unclear, we recommend that the minimum V-score should be used as this provides a superior selection for determining the optimal number and category of reference genes needed to generate accurate and reproducible gene expression results. If not, we fear that the issue generating superior and critically reproducible gene expression results accurately will remain problematic, continuing to mislead the bioengineering field in developing reliable working applications that clinically are needed to bioengineer lost or damaged tissue types.
In the present study, either osteogenic or chondrogenic induction was applied to three commonly used cell or tissue types in bone and cartilage bioengineering: rBMSCs, rat skeletal muscle and adipose tissue. Eight candidate reference genes and osteogenic or chondrogenic related target genes were examined by RT-qPCR. Subsequently, the stability and pairwise variance of candidate reference genes were analyzed using geNorm . The reference gene sets identified by different schemes would produce different NF values and different normalization results. The correlation of NF values and the variance of CNRQ values were performed to prove the variance between different schemes, then to define the optimal one for reference genes identification (Fig. 7).
Cell and tissue specimens
RBMSCs (passage 0, Sciencell, Carlsbad, CA, USA), skeletal muscle and adipose tissue from four F-344 adult female rats (Charles River Laboratories Wilmington, MA, USA) were used in this study. A total of 12 specimens per cell or tissue type were used with 4 specimens acting as the untreated control and the remaining 8 specimens being treated to either undergo chondrogenic (n = 4) or osteogenic differentiation (n = 4).
Skeletal muscle tissue and adipose tissue harvest
For the tissue part of the study, a single F-344 adult female rat was sacrificed using an overdose of isoflurane (Abbot, Chicago, IL, USA). All practical experimental steps were performed in keeping with the rules and regulations of the Animal Protection Laboratory Animal Regulations (2013) and approved by the Animal Care Committee of Renji Hospital (Shanghai, China, No.: 201606230235). Under sterile conditions, fresh abdominal muscle and subcutaneous adipose tissue was harvested and placed temporarily in Dulbecco’s modified Eagle’s medium (DMEM; Biochrom Ltd., Cambridge, United Kingdom) containing high concentrations of Penicillin/Streptomycin (9%, P/S, Biochrom GmbH). Muscle (n = 8) and adipose (n = 8) tissue fragments were then collected using a 5 mm diameter biopsy punch (PFM medical, Cologne, Germany) and transferred into 24-well Nunc well culturing plates (Thermo Fisher Scientific, Waltham, MA, USA) in recovery medium consisting of DMEM supplemented with 3% P/S for 48 h at 37 °C containing 5% CO2 being treated to either undergo chondrogenic (n = 4) or osteogenic differentiation (n = 4). Fresh muscle and adipose tissue fragments (n = 4) were also collected as these would serve as the endogenous normalization control to which all tissue samples would be compared to.
RBMSCs were used for the cellular culturing part, which were seeded at a density of 2*104 per monolayer flask (Thermo Fisher Scientific) and cultured in DMEM supplemented with 3% P/S at 37 °C containing 5% CO2 for the primary culture. When cells reached 80% of confluence, they were detached using trypsin–EDTA (Biochrom Ltd), washed and submitted to new monolayer flasks at the same density for the sub-culture in the same manner thereafter. Cell morphology was observed under the light microscope, and photographs were taken. Cell numbers were counted at each cell passaged. Once 2nd passage cells reached 80% of confluence, they were passaged for the following induction of differentiation procedures. Some of the rBMSCs in 2nd passage pure without culturing were collected immediately as these would be used as the endogenous normalization control in downstream analysis procedures.
Chondrogenic and osteogenic differentiation
To stimulate chondrogenic or osteogenic differentiation in both tissue and cell types the relevant media were utilized. The chondrogenic differentiation medium consisted of normal growth medium supplemented with 10 ng/mL recombinant BMP-6 (R&D Systems, Minneapolis, MN, USA), 10 ng/mL recombinant TGF-β3 (R&D Systems), 100 nM dexamethasone (Sigma-Aldrich, St. Louis, MI, USA), 50 μg/mL L-ascorbic acid-2-phosphate (Sigma-Aldrich), 40 μg/mL L-proline (Sigma-Aldrich), ITS+ 1(10 mg/L insulin, 5.5 mg/L transferrin, 4.7 μg/mL linoleic acid, 0.5 mg/mL bovine serum albumin, and 5 μg/L selenium) (Sigma-Aldrich) [32, 47,48,49]; the osteogenic differentiation medium consisted of normal growth medium supplemented with 50 μg/mL L-ascorbic acid-2-phosphate (Sigma-Aldrich), 1 mM L-glutamine (Sigma-Aldrich), 10 mM β-glycerophosphate (Sigma-Aldrich) and 100 nM dexamethasone (Sigma-Aldrich) ; the normal medium was DMEM supplemented with 3% P/S.
After 48 h recovery, harvested tissue specimens were collected, allocated randomly and then cultured in either chondrogenic differentiation (n = 4 per tissue type) or the osteogenic medium (n = 4 per tissue type), with normal recovery medium (n = 4 per tissue type) acting as the experimental control group. Tissue fragments were cultured for 7 days, medium changed every 2 days, collected and snap frozen in liquid nitrogen and then stored at − 80 °C until further use within 4 weeks.
Similarly, rBMSCs (passage 2) once having reached 80% confluence were trypsinized and seeded at 2*104 cells per culture flask. Chondrogenesis (n = 4) or osteogenesis (n = 4) was then induced by utilizing the corresponding chondrogenic or osteogenic differentiation medium, respectively. Normal medium (n = 4) acted as the experimental control. The medium was changed every 2 days, and 7 days later the cells were harvested, immersed in trizol (Ambion, Carlsbad, CA, USA) and stored at − 80 °C for downstream analysis procedures.
Primer design and optimization
Candidate reference genes were selected out of a gene library pool, known to be suitable for the optimization of reference genes in RT-qPCR, all with a standard deviation of the average amplification threshold cycle quantification value (Cq) less than 1 across 35 in rat tissues [51, 52]. Out of the candidate reference genes pool, the following eight genes were selected as candidates: Rpl13α, Gapdh, Tbp, Rna28s4, Polr2e, Actb, Rplp0, and Sdha. To study mRNA expression of the genes implicated in chondrogenesis, four chondrogenic-related genes were selected including Acan, Sox9, Tgf-β1 and Tgf-β3. Meanwhile, osteogenic-related genes included Bmp-2, Bmp-6, Ocn and Runx2. Primer sequences were designed utilizing PrimeQuest in conjunction with OligoAnalyzer 3.1 (https://eu.idtdna.com/site) and cross-referenced using the Basic Local Alignment Search Tool program (https://blast.ncbi.nlm.nih.gov/Blast.cgi) (see Additional file 1). All the primer sequences with relevant parameters were presented in Table 1.
As previously established , primers were then stringently assessed for sequence amplification specificity with the annealing temperature predetermined to best function at 60 °C. A melt curve was included in each run to confirm amplification of a single product. After PCR amplification wells identified with positive amplicons underwent purified using the Mini Elute PCR Purification Kit (Qiagen, Crawley, UK) and analyzed, after Sanger sequencing (GATC Biotech, Cologne, Germany) utilizing BLASTN (https://blast.ncbi.nlm.nih.gov/Blast.cgi? PROGRAM = blastn&PAGE_TYPE = BlastSearch&.
Cells and muscle tissue samples were homogenized by liquid nitrogen in conjunction with a mortar and pestle. For all adipose tissue samples, a Micro-Dismembrator S (Sartorius Stedim Biotech, Göttingen, Germany) was utilized to homogenize. Subsequently, the RNeasy Fibrous Tissue Mini Kit (Qiagen, Hilden, Germany) was used to extract total RNA following the manufacturer’s protocol and eliminate DNAse and RNAse. Total extracted RNA concentration was determined spectrophotometrically at A260/280 with a NanoDrop™ Lite (Thermo Scientific, Waltham, USA) and RNA quality was assessed with a Pico6000 RNA kit (Agilent Technologies, Santa Clara, CA, USA) on a Bioanalyzer 2100 (Agilent Technologies). RNA integrity numbers were 8 for the cells and 7.8 for the muscle tissue and 7.5 for the adipose tissue. Reverse transcription was conducted using the QuantiTect Reverse Transcription cDNA Synthesis Kit (Qiagen, Hilden, Germany). Negative results of non-reverse transcription control run were confirmed and cDNA were stored at − 20 °C until use within 4 weeks.
The RT-qPCR was then performed in duplicate with FastStart Essential DNA Green Master (Roche, Basel, Switzerland) in a LightCycler® 96 thermocycler (Roche, Basel, Swiss). The total volume per reaction was 10 μL containing 2 μL cDNA (5 ng/μL), 5 μL FastStart Essential DNA Green Master (Roche), 0.6 μL forward primer and 0.6 μL reverse primer (10 μmol/L stock) and 1.8 μL RNase-free water. Cycling parameters including a pre-incubation of 3 min at 95 °C, followed by a three-step amplification program of 40 cycles consisting of a denaturation, annealing and extension step set at 95 °C for 10 s, 60 °C for 15 s and 72 °C for 30 s, respectively; each gene in either BMSCs, muscle tissue or adipose tissue also included a standard curve for quality purposes (see Supplementary Table 1, Additional file 3). No cq value for the no template control was detected.
The relative quantity of all the candidate reference genes were detected in all samples including the rBMSCs, adipose and muscle tissue with or without chondrogenic or osteogenic induction. The geNorm algorithm (http://medgen.ugent.be/wjvdesomp/geNorm/) was used to evaluate the stability and priority of these candidate reference genes . The raw Cq values of each genes in each sub-study were pre-processed by 2ΔCq algorithm, then the generated data was inputted into geNorm. After the matrix was loaded, a table containing NF of each reference gene was produced, followed by two charts. The first chart showed the sequence of gene stability, in which the stability was improved from left to right, as shown by the decrease of M value. A gene with M < 1.5 is considered as a stable reference gene . The second chart determined the recommended number of the reference genes being used for a specific study, which was indicated by the Vn/n + 1-score. Here, two schemes were compared. Firstly, according to the geNorm algorithm , the value of Vn/n + 1 under 0.15 indicating that no additional reference genes are required for normalization was set as the control scheme. In certain cases, where no Vn/n + 1-score was less than 0.15, the Opt3 or V0.20 were considered as alternatives. Secondly, Vmin was set as the cut-off for choosing the optimal quantity of reference genes.
The relative quantity of osteogenic- or chondrogenic-related target genes
The normalization of each target gene was accomplished by qbase plus software version 3.0 (Biogazelle, Zwijnaarde, Belgium-www.qbaseplus.com), and the results were presented as CNRQ value, which reflect the relative quantity of each target gene based on the selected reference gene set. Upon different schemes, different reference gene sets were used and subsequently different relative quantities of a certain target gene were obtained. All CNRQ values were scaled to the endogenous control that were pure untreated muscle and adipose tissue including rBMSCs.
Normalization factors obtained by different schemes from geNorm were analyzed in GraphPad Prism (GraphPad software Version 5, San Diego, CA) using Spearman rank correlation (correlations with P < 0.05 were considered significant; correlations were very strong when Spearman’s rank correlation coefficient (r) was greater than 0.9). A two-tailed unpaired t-test in GraphPad Prism was used to determine whether different selection schemes of reference genes had significant effects on the normalization of relative expression levels of a certain gene. P < 0.05 values were considered significantly different.
Bhosale AM, Kuiper JH, Johnson WE, Harrison PE, Richardson JB. Midterm to long-term longitudinal outcome of autologous chondrocyte implantation in the knee joint: a multilevel analysis. Am J Sports Med. 2009;37(Suppl 1):131S–8S.
Taylor SC, Nadeau K, Abbasi M, Lachance C, Nguyen M, Fenrich J. The ultimate qPCR experiment: producing publication quality, reproducible data the first time. Trends Biotechnol. 2019;37(7):761–74. https://doi.org/10.1016/j.tibtech.2018.12.002.
Vitorino Carvalho A, Courousse N, Crochet S, Coustham V. Identification of Reference Genes for Quantitative Gene Expression Studies in Three Tissues of Japanese Quail. Genes (Basel). 2019;10(3):197. https://doi.org/10.3390/genes10030197.
Git A, Dvinge H, Salmon-Divon M, Osborne M, Kutter C, Hadfield J, Bertone P, Caldas C. Systematic comparison of microarray profiling, real-time PCR, and next-generation sequencing technologies for measuring differential microRNA expression. RNA. 2010;16(5):991–1006. https://doi.org/10.1261/rna.1947110.
Zhang WX, Fan J, Ma J, Rao YS, Zhang L, Yan YE. Selection of Suitable Reference Genes for Quantitative Real-Time PCR Normalization in Three Types of Rat Adipose Tissue. Int J Mol Sci. 2016;17(6):968. https://doi.org/10.3390/ijms17060968.
Murphy J, Bustin SA. Reliability of real-time reverse-transcription PCR in clinical diagnostics: gold standard or substandard? Expert Rev Mol Diagn. 2009;9(2):187–97. https://doi.org/10.1586/14737188.8.131.52.
Andersen CL, Jensen JL, Orntoft TF. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004;64(15):5245–50. https://doi.org/10.1158/0008-5472.CAN-04-0496.
Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Speleman F. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002;3(7):RESEARCH0034.
Guenin S, Mauriat M, Pelloux J, Van Wuytswinkel O, Bellini C, Gutierrez L. Normalization of qRT-PCR data: the necessity of adopting a systematic, experimental conditions-specific, validation of references. J Exp Bot. 2009;60(2):487–93. https://doi.org/10.1093/jxb/ern305.
Farrokhi A, Eslaminejad MB, Nazarian H, Moradmand A, Samadian A, Akhlaghi A. Appropriate reference gene selection for real-time PCR data normalization during rat mesenchymal stem cell differentiation. Cell Mol Biol (Noisy-le-grand). 2012;58 Suppl:OL1660–70.
Frampton GM, Fichtenholtz A, Otto GA, Wang K, Downing SR, He J, Schnall-Levin M, White J, Sanford EM, An P, Sun J, Juhn F, Brennan K, Iwanik K, Maillet A, Buell J, White E, Zhao M, Balasubramanian S, Terzic S, Richards T, Banning V, Garcia L, Mahoney K, Zwirko Z, Donahue A, Beltran H, Mosquera JM, Rubin MA, Dogan S, Hedvat CV, Berger MF, Pusztai L, Lechner M, Boshoff C, Jarosz M, Vietz C, Parker A, Miller VA, Ross JS, Curran J, Cronin MT, Stephens PJ, Lipson D, Yelensky R. Development and validation of a clinical cancer genomic profiling test based on massively parallel DNA sequencing. Nat Biotechnol. 2013;31(11):1023–31. https://doi.org/10.1038/nbt.2696.
Pombo-Suarez M, Calaza M, Gomez-Reino JJ, Gonzalez A. Reference genes for normalization of gene expression studies in human osteoarthritic articular cartilage. BMC Mol Biol. 2008;9(1):17. https://doi.org/10.1186/1471-2199-9-17.
Schulze F, Malhan D, El Khassawna T, Heiss C, Seckinger A, Hose D, Rosen-Wolff A. A tissue-based approach to selection of reference genes for quantitative real-time PCR in a sheep osteoporosis model. BMC Genomics. 2017;18(1):975. https://doi.org/10.1186/s12864-017-4356-4.
Okamura K, Inagaki Y, Matsui TK, Matsubayashi M, Komeda T, Ogawa M, Mori E, Tanaka Y. RT-qPCR analyses on the osteogenic differentiation from human iPS cells: an investigation of reference genes. Sci Rep. 2020;10(1):11748. https://doi.org/10.1038/s41598-020-68752-2.
An Y, Reimers K, Allmeling C, Liu J, Lazaridis A, Vogt PM. Validation of differential gene expression in muscle engineered from rat groin adipose tissue by quantitative real-time PCR. Biochem Biophys Res Commun. 2012;421(4):736–42. https://doi.org/10.1016/j.bbrc.2012.04.073.
Lee WJ, Jeon RH, Jang SJ, Park JS, Lee SC, Baregundi Subbarao R, Lee SL, Park BW, King WA, Rho GJ. Selection of reference genes for quantitative gene expression in porcine Mesenchymal stem cells derived from various sources along with differentiation into multilineages. Stem Cells Int. 2015;2015:235192.
Zhang G, Zhao M, Song C, Luo A, Bai J, Guo S. Characterization of reference genes for quantitative real-time PCR analysis in various tissues of Anoectochilus roxburghii. Mol Biol Rep. 2012;39(5):5905–12. https://doi.org/10.1007/s11033-011-1402-1.
He T, Hausdorf J, Chevalier Y, Klar RM. Trauma induced tissue survival in vitro with a muscle-biomaterial based osteogenic organoid system: a proof of concept study. BMC Biotechnol. 2020;20(1):8. https://doi.org/10.1186/s12896-020-0602-y.
Huang Y, Seitz D, König F, Müller PE, Jansson V, Klar RM. Induction of articular Chondrogenesis by chitosan/hyaluronic-acid-based biomimetic matrices using human adipose-derived stem cells. Int J Mol Sci. 2019;20(18):4487. https://doi.org/10.3390/ijms20184487.
Dheda K, Huggett JF, Chang JS, Kim LU, Bustin SA, Johnson MA, Rook GA, Zumla A. The implications of using an inappropriate reference gene for real-time reverse transcription PCR data normalization. Anal Biochem. 2005;344(1):141–3. https://doi.org/10.1016/j.ab.2005.05.022.
Goossens K, Van Poucke M, Van Soom A, Vandesompele J, Van Zeveren A, Peelman LJ. Selection of reference genes for quantitative real-time PCR in bovine preimplantation embryos. BMC Dev Biol. 2005;5(1):27. https://doi.org/10.1186/1471-213X-5-27.
Smits K, Goossens K, Van Soom A, Govaere J, Hoogewijs M, Vanhaesebrouck E, Galli C, Colleoni S, Vandesompele J, Peelman L. Selection of reference genes for quantitative real-time PCR in equine in vivo and fresh and frozen-thawed in vitro blastocysts. BMC Res Notes. 2009;2(1):246. https://doi.org/10.1186/1756-0500-2-246.
Tong Z, Gao Z, Wang F, Zhou J, Zhang Z. Selection of reliable reference genes for gene expression studies in peach using real-time PCR. BMC Mol Biol. 2009;10(1):71. https://doi.org/10.1186/1471-2199-10-71.
Erkens T, Van Poucke M, Vandesompele J, Goossens K, Van Zeveren A, Peelman LJ. Development of a new set of reference genes for normalization of real-time RT-PCR data of porcine backfat and longissimus dorsi muscle, and evaluation with PPARGC1A. BMC Biotechnol. 2006;6(1):41. https://doi.org/10.1186/1472-6750-6-41.
Ragni E, Vigano M, Rebulla P, Giordano R, Lazzari L. What is beyond a qRT-PCR study on mesenchymal stem cell differentiation properties: how to choose the most reliable housekeeping genes. J Cell Mol Med. 2013;17(1):168–80. https://doi.org/10.1111/j.1582-4934.2012.01660.x.
Lu Y, Yuan M, Gao X, Kang T, Zhan S, Wan H, Li J. Identification and validation of reference genes for gene expression analysis using quantitative PCR in Spodoptera litura (Lepidoptera: Noctuidae). PLoS One. 2013;8(7):e68059. https://doi.org/10.1371/journal.pone.0068059.
Ren B, Betz VM, Thirion C, Salomon M, Jansson V, Muller PE, Betz OB. Osteoinduction within BMP-2 transduced muscle tissue fragments with and without a fascia layer: implications for bone tissue engineering. Gene Ther. 2019;26(1–2):16–28. https://doi.org/10.1038/s41434-018-0047-2.
Schwanhausser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, Chen W, Selbach M. Global quantification of mammalian gene expression control. Nature. 2011;473(7347):337–42. https://doi.org/10.1038/nature10098.
Clements DN, Vaughan-Thomas A, Peansukmanee S, Carter SD, Innes JF, Ollier WE, Clegg PD. Assessment of the use of RNA quality metrics for the screening of articular cartilage specimens from clinically normal dogs and dogs with osteoarthritis. Am J Vet Res. 2006;67(8):1438–44. https://doi.org/10.2460/ajvr.67.8.1438.
Wehrli BM, Huang W, De Crombrugghe B, Ayala AG, Czerniak B. Sox9, a master regulator of chondrogenesis, distinguishes mesenchymal chondrosarcoma from other small blue round cell tumors. Hum Pathol. 2003;34(3):263–9. https://doi.org/10.1053/hupa.2003.41.
Yang Q, Peng J, Guo Q, Huang J, Zhang L, Yao J, Yang F, Wang S, Xu W, Wang A, Lu S. A cartilage ECM-derived 3-D porous acellular matrix scaffold for in vivo cartilage tissue engineering with PKH26-labeled chondrogenic bone marrow-derived mesenchymal stem cells. Biomaterials. 2008;29(15):2378–87. https://doi.org/10.1016/j.biomaterials.2008.01.037.
Dheda K, Huggett JF, Bustin SA, Johnson MA, Rook G, Zumla A. Validation of housekeeping genes for normalizing RNA expression in real-time PCR. Biotechniques. 2004;37(1):112–4, 116, 118-119. https://doi.org/10.2144/04371RR03.
Conception and design: T.H. and F.X; Administrative support: C.Z. and R.M.K; Provision of study materials or patients: T.H. and F.X; Collection and assembly of data: T.H., F.X. and X.C; Data analysis and interpretation: All authors; Manuscript writing: All authors; Final approval of manuscript: All authors.
All practical experimental steps were performed in keeping with the rules and regulations of the Animal Protection Laboratory Animal Regulations (2013) and approved by the Animal Care Committee of Renji Hospital (Shanghai, China, No.: 201606230235). In addition, the study was carried out in accordance with Arrive guidelines.
Consent for publication
The authors have no conflicts of interest to declare.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.