Optimizations for Identifying Reference Genes in Bone and Cartilage Bioengineering

doi:10.21203/rs.3.rs-96393/v1

Download PDF

Research Article

Optimizations for Identifying Reference Genes in Bone and Cartilage Bioengineering

https://doi.org/10.21203/rs.3.rs-96393/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 16 Mar, 2021

Read the published version in BMC Biotechnology →

Version 1

posted

You are reading this latest preprint version

Background：Quantitative reverse transcription polymerase chain reaction (qRT-PCR) is commonly considered as the best-established technique for gene expression assay. However, appropriate reference gene set selection remains the critical and challenging subject for the proper understanding of gene expression pattern. Mixed opinions pertain in how to choose optimal reference gene set according 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 was the most feasible for the identification of reference genes in bone and cartilage bioengineering experiments. 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 ensuring the stable and accurate interpretation of gene expression in bone and cartilage bioengineering.

Results：The stability and pairwise variance of eight candidate reference genes were analyzed using geNorm. The V_0.15- vs. V_min-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 V_0.15 and V_min schemes and the variance of target genes expression levels generated by these two schemes showed different schemes may 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 (Opt₃), V_0.20 and V_min schemes, and the analysis of NF values and eventual normalization of target genes showed that the different selection schemes have significant effect on the eventual normalization of target genes.

Conclusions：Based on the results, the proposed cut-off value of V_n/n+1 under 0.15, according to geNorm algorithm, should be considered with caution, especially in certain tissue types such as skeletal muscle and adipose tissue. Instead, the minimum V_n/n+1 should be used as the cut-off for choosing the optimal reference gene set especially when the stability and variation of candidate reference genes in a specific study are unclear.

Biotechnology and Bioengineering

qRT-PCR

reference gene set

osteogenesis

chondrogenesis

bioengineering

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 obviate the necessity of autograft utilization [1], inconsistent clinical outcomes of various methods including autologous chondrocyte transplantation, an intervention approach being used for decades [2], have shown the scientific community that bioengineering techniques and technology is still flawed. 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. Quantitative reverse transcription real-time polymerase chain reaction (qRT-PCR) has emerged as one of those fundamental and best-established technique, to date, critical for bioengineering principles [3] that if used properly, can generate reliable, comparable and reproducible gene expression data with high sensitivity and accuracy [4]. Thus, qRT-PCR is not only regarded as the benchmark for gene analysis [3, 5] but also is used to validate the results of NextGen sequencing and microarray assays [6, 7]. Nevertheless, improper optimization and standardization may significantly increase the variability of gene expression results thereby impairing the reproducibility and subsequently compromising the translation efficiency of present bioengineering techniques [3, 8-10].

