Integrated Analysis of Hub Genes and miRNAs in Osteoporosis

Background and objective: Osteoporosis (OP) is a systemic disease of bone metabolism, characterized by decreasing bone mass, increasing bone microstructure damages and fracture risk. It affects the quality of life of nearly 200 million people worldwide and is a major burden on the public health systems. We want to identify hub genes and miRNAs by the miRNA-mRNA interaction network in osteoporosis disease so that further understand the pathogenesis of this disease. Materials and methods: The differentially expressed miRNAs (DEMis) and mRNAs (DEMs) were selected using data of OP patients downloaded from the GEO database (GSE93883, GSE74209 and GSE35959). Gene Ontology (GO) pathway analysis and transcription factor enrichment analysis were used for selecting DEMis, and the target mRNAs of DEMis were ltered by using FunRich. Cytoscape software was used to visualize the network between miRNAs and mRNAs and calculate the hub genes. GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were used to analyze the mRNAs in the regulatory network. Results: A total of 17 DEMis and 655 DEMs were selected, from which we reconstructed the miRNA-mRNA network consisting of 6 miRNAs and 37 mRNAs. The top 10 nodes, hsa-let-7a-5p, hsa-miR-92a-3p, hsa-miR-92b-3p, hsa-miR-223-3p, hsa-miR-320c, SLC2A3(Solute Carrier Family 2 Member 3), LBX1(ladybird homeobox 1), HCN2(hyperpolarization-activated cyclic nucleotide-gated ion channel 2), DAB2IP(DAB2 Interacting Protein) and CIC(capicua transcriptional repressor), were identied as important regulators. Conclusions: The study uncovered several important hub genes and miRNAs involved in the pathogenesis of OP, among which, the hsa-let-7a-5p, hsa-miR-92a-3p and hsa-miR-92b-3p may play an important role in osteoporosis, which can help us provide potential therapeutic


Introuduction
Osteoporosis(OP) is a systemic disease of bone metabolism, characterized by decreasing bone mass, increasing bone microstructure damages and fracture risk [1]. Bone remodeling is maintained by osteoblastic regulation of bone formation and osteoclast mediated bone resorption. When bone resorption is greater than bone formation, it will lead to bone loss and lead to osteoporosis [2]. Osteoporosis affects the quality of life of nearly 200 million people worldwide and is a major burden on the public health systems [3]. Despite some advance in the treatment and diagnosis, the morbidity of osteoporotic fractures was not signi cantly reduced.
Besides, many studies have revealed the pathogenesis of osteoporosis at the level of cell and animal models, but few have revealed the pathogenesis of osteoporosis at the level of gene interaction.
The interactions between microRNAs (miRNAs) and mRNAs are involved in almost all biological processes, which play a critical role in the progression of various diseases. MiRNA research is expanding because miRNAs are crucial regulators of gene expression and promising candidates for biomarker development. MiRNA mimics and miRNA inhibitors currently in preclinical development have shown promise as novel therapeutic agents. [4,5] Previous studies have identi ed miRNA-182 and miRNA-365a-3p as regulators that are associated with osteoporosis [6,7]. Bioinformatics is an interdisciplinary eld that combines computer science and biology to study, analyze and interpret large sets of biological data. Computational tools are routinely used to model biological processes, predict disease mechanisms and generate experimental hypotheses [8]. Bioinformatics tools are also widely used to screen potential genes related to several human diseases, which are then validated by comprehensive follow-up experimental studies [9].
In the present study, bioinformatics tools were used to analyze the OP expression pro le chips in a public gene chip database, which provided a theoretical basis for the biological functions of related genes and their molecular mechanisms that were involved in the occurrence and development of OP. We reconstructed the miRNA-mRNA network according to the miRNA sponge theory [10]. The study distinguished human OP related miRNAs with high credibility and provided a novel approach to identify pathological mechanisms and potential targets for OP.

