Bioinformatic Analysis of MicroRNAs and Signaling Networks that Orchestrates the hMSC Osteogenic Differentiation

The drive to understand how microRNAs regulate the hMSCs fate and their osteogenic differentiation has inspired stem cells biology and regenerative biology eld development. Numerous experimental data shows crucial involvement of miRs in skeletal development in embryos, osteogenic differentiation, and maturation. However, molecular mechanisms of miRs action, in other words their target signaling pathways and transcriptional factors that specic drives osteogenic differentiation is far from being understood. With meta-analysis we identied miRs signicantly involved in hMSCs osteogenic differentiation. Statistical analysis revealed signicant trend of upregulation of let-7a, mir-21, mir-26a, mir-29b, mir-101, mir-143, mir-218 during hMSCs differentiation into osteoblast. And the opposite trend was shown for mir-17, mir-31, mir-138 and mir-222: their content was signicantly lower during osteogenic differentiation. Using bioinformatics approaches we identied predictable genes-target for each miRs, analyzed signaling networks and biological process enriched by these genes. is likely involved into osteoblast homeostasis via Hedgehog signaling. Our data expands knowledge in the eld of hMSCs fate and osteogenesis orchestration by miRs, points to pro-osteogenic and anti-osteogenic miRs and their potential molecular pathways.