To prevent the misinterpretation of results caused by non-standard experimental procedures and details [11], 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 [8, 12-14]. Particularly, the selection of reference genes has drawn considerable attention of researchers, because the expressions of these so-called “ideal” non-variant genes are unstable under certain conditions [13, 15]. Hence, the proper combination of multiple reference genes was established to minimize the instability and variation [8]. 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 this technique. [16, 17]. 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 each specific experiment [18], 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 (V_n/n+1-score ) [13]. The value of V_n/n+1with 0.15 (V_0.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 [13]. Besides the V_n/n+1-score-based scheme, the GeNorm algorithm also provide an alternative to select the best three reference genes based on geNorm results. Moreover, other schemes were raised due to the limitation of the V_0.15 method [19-21]. In some cases, 0.20 was set as a trade-off when in some tissue types (e.g. adipose tissue) the minimum V_n/n+1-score (V_min) was higher than 0.15 [19, 20]. In another group, Lee et al. [22] evaluated the potential benefit of using the V_min as threshold compared with the optimal three genes. Hence, debates persist on the schemes for choosing reference genes to improve the accuracy of qRT-PCR assays [19, 21-24].

Therefore, the purpose of this study was to investigate which schemes were the most feasible for the identification of reference genes in bone and cartilage bioengineering experiments. Our hypothesis was that identification of reference genes based on V_min was the optimal scheme for the normalization of qRT-PCR, gaining an accurate and reliable gene expression data in bone and cartilage bioengineering.

2.1 V_0.15- vs. V_min-based normalization scheme for analyzing gene expression data in rBMSCs

The correlation of NF values between the V_0.15 and V_min 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 V_0.15 was 0.078, while the V_min was 0.055 (Figure 1A). Hench, combining the sequencing of eight candidate reference genes based on M value, the V_0.15-based reference gene set contained ribosomal protein L13α (RPL13α) and actin beta (ACTB), while the V_min-based reference gene set contained RPL13α, ACTB and RNA polymerase II subunit e (POLR2e) (Figure 1B). The r value in the Spearman rank correlation analysis between NF_V0.15 and NF_Vminwas 0.9762 (Figure 1C), which showed a very strong correlation between these two schemes in terms of rBMSCs. Furthermore, the variance of target genes expression levels using these two schemes were compared. The osteogenic-related target genes, including bone morphogenetic protein-2 (BMP-2), BMP-6, osteocalcin (OCN) and runt-related transcription factor 2 (RUNX-2) were normalized by 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 V_0.15 or V_min scheme, did not show a significant difference (P > 0.05) (Figure 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, V_0.15 was 0.068, while V_min was 0.058 (Figure 2A). Hench, combining the sequences of eight candidate reference genes based on M value, the V_0.15-based reference gene set contained ACTB and RPL13α, and The V_min-based reference gene set contained RPL13α, ACTB and POLR2e (Figure 2B). The r value in the Spearman rank correlation analysis between NF_V0.15 and NF_Vmin was 0.9524 (Figure 2C). This result showed that there was no 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-β₁) and TGF-β₃, were normalized by two reference gene sets, and the CNRQ values were obtained. The relative expression levels of ACAN, SOX9, TGF-β₁ and TGF-β₃ did not show a significant deference between the V_0.15- and V_min- based CNRQ values (P > 0.05), which revealed that different selection schemes have no significant effect on the eventual normalization of target genes (Figure 2D).

2.2 V_0.15- vs. V_min-based normalization scheme for analyzing gene expression data in skeletal muscle tissue

The correlation of NF values between the V_0.15 and V_min 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 V_0.15 was 0.108, while V_min was 0.067 (Figure 3A). Hench, combining the sequencing of eight candidate reference genes based on M value, the V_0.15-based reference gene set contained POLR2e and RPL13α, while the V_min-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) (Figure 3B). The r value in the Spearman rank correlation analysis between NF_V0.15and NF_Vmin was 0.6429 (Figure 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 RUNX-2, were normalized by the two reference gene sets, and the CNRQ values were obtained. The relative expression levels of OCN and RUNX-2 showed no significant deference between the V_0.15- and V_min- 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 target genes (Figure 3D).

In the chondrogenic sub-study, the V_0.15 was 0.079, while the V_min was 0.057 (Figure 4A). Hench, combining the sequencing of eight candidate reference genes based on M value, the V_0.15-based reference gene set contains GAPDH and RPLP0, and the V_min-based reference gene set contains RPLP0, GAPDH, ACTB, POLR2e and TBP (Figure 4B). The correlation coefficient in the Spearman rank correlation analysis between NF_V0.15 and NF_Vmin was 0.833. This result showed that there was rare correlation between these two selection schemes (Figure 4C). The chondrogenic-related target genes, including ACAN, SOX9, TGF-β₁ and TGF-β₃, were normalized by the two reference gene sets, and the CNRQ values were obtained. The relative expression levels of ACAN, SOX9, TGF-β₁ and TGF-β₃ showed significant deference between the V_0.15- and V_min- based CNRQ values (P < 0.05), which revealed that different selection schemes have a significant effect on the eventual normalization of target genes (Figure 4D).

2.3 V_0.20/Opt3 vs. V_min-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 optimal three reference genes (Opt₃), V_0.20 and V_minschemes. The correlation of NF values between the V_0.20 and V_minschemes, as long as Opt3 and V_min schemes was analyzed and the variance of target genes expression levels generated by these three schemes was compared in terms of the muscle tissue.

In the osteogenic sub-study, V_0.20 was 0.194, while the V_min was 0.170 (Figure 5A). Hench, combining the sequencing of eight candidate reference genes based on M value, the Opt3-based reference gene set contained RNA 28S ribosomal 4 (RNA28S4), GAPDH and ACTB, the V_0.20-based reference gene set contained RNA28S4, GAPDH, ACTB and RPL13α, while the V_min-based reference gene set contained RNA28S4, GAPDH, ACTB, RPL13α, POLR2e and TBP (Figure 5B). The r value in the Spearman rank correlation analysis between NF_opt3 and NF_Vmin was 0.833, with that of NF_V0.20 and NF_Vmin was 0.762 (Figure 5C, E). These results showed that there was neither a correlation between V_min and V_0.20 schemes nor between V_min and Opt3 schemes. The osteogenic-related target genes, including BMP-2, BMP-6, OCN and RUNX-2, were normalized by three reference gene sets, and the CNRQ values were obtained. The relative expression levels of BMP-2, BMP-6, OCN and RUNX-2 all showed significant deference between the Opt3- and V_min- based CNRQ values or between the V_0.20- and V_min- based CNRQ values (P < 0.05), which revealed that different selection schemes have significant effects on the eventual normalization of target genes (Figure 5D, F).

In the chondrogenic sub-study, V_0.20was 0.181, while the V_min was 0.168 (Figure 6A). Hench, combining the sequencing of eight candidate reference genes based on M value, the Opt3-based reference gene set contained RPLP0, GAPDH and POLR2e, the V_0.20-based reference gene set contained RPLP0, GAPDH, POLR2e and RPL13α, while the V_min-based reference gene set contained RPLP0, GAPDH, POLR2e, RPL13α, RNA28S4, TBP and succinate dehydrogenase complex flavoprotein subunit A (SDHA) (Figure 6B). The r value in the Spearman rank correlation analysis between NF_opt3 and NF_Vmin was 0.4286, and meanwhile between NF_V0.20 and NF_Vmin was 0.6429 (Figure 6C, E). These results showed that there was neither a correlation between V_min and V_0.20 schemes nor between V_min and Opt3 schemes. The chondrogenic-related target genes, including ACAN, SOX9, TGF-β₁ and TGF-β₃, were normalized by three reference gene sets, and the CNRQ values were obtained. The relative expression levels of ACAN, SOX9, TGF-β₁ and TGF-β₃ showed significant deference between the V_opt3- and V_min- based CNRQ values or between the V_0.20- and V_min- based CNRQ values (P < 0.05), which revealed that different selection schemes have significant effect on the eventual normalization of target genes (Figure 6D, F).

Accurate gene analysis remains essential for objective evaluation of the efficacy of bioengineering [25, 26]. Although the establishment of MIQE and the advocacy of multiple reference genes increased the reliability of the results [8, 13], there are still no unified solutions for the choice of the quantity and priority of reference genes, which make the results across the bioengineering spectrum inconsistent [27, 28]. Several mathematical algorithms to solve this dilemma aimed at detecting the stability of reference genes have been developed [12, 29, 30], such as geNorm [13], NormFinder [12], and BestKeeper [31]. In this study, we attempted to optimize the scheme of using geNorm for reference genes determination in which we recommended that 0.15, commonly regarded as the cut-off point in the pairwise variation analysis, should be adjusted to the minimum V-score, especially in the analysis of RNA extracted from tissues.

By integrating the results of GeNorm, NormFinder, and BestKeeper, Tong et al. [32] and Zhang et al.[33] identified suitable reference gene(s) for normalization of gene expression in their own specific experiments. It seems that the conjunction of different methods is popular [12, 31-33], however, in our previous study, we indicated that the programs utilized were irrelevant with regards to reference genes evaluation of stability and the determination of quantity, which should not affect which reference genes are selected and the optimum number to use in an experiment [18]. Compared with the other algorithms, a more user-friendly interface, a more intuitive data presentation and the infinite number of samples including candidate genes, allows GeNorm to quite accessible [13]. 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 [18, 34], 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 [13, 18].

In terms of GeNorm, the M value reflects the stability of each reference gene, which the lowest value embodies as the most stable expression, while the determination of the optimal number of reference genes for accurate gene expression analysis relied on the pattern of pairwise variation V-score [8]. The usual V-curve showed that adding reference genes was a double-edged sword in the normalization process. The nonspecific variation was eliminated by adding the stable reference genes proven by the decrease in V-score, but 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 did not contribute to minimize 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 (V_n/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 [13]. Vandesompele et al. [13] 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 V_0.15 and V_min (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 NF_V0.15 with NF_Vmin(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 V_min was significant. Erkens et al. [35] and Pérez et al. [24] all attempted to evaluate the V_min-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.[36] 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 (V_min- or V_0.15-based), but not in the milieu of all target gene in rBMSCs. The results revealed that the variation erased by V_min based sets cannot be neglected in muscle tissue, which was consistent with the outcomes in Spearman rank correlation analysis. On the base of the current data, regarding rBMSCs, the selection of reference genes based on V_0.15 was still applicable and reasonable, but for muscle tissue, the conditions need be stricter, where we recommend using the V_min 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 V_0.15 based scheme inapplicable. Facing a similar dilemma, Ayers et al. [19] 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. [20] also drew on this new cut-off (V_0.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 [20]. In our study, the correlation between two NFs which corresponded to the V_min and V_0.2 schemes were also weak, suggesting that V_min, 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. [13] and Lu et al. [37] 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 V_min 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. [13] 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 [38, 39]. Studies have shown that for different cell types, the expression of reference genes was not as stable [18, 40], 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. [38] 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 [41]. 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 [42].

From the present study, we recommend that the proposed 0.15 cut-off value, according to the GeNorm algorithm, should not be taken as a too strict cut-off, especially when utilizing certain tissue types. 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. Indeed, if the stability and variation of candidate reference genes, in a specific study are unclear, the minimum V-score should be used as this provides a superior selection for the optimal number and category of reference genes needed to generate the correct 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.

5.1 Study Setup

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 qRT-PCR. Subsequently, the stability and pairwise variance of candidate reference genes were analyzed using geNorm [13]. 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 (Figure 7).

5.2 Cell and Tissue Specimens

RBMSCs (passage 0, Sciencell, Carlsbad, CA, USA), skeletal muscle and adipose tissue from F-344 adult female rat (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).

5.3 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). Under sterile conditions, fresh abdominal muscle and subcutaneous adipose tissue was harvested and placed temporarily in Dulbecco’s modified eagles 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 5mm 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 48h at 37 °C containing 5 % CO₂. 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.

5.4 Cell culture

RBMSCs were used for the cellular culturing part, which were seeded at a density of 2*10⁴ per monolayer flask (Thermo Fisher Scientific) and cultured in DMEM supplemented with 3% P/S at 37 °C containing 5 % CO₂ 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 2^nd passage cells reached 80 % of confluence, they were passaged for the following induction of differentiation procedures. Some of the rBMSCs in 2^nd passage pure without culturing were collected immediately as these would be used as the endogenous normalization control in downstream analysis procedures.

5.5 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 human BMP-6 (R&D Systems, Minneapolis, MN,USA), 10 ng/mL recombinant human TGF-β₃(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) [43-46]; 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) [47]; the normal medium was DMEM supplemented with 3% P/S.

After 48h 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 two days, collected and then stored at -80 °C until further use.

Similarly, rBMSCs (passage 2) once having reached 80% confluence were trypsinized and seeded at 2*10⁴ 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 two 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.

5.6 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 qRT-PCR, all with a standard deviation of the average amplification threshold cycle quantification value (Cq) less than 1 across 35 in rat tissues [48, 49]. 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-β₁ and TGF-β₃. Meanwhile, osteogenic-related genes included BMP-2, BMP-6, OCN and RUNX-2. 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). All the primer sequences were presented in Table 1.

As previously established [18], 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&LINK_LOC=blasthome), against the GenBank database (https://www.ncbi.nlm.nih.gov/genbank/) to validate primer reference gene sequence amplification specificity.

5.7 QRT-PCR and GeNorm assessment

Total RNA was extracted from cells and tissue samples by first homogenizing the material either with a Micro-Dismembrator S (Sartorius Stedim Biotech, Göttingen, Germany) or liquid nitrogen in conjunction with a mortar and pestle, respectively, and then using the RNeasy Fibrous Tissue Mini Kit (Qiagen, Hilden, Germany) following the manufacturer ́s protocol. Total extracted RNA was assessed using a NanoDropTM Lite (Thermo Scientific, Waltham, USA). Reverse transcription was conducted using the QuantiTect Reverse Transcription cDNA Synthesis Kit (Qiagen, Hilden, Germany), and cDNA were stored at -20°C until use.

The qRT-PCR 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; and a final extension at 72 °C for 5 min.

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 [50]. The raw Cq values of each genes in each sub-study were pre-processed by 2^Δ^C^q 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 [13]. The second chart determined the recommended number of the reference genes being used for a specific study, which was indicated by the V_n/n+1-score. Here, two schemes were compared. Firstly, according to the GeNorm algorithm [13], the value of V_n/n+1under 0.15 indicating that no additional reference genes are required for normalization was set as the control scheme. In certain cases, where no V_n/n+1-score was less than 0.15, the Opt₃ or V_0.20 were considered as alternatives. Secondly, V_min was set as the cut-off for choosing the optimal quantity of reference genes.

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

5.9 Statistics

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.

qRT-PCR	quantitative reverse transcription real-time polymerase chain reaction
MIQE	the Minimum Information for Publication of Quantitative Real-Time PCR Experimen
NF	normalization factor
rBMSCs	rat bone marrow mesenchymal stem cells
RPL13α	ribosomal protein L13α
ACTB	actin beta
POLR2e	RNA polymerase II subunit e
BMP	bone morphogenetic protein
OCN	osteocalcin
RUNX-2	runt-related transcription factor 2
CNRQ	calibrated normalized relative quantity
ACAN	aggrecan
SOX9	sex determining region Y-box 9
TGF-β	transforming growth factor beta
GAPDH	glyceraldehyde 3-phosphate dehydrogenase
TBP	TATA-binding protein
RPLP0	ribosomal protein lateral stalk subunit P0
Opt3	optimal three reference genes
RNA28S4	RNA 28S ribosomal 4
SDHA	succinate dehydrogenase complex flavoprotein subunit A
DMEM	Dulbecco’s modified eagles medium
P/S	penicillin/streptomycin
Cq	cycle quantification value

Ethics approval and consent to participate: 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).

Consent for publication: All authors consent the article for publication.

Data and materials availability: The necessary algorithmic codes of the program GeNorm are readily available at (http://medge n.ugent.be/wjvdesomp/genorm/). All data, raw and processed, is readily available from the corresponding author on request.

Conflicts of Interest: The authors declare that they have no competing interests.

Funding: This research was funded internally and by the CHINA SCHOLARSHIP COUNCIL, grant number 201703920036 (F.X.), 201708140085 (X.C) and 201606230235 (T.H.) and the APC was funded internally.

Author Contributions: 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.

Acknowledgements: Not applicable.

Richardson SM, Hoyland JA, Mobasheri R, Csaki C, Shakibaei M, Mobasheri A: Mesenchymal stem cells in regenerative medicine: opportunities and challenges for articular cartilage and intervertebral disc tissue engineering. J Cell Physiol 2010, 222(1):23-32.
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-138S.
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.
Bustin SA: Why the need for qPCR publication guidelines?--The case for MIQE. Methods 2010, 50(4):217-226.
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).
Fang Z, Cui X: Design and validation issues in RNA-seq experiments. Brief Bioinform 2011, 12(3):280-287.
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.
Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, Mueller R, Nolan T, Pfaffl MW, Shipley GL et al: The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem 2009, 55(4):611-622.
Sanders R, Mason DJ, Foy CA, Huggett JF: Considerations for accurate gene expression measurement by reverse transcription quantitative PCR when analysing clinical samples. Anal Bioanal Chem 2014, 406(26):6471-6483.
Sanders R, Bustin S, Huggett J, Mason D: Improving the standardization of mRNA measurement by RT-qPCR. Biomol Detect Quantif 2018, 15:13-17.
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-197.
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-5250.
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.
Pfaffl MW, Horgan GW, Dempfle L: Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res 2002, 30(9):e36.
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-493.
Frampton GM, Fichtenholtz A, Otto GA, Wang K, Downing SR, He J, Schnall-Levin M, White J, Sanford EM, An P et al: Development and validation of a clinical cancer genomic profiling test based on massively parallel DNA sequencing. Nat Biotechnol 2013, 31(11):1023-1031.
D'Haene B, Vandesompele J, Hellemans J: Accurate and objective copy number profiling using real-time quantitative PCR. Methods 2010, 50(4):262-270.
He T, Huang Y, Chak JC, Klar RM: Recommendations for improving accuracy of gene expression data in bone and cartilage tissue engineering. Sci Rep 2018, 8(1):14874.
Ayers D, Clements DN, Salway F, Day PJ: Expression stability of commonly used reference genes in canine articular connective tissues. BMC Vet Res 2007, 3:7.
Hosseini A, Sauerwein H, Mielenz M: Putative reference genes for gene expression studies in propionate and beta-hydroxybutyrate treated bovine adipose tissue explants. J Anim Physiol Anim Nutr (Berl) 2010, 94(5):e178-184.
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-742.
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-5912.
Perez R, Tupac-Yupanqui I, Dunner S: Evaluation of suitable reference genes for gene expression studies in bovine muscular tissue. BMC Mol Biol 2008, 9:79.
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 Biotechnology 2020, 20(1):8.
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. International journal of molecular sciences 2019, 20(18):4487.
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-143.
Eisenberg E, Levanon EY: Human housekeeping genes, revisited. Trends Genet 2013, 29(10):569-574.
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: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:246.
Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP: Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper--Excel-based tool using pair-wise correlations. Biotechnol Lett 2004, 26(6):509-515.
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:71.
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).
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 molecular biology 2009, 10:71-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: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-180.
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.
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.
Bach AD, Beier JP, Stern-Staeter J, Horch RE: Skeletal muscle tissue engineering. J Cell Mol Med 2004, 8(4):413-422.
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-342.
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-1444.
Hunter MJ, Commerford SL: Pressure homogenization of mammalian tissues. Biochim Biophys Acta 1961, 47:580-586.
Huang Y, Seitz D, Konig F, Muller 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).
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-269.
Estes BT, Guilak F: Three-dimensional culture systems to induce chondrogenesis of adipose-derived stem cells. Methods Mol Biol 2011, 702:201-217.
Yang Q, Peng J, Guo Q, Huang J, Zhang L, Yao J, Yang F, Wang S, Xu W, Wang A et al: 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-2387.
Sottile V, Thomson A, McWhir J: In vitro osteogenic differentiation of human ES cells. Cloning Stem Cells 2003, 5(2):149-155.
Warrington JA, Nair A, Mahadevappa M, Tsyganskaya M: Comparison of human adult and fetal expression and identification of 535 housekeeping/maintenance genes. Physiol Genomics 2000, 2(3):143-147.
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-114, 116, 118-119.
VanGuilder HD, Vrana KE, Freeman WM: Twenty-five years of quantitative PCR for gene expression analysis. Biotechniques 2008, 44(5):619-626.