Data Download and Screening Strategy
Microarray datasets of OP were downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geopro les). The miRNA expression pro le data GSE93883 performed on Platform GPL18058 contained plasma from six healthy people and twelve osteoporosis patients. And the miRNA expression pro le data GSE74209 [11] performed on Platform GPL20999 contained bone tissue from six healthy people and six osteoporosis patients. The gene/mRNA expression pro le GSE35959 [12] was performed on Platform GPL570, with human bone mesenchymal stem cells from ve OP patients and fourteen healthy people.
We used the limma package in R language to get the differentially expressed miRNAs (DEMis) and mRNA (DEM) between the OP groups and healthy controls. Bayesian methods corrected the batch effect. If multiple probes are located within the same gene, the average expression value used is equal to the expression value of the gene. The t-test was used to screen the differentially expressed genes. The DEMis were screened by the P values < 0.05 and log|FC|>1. And the DEMs were screened by the P values < 0.05 and log|FC|>1. To show the differential expression of DEMis and DEMs in different samples, volcano maps and heat maps were drawn by applying the plot and pheatmap packages in the R language. Then, the overlapped DEMis of the GSE93883 and GSE74209 were identi ed with Venn diagram analysis.

GO Enrichment Analysis for the Targets of Transcription Factors
DEMis had uploaded FunRich (3.1.3), which is a common tool for GO functional enrichment analysis for transcription factor pathway enrichment targets (gene/mRNA) [13]. GO enrichment analysis was also used to develop the interaction network analysis between miRNAs, gene/mRNA, and transcription factors (already integrated the explored information of miRNA and potentially targeted gene/mRNAs) [14].

Prediction of the Targeted mRNAs of DEMis
The miRNA enrichment function of FunRich could be applied to perform miRNA enrichment analysis, to predict targets of miRNAs, or to nd miRNAs through given target genes. Functional analysis of DEMis target genes was conducted to predict target genes with FunRich [15]. Subsequently, the predicted mRNAs of overlapped DEMis were further ltered by matching the DEM previously selected, and then the DEMi-DEM pairs were obtained.

Construction of the miRNA-mRNA Regulatory Network
The miRNA-mRNA network was constructed by putting all the DEMi-DEM pairs selected above together, and Cytoscape software [16] (version 3.8.0) was applied to visualize it simultaneously. All the node degrees, closeness, and betweenness of the regulatory network were calculated at the same time.

GO and KEGG Enrichment Analyses on mRNAs in the Network
We used enrichplot and ggplot2 packages in the R language to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses on mRNAs in the network. P value < 0.05 was considered statistically signi cant.

Results
The overlapped DEMi Screening Results in OP According to the ltering criterion described before, the cutoff value of log FC of miRNAs was 1. We identi ed 90 and 137 DEMis in GSE93883 and GSE74209, respectively. The distribution of differential miRNA expressions (GSE93883 and GSE74209) between OP and healthy controls was visually illustrated by the volcano map on the correlation of − log10 (P value) and log (FC) ( Fig. 1(a)). The heatmap was also painted to reveal the differences between OP and healthy groups ( Fig. 1(b)). And then we obtained 17 overlapped DEMis between the two datesets identi ed with Venn diagram analysis. (Fig. 1(c)) The 17 overlapped DEMis, including 7 upregulated miRNAs and 10 downregulated miRNAs, are presented in Table 1. A total of 1959 genes were mapped into 200 transcription factors. By exploring the enrichment of targets of transcription factors, we screened the top 10 transcription factors which had strong closeness to miRNAs, including SP1, EGR1, SP4, POU2F1, NKX6-1, LHX3, MEF2A, FOXA1, RORA, HOXA7 ( Fig. 2(a)), indicating these transcription factors are in regulatory relationships with DEMis.
As to the molecular function (MF) terms by GO enrichment analysis, most of the genes were involved in Transcription factor activity, Protein serine/threonine kinase activity, Ubiquitin-speci c protease activity, Receptor signaling complex scaffold activity, Auxiliary transport protein activity. (Fig. 2(b)). GO enrichment analysis also revealed that the top 5 biological progress (BP) terms with the most enriched targets of DEMis, including Regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolism, Signal transduction, Cell communication, Regulation of translation, Transport. (Fig. 2(c)). The top 5 of the cellular component (CC) terms with the most enriched targets of DEMis were proved to be the Nucleus, Cytoplasm, Lysosome, Cytosol, Cyclin-dependent protein kinase holoenzyme complex. (Fig. 2(d)).