Introduction
Mesenchymal stem cells (MSCs) are characterized by high plasticity and their ability to self-renewal and differentiation into various types of cells. Both, bone mesenchymal stromal cells (BMSCs), as well as adipose-derived mesenchymal stem cells (ASCs), can become osteoblasts, chondrocytes, and adipocytes.
MicroRNAs (miRs/miRNAs) are important posttranscriptional regulators of gene expression. Through the RNA interference pathway, miRs regulate negatively or positively the endogenous gene expression [6,7]. It is interesting that each miRNA has hundreds of targets, and vice versa a single mRNA may be regulated by diverse miRNAs. During last decades, the crucial function of some miRNAs in osteogenesis and bone formation was shown [8,38]. Moreover, deregulation of these miRNAs has been associated with skeletal and bone malformations.
The osteogenic differentiation is a highly organized and complex process, which is orchestrated by several transcription factors including RUNX2, OSTERIX, BMP, and others [9][10][11]. The expression of these transcription factors is negatively or positively regulated by canonical Wnt, JAK/STAT, MAPK and BMP signaling pathways [10,[12][13][14][15]. However, miRs are able to regulate mentioned signaling pathways activity or osteogenic transcriptional factors expression. Indeed, miR-125 targets VDR, ErbB2 and Osterix in mice, miR-133 recognize RANX2 and this way inhibit differentiation of osteoprogenitors [16][17][18]. Diminishing of osteogenic differentiation of ASCs by miR-140-5p was shown in rats, moreover authors demonstrate that miR-140-5p effects osteogenesis via its direct targets: TLR4 and BMP2 [19]. On the contrary, promotion of osteogenic differentiation was shown for miR-27a and miR-128, where both miRs activate Wnt signaling pathway through DKK2 targeting in rats BMSCs [20,21]. Regulation of humans BMSCs differentiation into osteoblast via targeting of canonical Wnt signaling was demonstrated for miR-381 and miR-889 [22, 23] Moreover, authors found higher level of miR-889 in osteoporosis patients [23]. Thus, our understanding of miRs regulation of osteogenesis continues to expand, however role of miRs in hMSCs fate is far from understood. Thus, all published works are focusing on different miRs function and targets analysis in humans, rats, or mice mesenchymal stem cells. Considering the complex nature of miRs regulation (i.e., many targets for each known miR), identi cation of the molecular mechanisms of miRs in osteogenic signaling pathways and transcriptional factors regulation can be bene cial for the future clinical applications of MSCs.
In this study, we used publicly available experimental data. First, we conducted a meta-analysis and identi ed a list of miRNAs that signi cantly change during hMSCs differentiation into osteoblast. Using bioinformatic approaches, we predict the signaling pathways and genes targeted by identi ed miRNAs during osteogenic differentiation of human mesenchymal stem cells.

Results
First, we have analyzed all available experimental data as described in the Materials and Methods section (Fig. 1). Basing on results of this meta-analysis, we identi ed a line of microRNAs that are highly likely to be up-or downregulated during osteogenic differentiation of hMSCs (Fig. 2). Further statistical analysis revealed signi cant trend to upregulation of let-7a, mir-21, mir-26a, mir-29b, mir-101, mir-143, mir-218 during hMSCs differentiation into osteoblast. At the same time, the opposite trend was shown for mir-17, mir-31, mir-138 and mir-222: their content in such conditions was signi cantly lower.
The results of Gene Ontology (GO) enrichment analysis of identi ed miRs target genes revealed the overrepresentation of genes involved into regulation of transcription from RNA polymerase II promoter, transcription DNA templated biological processes (Fig. 3b, 4b). However, as listed in Fig. 3b and 4b identi ed miRs were enriched at the same time in more than one biological processes, so, the overlapping has taken place.
Mentioned above stem cell-speci c biological processes -stem cell differentiation, stem cell division, positive regulation of mesenchymal cell proliferation, somatic stem cell population maintenance and mesodermal cell differentiation, -were enriched with gene targets of let-7a, mir-17, mir-21, mir-29b, mir-101, mir-222 and mir-218, whereas in case of mir-21 the results of GO data analysis prove the negative effect on stem cell differentiation processes. Speci c biological process for cells division, proliferation and apoptosis are greatly enriched by gene targets of mir-17, mir-26a, mir-29b, mir-101, mir-138 and mir-222. However, analysis of experimental data proves positive regulation of cell proliferation by mir-222, whilst there was no unambiguous effect shown for mir-17. Then, data analysis con rms the hypothesis of apoptosis affection related to high mir-26a expression, but mir-138 and mir-222 are both seem to regulate the apoptotic process negatively.
Further analysis reveals speci c biological processes that correspond to bone development and morphogenesis, namely osteoblast differentiation, bone morphogenesis, bone development, endochondral ossi cation, and bone mineralization to be enriched by mir-218, mir-143, mir-101, mir-21 and mir-17 target genes. The obtained results shown mir-101 and mir-218 to affect almost all mentioned process of bone development and morphogenesis. Morphogenesis process that is extremely close to bone development also is enriched by target genes of analyzed miRs. GO enrichment analysis points that chondrocyte development and differentiation is regulated by mir-26a and mir-21, where latter has a prochondrogenic effect. Cartilage development process is enriched by mir-17 and mir-138 gene targets.
Cartilage condensation and cartilage morphogenesis are driven by mir-29b and mir-21 targets genes, respectively. Hence, it is not surprising that the majority of these miRs involves in embryonic skeletal system development and limb formation (i.e., mir-29b, mir-101, mir-218 and mir-17). Additionaly, let-7a and mir-222 take part in embryonic skeletal system development and embryonic forelimb morphogenesis regulation.
Indeed, we found that speci c biological processes such as wound healing is regulated by mir-26a, mir-17 and mir-138. Additionaly, some of miRs level are linked with regulation of macrophage activation and macrophage differentiation; cytokine and interleukin biosynthetic; in ammatory response regulation (Fig.  3b, 4b).
KEGG pathway analysis revealed 17 potential pathways signi cantly enriched by analyzed miRs target genes: TGF-beta signaling pathway, Wnt signaling pathway, mTOR signaling pathway, Jak-STAT signaling pathway, Signaling pathways regulating pluripotency of stem cells, Rap1 signaling pathway, FoxO signaling pathway, Insulin resistance and signaling, p53 signaling pathway, Cell cycle, Ras signaling pathway, Hippo signaling pathway, ErbB signaling pathway and MAPK signaling pathway (Fig. 3с, 4с). It is interesting that, cGMP-PKG signaling pathway, Calcium signaling pathway, Chemokine signaling pathway, Hedgehog signaling pathway, Platelet activation, ECM-receptor interaction pathway, VEGF signaling pathway, T cell receptor signaling pathway and B cell receptor signaling are greatly enriched by up-regulated miRs targets (Fig. 3с). Further analysis show autophagy, endocytosis, estrogen signaling pathway and osteoclast differentiation pathways were found to be targeted by drop expressed miRs during hMSC differentiation (Fig. 4с).

Discussion
Differentiation of osteogenic lineage cells is highly organized biological process which underlies in skeletal elements developing, in embryo, bone tissue remodeling and healing in adult. Many signaling networks and transcription factors as well as co-factors gathered MCSs differentiation into mature osteocyte [26]. Numerous research also pointed out that mircoRNAs has crucial function in osteogenesis and bone tissue remodeling targeting the key signaling pathways and transcriptional factors (Tabl. 1). However, different paper show that different miR or set of miRs are driving the hMSCs osteogenic differentiation. On one hand, a different biological model has been used to explain this variety of results, on the other -nature of genes expression regulation by miRs. Indeed, bioinformatic analysis revealed that a single miRNA can regulate up to several hundred target genes and vice versa, the same 3′-UTR of a gene can be targeted by many different miRs [30,31]. In our opinion, it is depending on miRs or its targets tissue speci c expression, hence on cells type. Here we provide meta-analysis of experimental data and identify the line of miRs which are signi cantly involved in hMCSs osteogenic differentiation. The results of our analysis demonstrate complex involvement of different miRs in osteogenic hMSCs differentiation and show that each step of osteoblasts differentiation is orchestrated by a few miRs via similar molecular mechanisms (Fig. 5). Vice versa, the same miRs regulate different signaling pathways, regulating osteogenesis.
With KEGG pathway and GO analysis we found that let-7a, mir-17, mir-21 and mir-101 in hMSC are likely to be involved in transforming growth factor β (TGF-β) signaling where mir-17 is predicted to inhibit the TGF-β cytoplasmic mediator activity. At the same time, mir-29b, let-7a and mir-101 were predicted to regulate the platelet-derived growth factor (PDGF). Aadditionally, mir-101 was predicted to target the broblasts growth factor (FGF). Mir-29b was predicted to activate the vascular endothelial growth factor (VEGF) and maintains the cellular response to TGF-β and FGF. Together up-regulated PDGF, TGF-β, FGF and VEGF create appropriate microenvironment for hMSCs to enter in osteogenic lineage differentiation, namely promote collagen synthesis and bril organization. GO analysis additionally shows important regulatory role of let-7a and mir-29b in collagen synthesis and bril organization. Having summarized our data, we suppose that balance between high level of let-7a, mir-21, mir-29b, mir-101 and low level of mir-17 is necessary for microenvironment formation and hMSCs turn on into osteogenic lineage (Fig. 5).
Widely known TGF-β/BMP activates SMAD-dependent and non-SMAD-dependent signaling which results in main osteogenic transcriptional factor Runx2 expression [12]. TGF-β-SMAD signaling raises MSCs proliferation, chemotaxis, and early differentiation of osteoprogenitors, where BMP-SMAD signaling regulate almost every stage of osteoblast differentiation and maturation [12]. During osteoprogenitors development BMP-SMAD also activates the Runx2 and promotes hMSCs differentiation. We suggest that let-7a is mostly involves in TGF-β modulation, where mir-26a leads to all steps of osteoblast development and maturation via BMP-SMAD targeting. Additionaly, TGF-β/BMP-SMAD dependent pathway is supposed to be regulated by lowering in mir-17 which targets SMAD7, SMAD6 and SMAD1. But, surprisingly, overexpressed miRs mir-21 and mir-101 also targets SMAD proteins (SMAD7 and SMAD6) inhibiting SMAD dependent transcriptions.
However, non-SMAD-dependent pathway (JNK, ERK, p38 MAPK, Pi3K-Akt signaling) also promote osteogenesis and main transcriptional factors (Runx2, Osx, Dlx5 and other) expression [12][13][14][15]26]. Mentioned signaling can activate the Runx2 expression in hMSCs. Moreover, all these signaling together with mentioned above BMP pathway regulates the Osterix expression Runx2 independently. With DAVID Bioinformatics Resources we found that non-SMAD-dependent pathway activity is orchestrated by interplaying between different up-and downregulated miRs. As we mentioned above let-7a targets TGFβ-SMAD signaling and is crucial for entering of hMSCs into osteogenesis, however this microRNA regulates the TGF-β-non-SMAD-dependent signaling including MAPK, PI3K-Akt, cGMP-PKG, Jak-STAT, Hippo and AMPK positively affected their activity. Our data showed that let-7a is versatile regulator of osteogenesis and high level of its expression maintain osteoprogenitors differentiation as well as osteoblast maturation and metabolism (via AMPK pathway).
Our data show that not only let-7a, but MAPK, ErbB/MAPK and Pi3K-Akt signaling are predictably regulated by mir-29b, mir-101, mir-143 and mir-218 that are also expressed at higher level during osteogenic differentiation. Additionaly MAPK and Pi3K-Akt signaling supposed to be regulated by let-7a and mir-26a via cGMP-PKG signaling pathway and mir-21, mir-101, mir-29b and mir-218 through Ras signaling pathway. However, GO analysis revealed that mir-21 can affect negatively the ERK1/1 and Ras signaling pathway via inhibiting of MAPK activity. As shown, Ras/MAPK signaling branch activity is associated with osteogenesis inhibiting [32], so mir-21 may inactivate Ras/MAPK signaling promoting in such way osteogenic differentiation. We found that mir-101 and mir-218 also are possible inactivators of MAPK-cascade and Jak-STAT signaling. Thus, our research indicates that high content of let-7a, mir-21, mir-26a, mir-29b, mir-143 can be speci cally involved in non-SMAD dependent signaling and promote osteogenic TF expression (Fig. 5). On the other side non-SMAD dependent signaling modulated by miRs shows signi cantly lower expression in osteogenesis. Namely, mir-17, mir-222, mir-138 and mir-31 are supposed to be involved in non-SMAD-dependent pathway (i.e., MAPK, Pi3K-Akt, and canonical Wnt signaling) regulation. However, GO analysis indicates that mir-17 positively regulates Wnt and p38-MAPK signaling cascades where inhibits TGF-β signaling. Hence, we suppose that mir-17 level diminishing is necessary for TGF-β/non-SMAD dependent signaling activity. Our data show that mir-31 also positively regulates MAPK-cascade, however low level of its content points that this microRNA does not require for MAPK-dependent expression of osteogenic TF.
GO analysis reviled that mir-222 and mir-138 can modulate canonical Wnt via noncanonical calcium and Hippo signaling, respectively. Moreover mir-138, based on KEGG pathway analysis, is supposed to negatively regulate canonical Wnt-dependent transcription through Notch signaling pathway. Hence, the lower level of mir-222 and mir-138 in hMSCs maintains the canonical Wnt signaling up-regulation, where the latter one regulates other osteogenic TF expression (Runx2, Osx, Dlx5 and ATF4), strongly stimulates osteoblastogenesis and inhibits adipogenic TF (C/ERB and PPARγ) [10,13]. Moreover, canonical Wnt involves in osteoblast maturation and metabolism [33]. GO analysis revealed that canonical Wnt activity can be positively regulated by mir-26a, mir-101 and mir-29b which expresses on higher level during osteogenesis, where mir-101 and mir-29b regulate β-catenin binding and β-catenin-TCF complex assembly. Summarizing these data we came to the conclusion that canonical Wnt signaling activity in osteogenesis depends on low content of mir-222, mir-17 and mir-138 and high expression of mir-26a, mir-101 and mir-29b.
Other crucial for osteogenesis signaling is Hippo/Yap driven by a few miRs. Notably, Yap is a co-regulator for main osteogenic TF and act via both SMAD-dependent and non-SMAD-dependent molecular mechanisms, so it is up-regulation is crucial for osteogenic lineage commitment as well as osteoblasts maturation [34,35]. GO and KEGG pathway analysis revealed that Yap dependent transcription can be modulated by up-regulated let-7a, mir-101 and downregulated mir-31, mir-138 and mir-222. Our nding points that mir-138 and mir-222 can regulate Yap/TCF/Lef dependent transcription. Mir-31 and mir-222 are supposed to inhibit Yap/Taz and prevent its transcription activity, where mir-17 modulate TGF-β-SMAD/Yap, BMP-SMAD/Yap and Yap/TCF/Lef dependent signaling. Summarizing we assume that lower level of mir-31, miR-138, mir-222 are important for Yap dependent transcription of main osteogenic TF in hMSC during all steps of osteogenesis (Fig. 5).
Osteoblast differentiation relates to metabolic changes, hence AMPK signaling pathway is one of crucial factors which controls different aspects of metabolism in cells including osteoprogenitors and osteoblast. Moreover, AMPK signaling is required for hMSC to enter the adipogenic and osteogenic lineages through Runx2 and PPARγ and to be able to maintain osteogenesis and osteoblast homeostasis [36]. With bioinformatic analysis we found that AMPK signaling during osteogenesis can be positively modulated by let-7a, mir-29b and mir-218 which are signi cantly up-regulated (Fig. 5).
Other important signaling involved in osteoblast homeostasis as well as differentiation is Hedgehog signaling pathway, which play dual role in bone homeostasis and required in both bone resorption and formation process. Additionally, Hedgehog signaling induce osteoblast differentiation in BM-MSC through Runx2 and Osx expression [37]. We found that during osteogenesis mir-101 can promote Hedgehog signaling activity.
Our data expands knowledge in the eld of hMSCs fate and osteogenesis orchestration by miRs, demostrates pro-osteogenic and anti-osteogenic miRs and their potential molecular pathways.

