Identication of Biomarkers Modulating Osteogenic Differentiation With Immune and Inammation Landscapes in Mesenchymal Stem Cells: A Bioinformatics-based Comprehensive Study

Background: Inducing osteogenic mesenchymal stem cell (MSC) differentiation brings a huge potential for clinical bone defect repair. This present study is aimed to identify novel key biomarkers for osteogenesis of MSCs, and analyze their possible regulation roles on immune and inammation during this process. Results: Seven datasets (GSE159137, GSE159138, GSE114117, GSE88865, GSE153829, GSE63754, GSE73087) were obtained from the gene expression omnibus database, and 200 differentially expressed genes (DEGs) were identied from the train datasets. Biomarkers were screened through Logical regression of the selection operator algorithm, whose expressions and diagnostic performances were veried in the test datasets. FKBP prolyl isomerase 5 (FKBP5), insulin like growth factor binding protein 2, prostaglandin E receptor 2, SAM and HD domain containing deoxynucleoside triphosphate triphosphohydrolase 1 (SAMHD1) and transmembrane O-mannosyltransferase targeting cadherins 1 were established to be osteogenesis-related biomarkers. Finally, the DEGs of immune and inammation in osteogenesis were detected, and their correlation with the hub biomarkers were analyzed. The ve selected biomarkers are strongly associated with immune and inammation related genes. Conclusion: Taken together, our manuscript identied ve biomarkers related to stem cells osteogenesis and discussed their potential perspectives in immunoregulation, inammation and tissue engineering during the process.


Background
Though omnipresent challenges exist in face of the clinical and surgical repairing of large-volume skeletal abnormalities, recent stem cell-based therapy (SCBT) has shown great potential for the treatment of severe bone defects (1,2). In biological conditions, healthy bone tissues serve limited capabilities in self-renewal and regeneration which is majorly controlled by the activities of osteoclasts and osteoblasts. Nevertheless, such properties can be signi cantly suppressed, even halted in case of severe skeletal trauma or defects, resulting in great setbacks when taken into clinical consideration. Mesenchymal stem cells (MSCs), also known as mesenchymal stromal cells, including adipose-derived stem cells (ASCs) and bone-marrow stromal cells (BMSCs), plays a major role in regenerative medicine due to their fundamental properties in cell-regeneration and differentiation. They have high osteoinductive and osteogenic potential, which is theoretically ideal for bone regenerations and shows a promising alternative (3). Therefore, the search for and identi cation of important biochemical molecules and genes that may initiate or control seed osteogenic differentiation is vital for the establishment of relevant defect treatment strategies.
Depending on speci c microenvironment, MSCs which have various properties in regulating immune and in ammation play crucial biological roles in ways of direct immune cell contact and the production of in ammatory regulatory molecules, namely, the paracrine effects (4,5). Though the effect of "sensor and switcher of the immune system", MSCs have the capabilities of both up-and down-regulating in ammatory processes, mainly by responding to danger indicators and releasing anti-in ammatory mediators through activated Toll-like receptors (TLRs) ligands, such as TLR9 TLR7 and TLR2, and expressing the immune suppressors PD-L1 as well as PD-L2 under the stimulation of interferon gamma (IFN-γ), respectively(6). MSCs also has various effects in coordinating the migration, proliferation, and activation of immune cells according to different stages of osteogenesis and different categories of immune cells (7). And consistent with changes in proin ammatory and anti-in ammatory cytokines presenting in the micro-environment, MSCs secrete cytokines including prostaglandin E2 (PGE2), transforming growth factor beta (TGF-β), histocompatibility locus antigen-G (HLA-G) that induce the formation of regulatory T cells, suppress neutrophil migration and more importantly, serve the role in immune-inhibition by monocyte/macrophage regulation(8-10). Therefore, exploring the potential roles of immunomodulation and in ammatory regulation which involved MSC osteogenic differentiation is crucial in discovering novel therapeutic strategies and improving currently ine cient bone remodeling.
In this study, we performed a comprehensive strategy to determine biomarkers related to osteogenic differentiation of MSCs, and further investigated the relationship between those biomarkers and genes involving immunoregulation and in ammation participating in this process. The objective of our research is to form a basis for further explorations of signi cant genes, biochemical pathways as well as vital bioactive molecules affecting the osteogenic role of mesenchymal stem cells.