Results of the DEM Screening in OP
In this study, the expression levels of mRNAs from 12 OP patients and 6 healthy controls were explored.  Reconstruction of the miRNA-mRNA Network in OP As shown in Fig. 4, a miRNA-mRNA regulatory network with 37 mRNAs and 6 miRNAs was constructed to further demonstrate the interaction between overlapped DEMis and DEMs, which is helpful to understand the role of miRNAs in OP. The parameters of degree, closeness, and betweenness in the network were calculated by the plugin cytoHubba in Cytoscape (version 3.8.0). The top 10 nodes, including hsa-let-7a-5p, hsa-miR-92a-3p, hsa-miR-92b-3p, hsa-miR-223-3p, hsa-miR-320c, SLC2A3, LBX1, HCN2, DAB2IP and CIC, could be selected as hub nodes (Table 3). Three miRNAs (hsa-let-7a-5p, hsa-miR-92a-3p, hsa-miR-92b-3p) were considered to have the most node degrees, indicating which might play critical roles in the genesis and development of OP as the key miRNAs. In the top 10 nodes, SLC2A3, LBX1, HCN2, DAB2IP, CIC were the target mRNAs of hsa-miR-92a-3p and hsa-miR-92b-3p. These 10 pathways could be involved in the development of OP. Functional Enrichment Analysis of mRNAs in the Regulatory Network GO analysis of the mRNAs in the regulatory network showed that biological progress (BP) terms were signi cantly enriched in positive regulation of oxidoreductase activity, cellular response to lipopolysaccharide, cellular response to molecule of bacterial origin, glycerolipid catabolic process, extrinsic apoptotic signaling pathway, cellular response to biotic stimulus, positive regulation of lipid biosynthetic process, extrinsic apoptotic signaling pathway via death domain receptors, regulation of phospholipid metabolic process, regulation of oxidoreductase activity. The cell component (CC) terms included main axon, mitochondrial intermembrane space, voltage-gated potassium channel complex, organelle envelope lumen, speci c granule membrane, secretory granule membrane, potassium channel complex, acrosomal vesicle, zymogen granule membrane, clathrin-sculpted vesicle. The molecular function (MF) terms were enriched in death receptor binding, tumor necrosis factor receptor superfamily binding, phospholipid binding, neurotrophin receptor binding, phosphatidylinositol 3-kinase regulatory subunit binding, SH3 domain binding, death receptor activity, intracellular cyclic nucleotide activated cation channel activity, glucose binding, cyclic nucleotide-gated ion channel activity (Fig. 5(a)). The most important pathways by KEGG are also shown in (Fig. 5(b)). Including Other types of O-glycan biosynthesis, Aldosterone synthesis and secretion, TNF signaling pathway, Apoptosis, Glycosaminoglycan biosynthesis -keratan sulfate, Vitamin digestion and absorption.