Search strategy and selection criteria
We carried out a comprehensive literature search in the PubMed, Medline and Google Scholar databases using the following keywords: miR, microRNA, osteogenesis, osteogenic differentiation, MSC. The articles were not limited by the publication date, but their language was restricted to English. A total of 632 papers were identi ed with the initial search. The inclusion criterion for article selection was the presence of original research data in the eld of miRNA expression in osteogenesis or osteogenic differentiation. The exclusion criteria were: 1) the paper did not contain original research data; 2) the paper did not contain numeric data or any formalized miRNA expression data; 3) data provided in the paper could not be compared with other research data.
As a result, 87 papers were included in this meta-analysis. (Fig. 1). The detailed information of relevant citations is listed in Supplement 1.

Data extraction and study assessment
Two researchers independently extracted data and reviewed the essence of the articles to determine whether they met the criteria for inclusion. All disagreements were discussed. Gathered data were analyzed by the following criteria: 1) miRNA expression has been measured; 2) process has been studied. In total 2480 different miRNA have been analyzed in included papers, but only miRNA mentioned by 3 different authors have been studied.

Statistical analysis
The Aggregate Data approach was used. First, we calculated the Pearson correlation coe cient using the standard procedures. Then we performed one-way ANOVA for the study of changes in levels of individual miRNAs during processes of osteogenesis, osteogenic differentiation, and in osteopathological conditions. For statistical analysis of the received data, the simple regression method was used.
We did not employ data weighting methodology, as data used for the purposes of this study did not include any clinical trial data and all experimental groups in the studies were rather small.
The effect size was calculated by Hedges and Olkin.
Both types of publications -those that demonstrated a signi cant changes of miRNA expression and/or experimental modulation of processes have been studied, and those that did not -were included into analysis; the fail-safe number was calculated by Rosenthal method.
All calculations were performed with Origin 8.1.

Declarations
Author Contributions: LLM, and OOP planned work, discussed the results, and wrote the manuscript. DSA, LLM extracted and analyzed data. LLM performed the statistical analysis. Bioinformatic analysis was done by OOP. All authors participated in the design of the study, interpretation, and analysis of the data, reviewed of the manuscript.

Funding:
This work was funded by the Ukrainian State Budget Program "Support for the Development of Priority Areas of Scienti c Research" (Code: 6541230).

Availability of data and materials
The datasets used and analyzed during the current study are available from the corresponding author on request.

Con icts of Interest:
The author declares no con ict of interest. The PRISMA ow diagram  Bioinformatic analysis of miRs downregulated in osteogenic differentiation of hMSCs. a) Venn diagram of the identi ed miRs target genes in Targetscan, miRDB, and miRanda databases; b) Gene Ontology of the target genes; c) KEGG pathway of the target genes.