Data Collection and Processing
The study owchart is presented in Figure 1. Three datasets which expression pro ling by high throughput sequencing (GSE159137, GSE159138, GSE114117) and four datasets which expression pro ling by array (GSE88865, GSE153829, GSE63754, GSE73087) were retrieved from the GEO database (http://www.ncbi.nlm.nih.gov/geo), and the characteristics of these included datasets were presented in Table 1. The three high throughput sequencing were selected as the train group containing a total of 7 uninduced MSCs samples and 12 osteogenic induced MSCs samples. And four datasets which expression pro ling by array, with a total of 10 uninduced MSCs samples and 18 osteogenic MSCs induced samples, were utilized for test validation. The RNA algorithm was applied to background correction and data normalization (39). Then, probes were transmuted into their corresponding gene symbols based on platform annotation information.

Enrichment Analyses of DEGs
Gene Ontology (GO) (http://geneontology.org/) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp/) were used for functional and pathway enrichment analyses, respectively. Disease Ontology (DO) (http://disease-ontology.org) comprises 8043 inherited, acquired, and developmental human diseases. Based on molecular signature database, Gene Set Enrichment Analysis (GSEA) (http://software.broadinstitute.org/gsea/index.jsp) is a statistical tool for interpreting gene expression data. GSEA 4.1.0 was utilized to predict the potential functions and downstream access of the two clusters. Prior to enrichment analyses, "org.Hs.eg.db" R package (http://www.bioconductor.org/packages/release/data/annotation/html/org.Hs.eg.db.html), a human genome annotation database was used for conversion of gene symbol codes into Entrez ID. FDR q-value < 0.05 were considered markedly enriched for GO terms, DO and GSEA analyses, while p-value < 0.05 was considered signi cant for KEGG analyses.
Screening for MSCs osteogenesis related biomarkers LASSO was used to investigate osteogenesis-related biomarkers via R package "glmnet" R package (https://mirrors.sjtug.sjtu.edu.cn/cran/web/packages/glmnetSE/index.html) in train group. while "pROC" R package (https://cran.rstudio.com/web/packages/pROC/index.html) was used for computing the area under the curve (AUC) in train and test group. Accuracy, sensitivity, and speci city of a receiver operating characteristic (ROC) curve to assess the diagnostic ability of the markers. p < 0.05 was the cut-off for signi cance.

Selection of immune and in ammation related genes
To make the research more comprehensive, we selected immune and in ammation related genes, and the differential expressed genes in osteogenesis and non-osteogenesis groups was visualized by using the methods mentioned above. PPI networks were established to assess the associations among genes in selected modules using the Search Tool for the Retrieval of Interacting Genes version 11 (STRING V11, https://string-preview.org/) with > 0.4 as the con dence level. Cytoscape version 3.8.2 was used for network visualization. Hub genes, which are presented by highly interconnected nodes, may play vital roles in the PPI network.

Interaction analysis of selected biomarkers
Co-expression analysis of selected biomarkers with immune-related genes and in ammation-related genes was performed using "limma" in R (http://www.bioconductor.org/packages/release/bioc/html/limma.html ), respectively. We then performed co-expression analysis, and the following parameters were used as lter conditions to select related genes: "correlation coe cient = 0.4" and "pvalueFilter = 0.05".

Identi cation of DEGs
Differential expression analysis was conducted based on the train datasets to screen for DEGs. There were clear differences between osteogenesis MSCs and the control ones. From our heatmap ( Figure 2A) and volcano plots ( Figure 2B), we observed there were 200 DEGs. Among them, 156 genes were upregulated while 44 genes were downregulated.

Functional enrichment analyses of DEGs
With regards to biological processes, GO analysis showed that robust DEGs were mainly enriched in cellular divalent inorganic cation homeostasis, epithelial cell proliferation, and cellular calcium ion homeostasis in the biological process (BP) part. The cellular component (CC) part, the DEGs mainly mainly concentrated in collagen-containing extracellular matrix, postsynaptic density, and asymmetric synapse. In Molecular functions included peptide binding and amide binding among other important functions ( Figure 3A). The KEGG analysis demonstrated that neuroactive ligand-receptor interaction, complement and coagulation cascades, as well as cytokine-cytokine receptor interaction were the most signi cant processes in MSCs osteogenic differentiation process ( Figure 3B). DO analysis showed that DEGs were enriched in kidney disease and urinary system disease ( Figure 4A).
Moreover, the GSEA analysis result indicated active pathways that varied between uninduced and osteogenic induced groups ( Figure 4B), including cellular ion homeostasis, divalent inorganic cation homeostasis, in ammatory response, and metal ion homeostasis were established to be active in the osteogenic induced group. While the active pathways of uninduced group ( Figure 4C) were chromosome segregation, polymeric cytoskeletal ber, postsynapse, postsynaptic membrane, and synaptic membrane.
To assess the accuracy for the biomarkers, we used 4 datasets (GSE88865, GSE153829, GSE63754, GSE73087) as veri cation sets to validate the e cacies of candidate biomarkers. Differential expressions of these ve genes, including FKBP5, IGFBP2, PTGER2, SAMHD1 and TMTC1, were veri ed in test group, and their expressions were all signi cantly higher in the osteogenic induced groups ( Figure 5B-H). The ve biomarkers were all with an AUC more than 0.7 in the test datasets, and SAMHD1 presented the highest AUC of 0.861, indicating the capabilities of the ve biomarkers to diagnose MSCs osteogenic differentiation with excellent speci city and sensitivity ( Figure 6).

Identi cation of DEGs related to immune and in ammation
Heatmaps revealed 38 differentially expressing immune-related regulators between osteogenic induced and control groups with FDR q-value < 0.05 ( Figure 7A), with 31 upregulated and 7 downregulated. The interrelation of immune-related regulators was constructed using the STRING database, the proteinprotein interaction (PPI) network analysis ( Figure 7B) indicated that IL6, LEP and ANGPT1 had the most interactions of the studied genes. GO and KEGG analyses were performed to the potential functions and pathways of the differentially expressing immune-related genes (Figure 8). The GO results revealed that the differentially expressing immune-related regulators were majorly concentrated in epithelial cell proliferation in BP, basolateral plasma membrane, receptor ligand activity in CC, receptor ligand activity and signaling receptor activator activity in MF. In KEGG, it concentrated in neuroactive ligand-receptor interaction, as well as cytokine-cytokine receptor interaction.
We screened differentially expressing in ammation-related regulators between osteogenic induced and control groups through p-value < 0.05 and 32 upregulated and 18 downregulated genes were presented in Heatmaps ( Figure 9A). The interaction as well as the expression number of in ammation-related regulators were shown in Figures 9B, CXCL8 and IL6 were treated as the central gene. GO and KEGG enrichment analyses found that response to molecule of bacterial origin and response to lipopolysaccharides in BP, secretory granule membrane and plasma membrane signaling receptor complex in CC, receptor ligand activity and signaling receptor activator activity in MF, while cytokinecytokine receptor interaction were the most important functions and pathways of the differentially expressing in ammation-related regulators according to the enrichment scores in KEGG ( Figure 10).

Potential roles of hub biomarkers in immune and in ammation
We analyzed RNA-seq transcriptome to explore the association between differentially expressed immunerelated genes and selected biomarkers ( Figure 11A) as well as the relationship between differentially expressed in ammation-related genes and the ve hub biomarkers ( Figure 11B). There seemed to be a regulating network of FKBP5, IGFBP2, PTGER2, SAMHD1, and TMTC1 with immune-related and in ammation-related genes. The correlation coe cients and positive and negative regulatory relationships between the genes we predicted are shown in the Table 2 and Table 3. In the process of stem cell osteogenesis, the biomarkers we screened showed a strong correlation with immunity and in ammation, which laid a foundation for us to further study the mechanism of stem cell osteogenesis.

Discussion
Although the signi cance of MSCs in bone tissue engineering have been evaluated, their e cacy in osteogenic differentiation is limited without extra interventions in vivo and in vitro (11,12). How to improve the e ciency of regenerative abilities of MSCs in damage restoration is still urgent. The search of novel key biomarkers in stem cell osteogenesis will be helpful in clarifying osteogenic potentials and mechanisms of stem cells, thus designing more appropriate bone lling materials to induce bone regeneration with increased pro ciency and biocompatibility.
In this study, we identi ed hub biomarkers of stem cells osteogenesis with a comprehensive strategy of LASSO logistic regression and ROC curve analysis. The ve hub biomarkers we identi ed are FKBP5, IGFBP2, PTGER2, SAMHD1 and TMTC1. Although there have been studies focusing on exploring stemcell-related osteogenic biomarkers via bioinformatics analyses (13,14), the hub markers SAMHD1 and FKBP5 we determined in this study have an average AUC over 0.8 in test datasets, indicating their reliable predictive abilities.
Among these 5 hub biomarkers, the effects of IGFBP2 and PTGER2 on osteogenesis have been reported, but remain controversial. IGFBP2, an important member of the insulin-like growth factor family, plays crucial roles in growth, development, and metabolism of human body (15,16). HamidoucheI et al. demonstrated that IGFBP2 has the property of promoting the osteogenic differentiation of MSCs by enhancing the osteoblast phenotype genes, including RUNX2, ALP and COI1A1B (17). And consistent with our results, the expression of IGFBP2 increased signi cantly when inducing osteogenesis(18). Moreover, a recent study revealed that IGFBP2 is also capable of activating the JNK/Akt pathway and inducing the adipogenic differentiation of MSCs (19). Further researches are needed to elucidate the mechanisms and speci c roles played by IGFBP2 in promoting osteogenesis or adipogenesis in MSCs. PTGER2 acts as a main receptor of prostaglandin E2, which has been shown to encourage osteogenic differentiation of MSCs, and stimulate the resorption and formation of bones (20)(21)(22). It is suggested in several current studies that local application of PTGER2 agonists to fractured rats in vitro can effectively stimulate the formation of bone marrow and periosteum, thus promoting fracture repairing (23)(24)(25). However, the direct effect of PTGER2 on the osteogenic differentiation of MSCs needs further research to con rm. TMTC1, a glycosyltransferase, promotes O-mannosylation of cadherin and maintains intracellular calcium homeostasis(26-28). Glycosylation is a key posttranslational modi cation process that participates in controlling extracellular matrix formation during osteogenesis (29,30). GO analysis revealed that DEGs were mainly concentrated in extracellular matrix organization, and collagencontaining extracellular matrix. Previous studies have revealed that extracellular matrix modulates osteogenic effects by adhering to cytoskeletal proteins in MSCs (31,32). Together with our results, TMTC1 and glycosylation may play crucial roles in the extracellular matrix formation during osteogenic differentiation.
The strong relationship between stem cells osteogenesis and immunity has been well reported within the osteoimmunology scope, and a controlled in ammatory response is vital for favoring osteogenesis (33)(34)(35). In this study, we identi ed differentially expressing immune and in ammatory related regulators between osteogenic induced and uninduced groups, and numerous DEGs between groups were found, demonstrating a signi cant difference of immune and in ammation modulation between osteogenic and nonosteogenic status. And we further analyzed their association with hub biomarkers. Among the ve biomarkers, SAMHD1 and FKBP5 represented to be most correlated with immune-related and in ammation-related genes in osteogenic differentiation. SAMHD1 is a deoxynucleoside triphosphohydrolase that can suppress innate immune responses to viral infection through interacting with various key proteins in immune signaling pathways(36, 37). FKBP5 is a stress-responsive molecule which also possesses the property of modulating immune function. Earlier researches have demonstrated the effects of SAMHD1 and FKBP5 on immune response and in ammation, but whether these effects present during the process of osteogenesis have not been fully elucidated(38). Our present results suggest that the potential mechanisms of SAMHD1 and FKBP5 exerting in osteogenesis might be greatly linked to immune and in ammation regulation.
In summary, we combined bioinformatic analysis and ROC analysis to screen the biomarkers associated with osteogenesis of stem cells, and analyzed their potential roles with immune and in ammation during this process. Although these hub genes are highly associated with osteogenesis of stem cells based on this study and other previous studies mentioned above, their biological mechanisms in osteogenesis remain unclear. Therefore, due further researches can be focused on specifying the precise mechanisms as well as discovering novel, related signaling pathways and bioactive molecules.

Conclusion
In conclusion, by combining bioinformatic and analysis, the current study has screened several hub biomarkers which are related to the osteogenic differentiation of mesenchymal stem cells, and analyzed their potential roles in immunomodulation and in ammatory responses during this process. It is suggested, by our results, the possible properties of these markers act through signaling pathways that are related to their effects in immune and in ammation regulation, also lay a foundation for us to further study the mechanism of mesenchymal stem cell osteogenesis.

Declarations
Ethics approval and consent to participate Not applicable

Consent for publication
All authors have provided their consent for publication.

Availability of data and materials
The datasets generated and/or analyzed included in the manuscript are available in the GEO repository, http://www.ncbi.nlm.nih.gov/geo.  Figure 5 See image above for gure legend.

Figure 6
Page 23/23 See image above for gure legend.   See image above for gure legend.  See image above for gure legend.