Table 1. Gene primers for Rattus norvegicus with accession number and sequences.

	Gene	Accession Number	5’ sequence	3’ sequence
Reference genes	TBP	BC081939.1	TAACCCAGAAAGTCGAAGAC	CCGTAAGGCATCATTGGA
	GAPDH	BC083511.1	CATGGGTGTGAACCATGA	TGTCATGGATGACCTTGG
	POLR2e	BC158787.1	GACCATCAAGGTGTACTGC	CAGCTCCTGCTGTAGAAAC
	RPLP0	BC001834.2	CAACCCAGCTCTGGAGA	CAGCTGGCACCTTATTGG
	SDHA	NM_130428.1	GCGGTATGAGACCAGTTATT	CCTGGCAAGGTAAACCAG
	RPL13α	NM_173340.2	TTTCTCCGAAAGCGGATG	AGGGATCCCATCCAACA
	ACTB	NM_031144.3	AGCTATGAGCTGCCTGA	GGCAGTAATCTCCTTCTGC
	RNA28S4	NR_145822.1	GCGGCCAAGCGTTCATA	CCTGTCTCACGACGGTCTAA

Genes of interest	TGF-β₁	NM_021578.2	TTTAGGAAGGACCTGGGTT	ACCCACGTAGTAGACGATG
	TGF-β₃	NM_013174.2	AACCTAAGGGTTACTATGCC	ACCACCATGTTGGACAG
	SOX9	NM_080403.1	CCAGAGAACGCACATCAAG	ATACTGATGTGGCTGGTGG
	RUNX2	NM_001278484.2	CCCAAGTGGCCACTTAC	CTGAGGCGGTCAGAGA
	ACAN	NM_022190.1	CAAGTGGAGCCGTGTTT	TTTAGGTCTTGGAAGCGAG
	BMP-6	NM_013107.1	GGACATGGTCATGAGCTTTG	GTCAGAGTCTCTGTGCTGAT
	BMP-2	NM_017178.1	GGAAGTGGCCCACTTAGA	TCACTAGCAGTGGTCTTACC
	OCN	NM_013414.1	ACCTGGCAGGTGCAAA	CTCACACACCTCCCTGTG
TBP: TATA-binding protein; GAPDH: Glyceraldehyde 3-phosphate dehydrogenase; POLR2e: RNA polymerase II subunit e; RPLP0: Ribosomal protein lateral stalk subunit P0; SDHA: Succinate dehydrogenase complex flavoprotein sub-unit A; ROL13α: Ribosomal protein L13 α; ACTB: Actin beta; RNA28S4: RNA 28S ribosomal 4; TGF-β1: transforming growth factor beta 1; TGFß-3, transforming growth factor beta 3; SOX9: Sex determining region Y-box 9; RUNX2: Runt-related transcription factor 2; ACAN: Aggrecan; BMP-2: Bone morphogenetic protein 2; BMP-6: Bone morphogenetic protein 6; OCN: Osteocalcin.

Download PDF

Journal Publication

published 16 Mar, 2021

Read the published version in BMC Biotechnology →

Version 1

posted

You are reading this latest preprint version

Optimizations for Identifying Reference Genes in Bone and Cartilage Bioengineering

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Results

Discussion

Conclusions

Methods

Abbreviations

Declarations

References

Tables

Status:

Journal Publication

Version 1