Bioinformatics analysis of the biological changes involved in the osteogenic differentiation of human mesenchymal stem cells

Abstract The mechanisms underlying the osteogenic differentiation of human bone marrow mesenchymal stem cells (hBMSCs) remain unclear. In the present study, we aimed to identify the key biological processes during osteogenic differentiation. To this end, we downloaded three microarray data sets from the Gene Expression Omnibus (GEO) database: GSE12266, GSE18043 and GSE37558. Differentially expressed genes (DEGs) were screened using the limma package, and enrichment analysis was performed. Protein‐protein interaction network (PPI) analysis and visualization analysis were performed with STRING and Cytoscape. A total of 240 DEGs were identified, including 147 up‐regulated genes and 93 down‐regulated genes. Functional enrichment and pathways of the present DEGs include extracellular matrix organization, ossification, cell division, spindle and microtubule. Functional enrichment analysis of 10 hub genes showed that these genes are mainly enriched in microtubule‐related biological changes, that is sister chromatid segregation, microtubule cytoskeleton organization involved in mitosis, and spindle microtubule. Moreover, immunofluorescence and Western blotting revealed dramatic quantitative and morphological changes in the microtubules during the osteogenic differentiation of human adipose‐derived stem cells. In summary, the present results provide novel insights into the microtubule‐ and cytoskeleton‐related biological process changes, identifying candidates for the further study of osteogenic differentiation of the mesenchymal stem cells.


| INTRODUC TI ON
Bone defects caused by bone trauma, tumours and infection are common and refractory diseases encountered in orthopaedic clinical settings. The reconstruction of bone defects is challenging. Currently, the most common treatment strategy is repair by bone transplantation; however, limitations of the donor restrict its application, thereby affecting patient outcomes. Although allogeneic bone transplantation does not suffer limitations arising from a lack of suitable donors, this strategy is associated with antigenicity, which often results in transplant failure due to severe immune rejection. 1,2 In this context, bone tissue engineering is an indispensable strategy because it offers clinicians the choice of a variety of artificial bone substitutes made of metal, ceramic or polymer, thus providing great convenience for clinical treatment (eg reducing the transplant failure risk). 3,4 Bone tissue engineering requires three elements: seed cells, scaffold materials and bone inducers. Mesenchymal stem cells (MSCs) such as human bone MSCs (hBMSCs) and human adipose-derived stem cells (hASCs) exhibit the potential for self-proliferation and differentiation into adipogenic, osteogenic, chondrogenic and other lineages. Both cell types represent ideal seed cells for bone tissue engineering. 5,6 The osteogenic differentiation of stem cells is the basis for bone formation. Elucidation of the mechanisms underlying the osteogenic differentiation of seed cells should enable the development of techniques aimed at the osteogenic regulation of seed cells, with potential application in the construction of biomaterials for application in bone tissue engineering. 3,6 Both hBMSCs and hASCs are commonly used as seed cells in bone tissue engineering applications. Increasing evidence shows that the osteogenic differentiation of MSCs involves a series of signalling pathways [7][8][9] such as BMP-SMAD, 10,11 WNT/Catenin, 12 Notch 13,14 and MAPK,, 15 and complex regulatory networks formed by interactions between these pathways. 16,17 Transcription factors, such as TWIST and MSX2, may also be involved in the regulation of osteogenic differentiation. 18,19 In addition, certain physical stimuli promote the osteogenic differentiation of stem cells; for example, Heydari demonstrated that substrate stiffness and exposure to electromagnetic fields enhance the osteogenic potential of stem cells in the absence of chemical stimulation. [20][21][22] Therefore, elucidation of the precise molecular mechanisms underlying osteogenesis is crucial to the development of bone tissue engineering applications and treatment strategies for bone defects.
Microarray techniques and bioinformatics analysis have been widely used to screen for genome-level differences involved in osteogenic differentiation, enabling the identification of differentially expressed genes (DEGs) and functional pathways associated with osteogenic differentiation of stem cells. [23][24][25][26] However, the false-positive rates of independent microarray analysis data make it difficult to obtain reliable results. In this study, three mRNA microarray data sets were downloaded from Gene Expression Omnibus (GEO) for analysis to identify DEGs between the control group and the induction group. We then carried out gene ontology (GO), Kyoto genome and genome encyclopaedia (KEGG) pathway enrichment analysis and protein-protein interaction (PPI) network analysis to elucidate the underlying molecular mechanisms. In summary, a total of 240 DEGs, 10 hub genes and one important biological process (microtubule-related) were identified, and the changes in microtubules may be a key factor in osteogenic differentiation of mesenchymal stem cells.