Discussion
Various studies demonstrate that miRNAs play a critical role in the process of pathology relating to OP. Chen, R., et al proved that miRNA-19a-3p alleviates the progression of osteoporosis by targeting HDAC4 to promote the osteogenic differentiation of hMSCs(Human bone marrow mesenchymal stem cells) [17]. Wang, W. W., et al found that miR-218 and miR-618 play an important role in osteoclastogenesis via TLR-4/MyD88/NF-κB signaling pathway [18].
In this study, we ltered hsa-let-7a-5p as a hub node in the regulation of OP. The let-7a induced by TNF-α in osteoporosis, inhibited the expression of the Fas/FasL system via post-transcriptional regulation. Silencing let-7a expression can recover the immunosuppressive capacity of osteoporotic bone marrow mesenchymal stem cells (BMSCs) [19]. By regulating TGFBR1, the let-7a-5p might inhibit the osteogenic differentiation of BMSCs in postmenopausal osteoporosis (PMOP) mice [20]. The miR-92a inhibits the Smad6-mediated runt-related transcription factor 2 degradation, which can facilitate the osteogenic differentiation of BMSCs [21]. Penzkofer, D., et al. proved the deletion of miR-92a will induce skeletal defects [22]. Increasing the expression level of miR-92b will suppress the Cathepsin K expression and Nox4/NF-κB signaling pathway, which can alleviate the progression of osteoporosis [23]. In the one hand the upregulated miR-92b accelerates osteogenesis of BMSCs.
In the other hand, silence of miR-92b-5p inhibits the osteogenesis of BMSCs [24]. According to these literatures, let-7a and miR-92a/92b might play a critical role in the differentiation and proliferation of osteoblasts and osteoclasts in OP.
Hamam, D., et al. demonstrated that overexpression of miR-320c in hMSC enhanced adipocytic differentiation and accelerated formation of mature adipocytes in vitro, which might lead to impact the differentiation of hMSCs into osteoblasts [26].
SLC2A3, LBX1, HCN2, DAB2IP and CIC were selected as the key genes in the miRNA-mRNA network. SLC2A3, which was also named Glucose Transporter Type 3(GLUT3), is the principal transporters of glucose in osteoblasts. This means SLC2A3 might affect the ATP demand for osteoclast differentiation [27]. LBX1 is regulated by the Wnt/B-catenin pathway, which is a canonical pathway for osteoblast differentiation. Therefore, LBX1 might accelerate the differentiation and proliferation of osteoblasts. [28] HCN2 might promote TNF-induced degradation of type 2 collagen and aggrecan, which might cause a reduction in the number of osteoblasts and their differentiation. [29] He, J., et al proved DAB2IP expression in uence cell proliferation, apoptosis and cell cycle distribution. And it inhibited the migration and invasion of normal osteoblast cells. [30] CIC regulates a transcriptional network downstream of the RAS/MAPK signaling cascade. RAS/MAPK is an essential pathway for the activity and differentiation of osteoblasts. [31] To assess the strength of the miRNA-mRNA network analysis, we ltered 10 mRNAs and genes including hsa-let-7a-5p, hsa-miR-92a-3p, hsa-miR-92b-3p, hsa-miR-223-3p, hsa-miR-320c, SLC2A3, LBX1, HCN2, DAB2IP and CIC, which are believed to play critical roles in the development and pathological mechanisms of OP. Among which, hsa-let-7a-5p, hsa-miR-92a-3p and hsa-miR-92b-3p were considered to be the most important functional miRNA, being recognized as the most important signaling pathways in OP.
This study had some limitations. Firstly, due to the lack of experimental veri cation, these conclusions need to be further explored. In addition, the number of samples we obtained from GSE93883, GSE74209 and GSE35959 was small, generating some bias when analysing the miRNA-mRNA interaction network. Authors' Contributions HCB: study concept and design; acquisition of data; analysis and interpretation of data; drafting of the manuscript; critical revision of the manuscript for important intellectual content; statistical analysis. ZDY and XTH: acquisition of data; analysis and interpretation of data; added experiment results. XJ and JC: material support; analysis and interpretation of data; critical revision of the manuscript. YL: revision of the manuscript and study supervision. The author(s) read and approved the nal manuscript.

Availability of data and materials
All data are available in the manuscript and they are showed in gures and tables.

Competing interests
The authors declare that they have no competing interests