| Identification of DEGs
We combined samples from the control group and the induction group in three data sets, and removed or averaged the probe sets without corresponding gene symbols or the genes with multiple probe sets, respectively. Then, the limma package was used to remove batch effect and identify DEGs. LogFC > 0.584963 (ie FC > 1. 5) and P-value < .05 were considered to indicate statistical significance.

| KEGG and GO enrichment analyses of DEGs
Metascape (https://metas cape.org/gp/index.html#/main/step1) 30 is an analytical website that combines functional enrichment, interactome analysis, gene annotation and membership search to leverage over 40 independent knowledgebases within one integrated portal. KEGG is a database resource for elucidating high-level functions and effects of the biological system. 31,32 Gene Ontology (GO) is a major bioinformatics initiative that for high-quality functional gene annotation and analysing gene biological processes (BP), molecular functions (MF) and cellular components (CC). 33 Metascape was used for analysing the function of DEGs. Min overlap = 3 and Min Enrichment = 1.5 were the screening conditions. P < .01 was considered statistically significant.

| PPI network construction and module analysis
The PPI network was analysed using the Search Tool for the Retrieval of Interacting Genes (STRING; http://strin g-db.org). Analysis of functional interactions between proteins was performed in order to elucidate the mechanisms of osteogenesis and development. An interaction with a combined score > 0.4 was selected and used to construct a PPI network with Cytoscape software. Cytoscape (version 3.7.1) is an open source bioinformatics software platform for visualizing molecular interaction networks. 34,35 Dense connected regions were analysed using Cytoscape's plug-in molecular complex detection (MCODE). Our selection criteria were as follows: MCODE scores > 5, degree cut-off = 2, node score cut-off = 0.2, Max depth = 100 and k-score = 2. KEGG and GO analyses were then performed using Metascape.

| Selection and analysis of hub genes
The top 10 genes were obtained by MCC algorithm with Cytoscape's plug-in cytoHubba. Protein expression profiles of hub genes at tissue level (low, medium, high) and gene expression level (normalized expression, NX) were obtained from the HPA (Human Protein Atlas) database. The gene expression scores of hub genes at the tissue level were obtained from the Bgee database (https://bgee.org).

| Cell culture and osteogenic differentiation
Aseptic human adipose tissue was obtained from the Plastic Surgery Department of Nanfang hospital, digested with 0.15% type I collagenase for 40 minutes, terminated with growth medium and then centrifuged at 800 rpm for 5 minutes. The cells were plated in a 10-cm dish, incubated in a 37°C incubator with 5% CO 2 , and the medium was replaced 24 hours later. Cells cultured for 3-6 passages were used for seeding a 6-cm dish at a density of 8000 cells/cm 2 .
with the medium was replaced with osteogenic differentiation medium when cell confluence was ~80%. The medium was changed every 2 days. The osteogenic effects were identified by ALP and alizarin red staining at day14 and day 21, respectively.

| Immunofluorescence analysis
Cells were fixed with 4% paraformaldehyde for 10 minutes, and membranes were ruptured by 0.1% Triton X-100 for 5 minutes. The samples were blocked with 2% BSA for 1 hour at room temperature, then incubated with α-tubulin antibody (ab7291; Abcam), 4°C overnight. Incubation with the second antibody was performed at room temperature for 1 hour the next day. Nuclei were labelled with DAPI. Samples were sealed with glycerine gelatin. Images were collected using a confocal microscope (LSM 880 with Airyscan; Carl Zeiss) and analysed with software ZEN-blueedition (Carl Zeiss).

| Western blotting
Protein samples were extracted from cells using a protein extraction kit (KeyGEN Biotech), separated by sodium dodecyl sulphatepolyacrylamide gel electrophoresis (SDS-PAGE) and transferred to a polyvinylidene fluoride (PVDF) membrane (Millipore). The membranes were incubated with α-tubulin antibody (ab7291; Abcam) after blocking with 5% milk, at 4°C, overnight. Incubation with the secondary antibody was performed at room temperature for 1 hour.
FDbio-Dura ECL kit (Fudebio) was used to detect the signal. The results of Western blotting (WB) were analysed with ImageJ software and plotted with GraphPad Prism 7.0 software.

| Alizarin red staining and alkaline phosphatase staining
After fixation with 4% paraformaldehyde for 10 minutes, cells were incubated with fresh alizarin red solution (Cyagen) for 5 minutes.
Images were obtained using a Olympus microscope.
After fixation with 4% paraformaldehyde for 10 minutes, cells were incubated with alkaline phosphatase (ALP) staining liquid (Beyotime Biotechnology; Shanghai; China) for 30 minutes. Images were obtained using Olympus microscope.

| Statistical analyses
All experiments were independently repeated at least three times.
Statistical analysis was performed with GraphPad Prism 7.0 software. One-way analysis of variance (ANOVA) was used to identify significant differences, and P < .05 was considered to indicate statistical significance.

| Identification of DEGs during hBMSC osteogenesis
Because different osteogenic induction guidelines were used for the three data sets (GSE12266, GSE18043 and GSE37558 (Table S1)

| KEGG and GO enrichment analyses of DEGs
In order to analyse the biological classification of DEGs, we per-

| PPI network construction and module analysis
The PPI network of DEGs and most dense connected regions (48 nodes, 1056 edges) were obtained by Cytoscape ( Figure 3A-B).

| Selection and analysis of hub genes
Ten genes were identified as hub genes using the plug-in cy-  Table 1. According to the literature, osteogenic differentiation and adipogenic differentiation of stem cells tend to be the opposite of each other: stem cells are more likely to differentiate into adipogenic cells in an environment with lower substrate stiffness, and more likely to differentiate into osteogenic cells in an environment with greater substrate stiffness. Therefore, we compared the protein ( Figure 4A) and gene expression levels ( Figure 4B) of hub genes between human bone marrow and adipose tissue with the HPA database, and used this as a preliminary reference for identifying whether these genes were differentially expressed during osteogenesis. The results showed that, at the protein level, NUSAP1, KIF11, CCNB1 and TOP2A were highly expressed, while PBK was not detected, in bone marrow; in contrast, KIF11 was expressed at low levels, while expression of the other genes was not detected in adipose tissue ( Figure 4A). The gene expression levels of these 10 hub genes in bone marrow were all higher than in adipose tissue ( Figure 4B).
Subsequently, we compared the gene expression scores of hub genes in trabeculae bone tissue, bone marrow, subcutaneous adipose tissue and the omental fat pad using data obtained from the Bgee database. Data showed that the gene expression scores of NUSAP1, KIF11, CCNB1, CDCA8, TTK, CDC20, TOP2A, PBK and NCAPG in trabecular bone tissue and bone marrow were higher than that in subcutaneous adipose tissue and the omental fat pad.
PLK1 was the only gene whose expression score was higher in the subcutaneous adipose tissue and omental fat pad than in trabecular bone tissue and bone marrow. Therefore, we believed that the expression of these 10 hub genes might differ between bone tissue and adipose tissue, and speculated that they may represent key genes in the process of osteogenic differentiation.

| Microtubule changes during osteogenic differentiation of stem cells
Functional enrichment analysis showed that 10 hub genes were mainly concentrated in three biological processes (BP), namely sister chromatid segregation, microtubule cytoskeleton organization involved in mitosis, and chromosome condensation, as well as two cell components (CC) that is chromosome, centromeric region and spindle microtubule ( Figure 5A-C, Table 2).
However, microtubule activity is essential for each of these biological processes and cell components. Therefore, we spec- These results suggest that the microtubule-related changes are a representative biological process that occurs during the osteogenic differentiation of stem cells.

F I G U R E 3 PPI network construction and module analysis (A)
The PPI network of DEGs. The up-regulated genes are marked in red, while the down-regulated genes are marked in blue. The greater the difference in expression, the darker the colour. The size of nodes represents the difference in expression; the larger the size, the more significant the P value. B, The densest connected regions (48 nodes, 1056 edges) in the PPI network were identified with Cytoscape. C, Ten hub genes were identified in the densest connected regions with MCC algorithm, using cytoHubba. The score is indicated in red colour. Darker colour indicates a higher score However, the expression of TOP2A in osteogenesis is controversial.

TA B L E 1 Ten hub genes and their functions
Some studies have suggested that TOP2A is expressed in osteoblasts, and that parathyroid hormone can suppress the proliferation of osteoblasts partly by targeting TOP2A expression. 42 Yamagishi suggested that TOP2A plays a role in the formation of osteoclasts, 43 while Feister reported that TOP2A is not expressed in mature osteoblasts on the surface of trabeculae. 44 CCNB1 can regulate the proliferation of bone marrow stem cells 45,46 ; however, the relationship between CCNB1 and osteogenic differentiation remains poorly understood.
According to the HPA and Bgee online database, we found that the expression of these hub genes in bone marrow and trabeculae was different from that in soft tissues such as adipose tissue ( Figure 4A-C). GO analysis showed that these genes were highly and accurate chromosome segregation. [51][52][53][54] Mitosis results in the formation of two new cells with the same genetic information. 55 This involves deep remodelling of the microtubule network, including kinetochore microtubules that connect to the kinetochore of chromosomes, interpolar and central spindle microtubules that locate the furrow during cell division, and astral microtubules that anchor the furrow to the cellular cortex. 56,57 Microtubules are the core cytoskeletal filaments with a series of important cellular functions in eukaryotic cells, acting as structural scaffolds, cellular highways, force generators and signal platforms. 58 As a result, microtubule production is tightly regulated in cells. 59 The biological processes corresponding to the present hub genes involve extensive spatial, temporal and dynamic regulation of microtubules. [58][59][60] Therefore, microtubule-related changes in the osteogenic differentiation of stem cells were the focus of the present study.
To verify the universality of microtubule-related changes in the osteogenic differentiation of stem cells, we performed preliminary verification in hASCs. Western blotting results showed that microtubules showed significantly increased protein expression during osteogenesis. Immunofluorescence images indicated that microtubules underwent drastic morphological changes during osteogenic differentiation, from radial to parallel arrangement. Although more detailed analysis is necessary to confirm the observed changes in microtubule-related cell division and spindle during osteogenic differentiation, the present findings indicate that microtubule-related changes play a major role in the osteogenic differentiation of stem cells and are thus worthy of further exploration.
In summary, the purpose of this study is to identify DEGs that may be involved in the occurrence and progression of osteogenic differentiation of stem cells. A total of 240 DEGs and 10 hub genes were identified, as well as non-negligible biological changes-microtubule-related changes. These results provide the basis for further study of osteogenic differentiation of stem cells; in particular, further research should aim to clarify the biological functions of the hub genes in this context. In the future, we intend to determine which of the present hub genes are most closely related to microtubules in order to further explore the osteogenic differentiation of stem cells.

ACK N OWLED G EM ENTS
The authors thank Dr Feng Lu (Nanfang hospital, Southern Medical University, Guangzhou, China) for providing human adipose tissue for the purposes of this research.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest. Writing-review & editing (supporting).

E TH I C A L A PPROVA L
This study was approved by the Ethics Committee of the School of Basic Medical Sciences, Southern Medical University.

DATA AVA I L A B I L I T Y S TAT E M E N T
All the supporting data can be downloaded.