Abstract Angiogenesis-osteogenesis coupling is critical for proper functioning and maintaining the health of bones. Any disruption in this coupling, associated with aging and disease, might lead to loss of bone mass. Osteoporosis (OP) is a debilitating bone metabolic disorder that affects the microarchitecture of bones, gradually leading to fracture. Computational analysis revealed that normal angiogenesis is disrupted during the progression of OP, especially postmenopausal osteoporosis (PMOP). The genes associated with OP and PMOP were retrieved from the DisGeNET database. Hub gene analysis and molecular pathway enrichment were performed via the Cytoscape plugins STRING, MCODE, CytoHubba, ClueGO and the web-based tool Enrichr. Twenty-eight (28) hub genes were identified, eight of which were transcription factors (HIF1A, JUN, TP53, ESR1, MYC, PPARG, RUNX2 and SOX9). Analysis of SNPs associated with hub genes via the gnomAD, I-Mutant2.0, MUpro, ConSurf and COACH servers revealed the substitution F201L in IL6 as the most deleterious. The IL6 protein was modeled in the SWISS-MODEL server and the substitution was analyzed via the YASARA FoldX plugin. A positive ΔΔG (1.936) of the F201L mutant indicates that the mutated structure is less stable than the wild-type structure is. Thirteen hub genes, including IL6 and the enriched molecular pathways were found to be profoundly involved in angiogenesis/endothelial function and immune signaling. Mechanical loading of bones through weight-bearing exercises can activate osteoblasts via mechanotransduction leading to increased bone formation. The present study suggests proper mechanical loading of bone as a preventive strategy for PMOP, by which angiogenesis and the immune status of the bone can be maintained. This in silico analysis could be used to understand the molecular etiology of OP and to develop novel therapeutic approaches. Supplementary Information The online version contains supplementary material available at 10.1186/s12863-024-01269-z. Keywords: Hub genes, Osteoporosis, Postmenopausal osteoporosis, Protein-protein interaction network, Reactome, KEGG, Single nucleotide variations Introduction Osteoporosis (OP) is a serious disorder of bone metabolism. The symptoms of this condition are fragility fractures and low bone mineral density (BMD) [[36]1]. It is diagnosed by assessing the bone mineral density, for which dual-energy X-ray absorptiometry (DXA) is commonly used. A T-score of -2.5 or less implies OP [[37]2]. According to the International Osteoporosis Foundation, every 3 s, an osteoporotic fracture occur in the human population, resulting in 8.9 million fractures annually. In Europe, 22.1% of women and 6.6% of men aged greater than 50 years were affected by OP in 2019 [[38]3]. The health-related quality of life of postmenopausal women with OP is affected to variable extents [[39]4]. Postmenopausal women show a consequential decline in BMD from the age of 45–50, whereas in premenopausal women, no such decline is observed until the age of 55 [[40]5]. A person can also be affected by secondary OP, i.e. OP due to an underlying disease or a different metabolic condition. For example, the prevalence of OP in patients with chronic obstructive pulmonary disease (COPD) is 52.8% in men, as reported by a study carried out in Qazvin, whereas the prevalence of COPD in the community is 5–13% [[41]6]. Normal bone mass is maintained in the body via bone remodeling via signals such as mechanical loading, microfractures, stress, and hormones [[42]7]. Increased production of osteoclasts and decreased production of osteoblasts result in enhanced bone loss and OP. In patients with PMOP, a decrease in estrogen levels profoundly influences the induction of bone loss. Estrogen has the capacity to regulate the secretion of osteoclastogenic cytokines in the body. Loss of estrogen increases the production of proinflammatory cytokines viz. IL-1 in the bone remodeling circuitry that can induce bone resorption and result in bone loss [[43]8]. Understanding the molecular signals that regulate osteoclasts and osteoblasts is essential in the development of drugs and treatment strategies [[44]9]. Recently, increased focus has been given to the crosstalk between molecular signaling pathways that are involved in the maintenance of BMD. RANKL/RANK/OPG pathway plays a crucial role in osteoclast differentiation [[45]10] and the Wnt signaling pathway is integral for bone formation [[46]11]. Angiogenesis-osteogenesis coupling is essential for the development of bone. A recent study in juvenile mice revealed the importance of endothelial SMAD1/5 signaling in maintaining the integrity of metaphyseal vessels. The depletion of SMAD1/5 resulted in excessive sprouting of metaphyseal vessels with disrupted morphology and altered bone formation [[47]12]. Type H vessels in the bone marrow are a specialized capillary subtype that is involved in the coupling of angiogenesis and osteogenesis [[48]13]. Physical activity plays a crucial role in maintaining bone mineral density. Mechanical loading of bones through weight-bearing exercises can activate osteoblasts via mechanotransduction leading to increased bone formation. Studies have shown that individuals with regular physical activity have higher BMD and a lower risk of OP [[49]14]. Mechanical forces produced during physical activity, such as shear stress induced by blood flow, can lead to remodeling of the extracellular matrix in endothelial cells and influence angiogenesis [[50]15]. The impairment of angiogenesis due to changes in the ECM caused by cellular senescence during aging [[51]16] can negatively affect bone quality and lead to senile osteoporosis. Various treatment options are used for OP, such as calcium and vitamin D supplements, antiresorptive agents, anabolic drugs that promote bone formation and drugs with bidirectional regulation i.e., those that promote bone formation and inhibit bone resorption [[52]9]. A synthetic peptide analog of the parathyroid hormone related protein abaloparatide was recently approved by the US Food Drug Administration (FDA) for use in high risk PMOP patients [[53]17]. Teriparatide, a synthetic hormone, human recombinant PTH (1–34) is used to increase bone mineral density in PMOP patients [[54]18]. Romosozumab, an approved humanized monoclonal antibody, is used for the treatment of PMOP. This antibody inhibits sclerostin, stimulates bone formation and inhibits bone resorption. However, this antibody should not be used by patients with cardiovascular problems [[55]19]. Combination and sequential therapy are now considered as the most effective methods for the restoration of bone mass. Most commonly, an osteoanabolic agent is followed by an antiresorptive agent and at each level, a combination of different drugs is used [[56]20]. Through this strategy, OP patients benefit from the bone mass-increasing effects of anabolic drugs and the bone mass-preserving effects of antiresorptive drugs. Recent studies suggest that the gut microbiome plays a major role in bone metabolism. The gut microbiota influences bone physiology by regulating the immune system. Treatment with probiotics has been shown to increase bone mass in mice [[57]21]. The intestinal flora and their metabolites can be used to monitor patients with PMOP [[58]22]. Osteoporosis treatments can have various adverse side effects according to the therapy used. For example, the risks of osteonecrosis of the jawbone (ONJ) and atypical femoral fractures are greater in patients receiving high- and long-term use of bisphosphonates [[59]9]. The monoclonal antibody denosumab can lead to hypocalcemia and ONJ [[60]23]. Selective estrogen receptor modulators (SERMs) such as raloxifene, increase the risk of venous thromboembolism (VTE) and hormone replacement therapy is associated with the risk of cardiovascular events [[61]9]. Identifying biomarkers that can predict a patient’s response to OP treatments will help minimize adverse effects. The focus is now on improving the overall quality of the bone rather than just increasing the bone mass. Stem cell therapy, miRNA-based therapy and bone-specific targeting technology are gaining interest as the latest treatment options for OP, and studies of molecular pathways and miRNAs have become essential [[62]9]. The molecular etiology of OP is complex and has not been fully established [[63]24]. Several metabolic, immune, cell cycle and other pathways are perturbed during the progression of OP [[64]25, [65]26]. In recent years, pathway enrichment analysis and databases providing experimental data have become popular and reliable tools for the initial characterization of genes, molecules or pathways that are to be targeted for therapy [[66]27]. STRING, GeneMania, PathwayCommons, David, ClueGO and Panther are tools used to identify protein-protein interactions (PPIs), hub genes and pathways related to a set of genes. Novel treatments for osteoporotic patients could be deduced by utilizing newly identified genes and pathways [[67]28]. Hub genes play important roles in the gene regulatory network and they organize the behaviour of the network. These genes interact with many genes in the network and control key biological processes [[68]29]. Mutations in these genes can affect the function of many other genes, leading to different phenotypic consequences. Missense variations are single nucleotide variations that result in the substitution of one amino acid by another. This can affect protein structure and function. A change in the amino acid sequence can be benign or deleterious. Deleterious missense variations can result in loss of function, gain of harmful function or an altered function. If these missense variations occur in hub genes and if they are deleterious, they can affect the function of a range of other genes and molecular pathways to which the hub gene is associated with [[69]30]. Identifying missense mutations in hub genes would help to predict the pathogenic potential of these mutations. The primary goal of this study was to understand the disease mechanism of OP and PMOP. The focus of this work was to identify hub genes and molecular pathways, possibly dysregulated during the development and progression of this metabolic disorder. Hub genes are the highly interconnected genes in a gene regulatory network and play crucial roles in various biological processes [[70]31]. They act as key regulators that interact with and regulate multiple other genes and pathways. In coexpression and regulatory gene networks, they are associated with cellular functions. These genes may be associated with disease mechanisms when dysregulated and their absence can disrupt the pathway function. Hub genes are identified on the basis of their high connectivity with specific modules of a network [[71]30, [72]32]. Identification of the hub genes involved in osteoporosis would help to identify crucial genes which are involved in maintaining bone mineral density. These hub genes can act as biomarkers of this metabolic condition. Early diagnosis, therapeutic interventions and monitoring of disease progression could be performed with the help of these biomarkers and pathways. This would be a great contribution to the design of personalized medicine, novel diagnostic measures and preventive measures. Women, especially those nearing menopause or having reached menopause and the geriatric population, highly benefit from the design of proper preventive strategies. In this study, genes associated with OP and PMOP were retrieved from the DisGeNET database and was subjected to computational analysis (Fig. [73]1). Fig. 1. [74]Fig. 1 [75]Open in a new tab Flow diagram showing the work strategy Methods Protein–protein interaction network and gene clusters The complete set of genes associated with OP and PMOP were downloaded from the DisGeNET database (supplementary data). DisGeNET is a knowledge database, a large repository of human gene–disease associations and variant–disease associations. To download the complete set of genes, one should register and login to the database. The Cytoscape STRING app was used to construct a protein–protein interaction (PPI) network. The full lists of genes associated with OP and PMOP were pasted into the STRING identifier separately. The individual networks created for OP and PMOP were used for further analysis. Clusters of highly connected areas in the PPI network which signifies protein complexes, were identified via molecular complex detection (MCODE), a Cytoscape app. Cluster detection was performed on the previously identified PPI networks via STRING (supplementary data). The scores of different cluster modules and the genes in each module were noted. Hub gene identification Hub genes were identified via the Cytoscape plugin CytoHubba, in the PPI network of OP and PMOP obtained from Cytoscape STRING. The top 20 Hub nodes ranked by the 12 topological analysis methods were established (supplementary data). Genes present in at least half of the centrality analyses were identified via commonality analysis. The presence of hub genes was then confirmed in the cluster modules identified via MCODE analysis. Single nucleotide variations (SNVs) SNVs associated with the hub genes were identified via the DisGeNET, PubMed and Google Scholar databases. Missense variations were shortlisted by checking the rsIDs in NCBI-SNP. The rsIDs are unique labels used to identify an SNP. Loss-of-function effect of the variant was predicted via the tool gnomAD. The variants were given as inputs one by one in the gnomAD v4.1.0 browser. Only those missense SNPs that passed all variant filters in gnomAD were subjected to structure stability prediction via I-Mutant2.0 and MUpro. I-Mutant2.0 ([76]https://folding.biofold.org/i-mutant/) predicts protein stability changes when there is a single nucleotide mutation. The protein sequence, position and the new residue were given as inputs. The prediction is displayed as the free energy change value (DDG). If the value is less than 0, the substitution decreases stability and vice versa. The same method is also followed in MUpro too ([77]https://mupro.proteomics.ics.uci.edu/). MUpro also predicts changes in protein stability for single-nucleotide variants (SNVs) via amino acid sequences. The sequence without headers, the position of mutation, the original amino acid and the substitute amino acid were given as inputs. The result is obtained as a DDG value and is predicted as decreased or increased stability. The ConSurf webserver ([78]https://consurf.tau.ac.il/consurf_index.php) was used to identify whether the variation was in the functional region of the protein products of the hub genes. The FASTA sequence was used as input for the prediction. MusiteDeep ([79]https://www.musite.net/) was used for the prediction of posttranslational modifications (PTMs) at the site of variation. The FASTA sequence was submitted to this server and a prediction model was selected to identify the PTMs. Protein ligand binding site prediction was performed via the COACH server. It uses predictions from TM-SITE, S-SITE, COFACTOR, FINDSITE and ConCavity to generate final ligand binding site predictions. The native protein structure of the IL-6 gene with an SNV (F201L), identified as most deleterious was modeled via the SWISS-MODEL server and ΔΔG was calculated via the YASARA FoldX plugin. Molecular docking studies were carried out in AutoDock 1.5.7 to assess the impact of the F201L mutation. The mutant IL-6 protein was generated via SWISS-MODEL, which employs the wild-type IL-6 structure (UniProt ID: [80]P05231) as a template. Molecular pathways Genes retrieved from DisGeNET were used as inputs to discern molecular pathways and were visualized via ClueGO and Enrichr. Overlap analysis of the KEGG and Reactome pathways (p value ≤ 0.05) was performed via multiple list comparators and Venn diagrams of molbiotools to identify pathways unique to PMOP. Literature mining and the GenomeNet and Reactome databases were utilized to connect pathways to various functions. GO molecular functions of hub genes The GO molecular functions of the hub genes were retrieved from Enrichr. For this purpose, the genes were separately entered into the Enrichr tool as input data and the molecular functions of each gene were used to tabulate the key functions of the hub genes. Results Cluster modules from the PPI network From DisGeNET, 1098 genes associated with osteoporosis were retrieved for OP and 171 genes were retrieved for PMOP. A total of 33 MCODE clusters were identified from the STRING protein–protein interaction network for OP and these were visualized via Cytoscape. The highest score cluster module had 97 nodes and 3567 edges with a score of 74.312 and the second cluster with a score of 19.286 had 71 nodes and 675 edges. All 20 hub genes associated with osteoporosis were located in the first two MCODE clusters. Three gene clusters were identified in the PPI network created for PMOP. Cluster 1 had 46 nodes and 699 edges with the highest score being 31.067. The second cluster had 16 nodes and 39 edges with a score of 5.200. The third MCODE cluster had 6 nodes and 9 edges with a score of 3.600. Hub genes of OP and PMOP Overlap analysis was performed among the top 20 Hub nodes ranked by 12 topological analysis methods in CytoHubba and genes present in at least 6 topological analyses were selected. The twenty hub genes identified for OP are ACTB, AKT1, ALB, GAPDH, HIF1A, IL1B, IL6, INS, JUN, MAPK3, STAT3, TGFB1, TNF, TP53, CTNNB1, EGFR, ESR1, MYC, PPARG and SRC. The sixteen hub genes identified for PMOP are ESR1, CTNNB1, GAPDH, IGF1, IL1B, IL6, MAPK3, SPP1, TNF, VEGFA, RUNX2, TLR4, BGLAP, BMP2, HIF1A and SOX9. The identified hub genes were confirmed by their presence in the prominent MCODE clusters as shown in Fig. [81]2. When pooled, the total number of hub genes identified for OP and PMOP was 28. Fig. 2. [82]Fig. 2 [83]Open in a new tab Hub genes (Hub nodes) visualized in the MCODE clusters from OP- (A, B) and PMOP- (C, D) associated genes (highlighted yellow). A Hub nodes in the MCODE cluster module with a score of 74.312 B Hub nodes in the MCODE cluster module with a score of 19.286 C Hub nodes in the MCODE cluster module with a score of 31.067 D. Hub nodes in the MCODE cluster module with a score of 5.200 Commonality analysis of hub genes Venn diagrams were drawn from commonality analysis to identify hub genes that are unique to PMOP and OP and are shown in Table [84]1. The hub genes common to OP and PMOP as shown by commonality analysis are CTNNB1, ESR1, GAPDH, HIF1A, IL1B, IL6, MAPK3 and TNF (Fig. [85]3). Those unique to PMOP are IGF1, SPP1, VEGFA, RUNX2, TLR4, BGLAP, BMP2 and SOX9. Table 1. Results of commonality analysis of hub genes identified from OP-associated and PMOP-associated genes. *bold – TFs PMOP hub genes ∩ OP hub genes (common) PMOP hub genes minus OP hub genes (only in PMOP) OP hub genes minus PMOP hub genes (only in OP) CTNNB1 IGF1 ACTB ESR1 SPP1 AKT1 GAPDH VEGFA ALB HIF1A RUNX2 INS IL1B TLR4 JUN IL6 BGLAP STAT3 MAPK3 BMP2 TGFB1 TNF SOX9 TP53 EGFR MYC PPARG SRC [86]Open in a new tab Fig. 3. [87]Fig. 3 [88]Open in a new tab MCODE clusters and Hub genes A Cluster genes and Hub genes for OP-associated genes B Cluster genes and Hub genes for PMOP-associated genes Key functions of the hub genes The functions of the hub genes identified from the genes associated with OP and PMOP were obtained from GeneCards, the Human Gene Database, UniProt and the NCBI-gene database. The GO molecular functions of the hub genes that were retrieved from Enrichr are given in the supplementary data. These results were used to organize the key functions of the hub genes, as shown in Tables [89]2 and [90]3. Notably, eight of these genes were TFs and many others are closely associated with transcription. Furthermore, these genes play crucial roles in angiogenesis, endothelial function and immune signaling. Table 2. Key functions of the protein products of 16 hub genes of PMOP according to GeneCards, NCBI-gene, UniProt, Enrichr GO molecular function and literature mining Gene symbol Gene name Functions of protein product as per NCBI-gene and UniProt ESR1 Estrogen receptor 1 (alpha) Transcription factor, nuclear hormone receptor, modulate angiogenesis in adipose tissue CTNNB1 Catenin beta 1 Intracellular signal transducer, coactivator for transcription factors of the TCF/LEF family GAPDH Glyceraldehyde-3-phosphate dehydrogenase Enzyme, Component of OCA-S transcriptional coactivator complex IGF1 Insulin like growth factor 1 Hormone activity, growth factor activity IL1B Interleukin 1 Beta Pro-inflammatory cytokine, role in angiogenesis IL6 Interleukin 6 Cytokine, role in angiogenesis MAPK3 Mitogen-activated protein kinase 3 Kinase, phosphorylation of transcription factors SPP1 Secreted phosphoprotein 1 Osteopontin, cytokine, bone protein, extracellular matrix binding viz. hydroxyapatite binding TNF Tumor necrosis factor An adipokine and a cytokine, role in angiogenesis VEGFA Vascular Endothelial Growth Factor A Heparin-binding protein, cytokine, role in angiogenesis RUNX2 RUNX Family Transcription Factor 2 Transcription factor, osteoblastic differentiation, role in angiogenesis and bone vascularization TLR4 Toll Like Receptor 4 Single transmembrane cell-surface receptors, mediate cytokine secretion and inflammatory response, role in angiogenesis BGLAP Bone Gamma-Carboxyglutamate Protein Osteocalcin, bone protein, calcium ion binding, hydroxyapatite binding, role in angiogenesis BMP2 Bone Morphogenetic Protein 2 Secreted ligand of the TGF-beta receptors, cytokine, role in angiogenesis HIF1A Hypoxia Inducible Factor 1 Subunit Alpha Transcription factor, role in angiogenesis SOX9 SRY-Box Transcription Factor 9 Transcription factor [91]Open in a new tab Table 3. Key functions of 20 hub genes associated with OP according to GeneCards, NCBI-Gene, UniProt, Enrichr GO molecular function and literature mining Gene symbol Gene name Functions of protein product as per NCBI-gene and UniProt ACTB Actin Beta Actin protein, intercellular signaling AKT1 AKT Serine/Threonine Kinase 1 Kinase, role in angiogenesis ALB Albumin Regulation of the colloidal osmotic pressure of blood, major calcium transporter in plasma GAPDH Glyceraldehyde-3-phosphate dehydrogenase Enzyme, Component of OCA-S transcriptional coactivator complex HIF1A Hypoxia Inducible Factor 1 Subunit Alpha Transcription factor, role in angiogenesis IL1B Interleukin 1 Beta Pro-inflammatory cytokine, role in angiogenesis IL6 Interleukin 6 Cytokine, role in angiogenesis INS Insulin Peptide hormone, stimulates glucose uptake JUN Jun proto-oncogene Transcription factor MAPK3 mitogen-activated protein kinase 3 Kinase, phosphorylation of transcription factors STAT3 Signal Transducer and Activator of Transcription 3 Signal transducer and transcription activator, acts as a regulator of inflammatory response, role in angiogenesis TGFB1 Transforming growth factor beta 1 Protein ligand, role in bone remodeling, osteoblast differentiation TNF Tumor necrosis factor An adipokine and a cytokine, role in angiogenesis TP53 Tumor Protein P53 Transcription factor CTNNB1 Catenin beta 1 Intracellular signal transducer, coactivator for transcription factors of the TCF/LEF family EGFR Epidermal Growth Factor Receptor Protein kinase ESR1 Estrogen receptor 1 (alpha) Transcription factor, nuclear hormone receptor, modulate angiogenesis in adipose tissue MYC MYC proto-oncogene Transcription factor, promote VEGFA production, role in angiogenesis PPARG peroxisome proliferator activated receptor gamma Transcription factor, regulator of adipocyte differentiation SRC SRC proto-oncogene Tyrosine-protein kinase, role in osteoclastic bone resorption, mediates IL6 signaling [92]Open in a new tab Predictions of single-nucleotide variations Missense variations were detected in 7 hub genes: IL6, IGF1, TLR4, BGLAP, BMP2, TGFB1 and TP53. The list of SNPs collected and the results of the analysis with various tools are shown in Table [93]4. The position of the SNP rs2069849 was predicted as a conserved region by ConSurf as shown in Fig. [94]4. In the COACH server which predicts ligand binding sites, predictions from TM-SITE, S-SITE, COFACTOR, FINDSITE and ConCavity results are available. An image from the ConCavity results showing position 201 in IL6 as the ligand binding site is shown in Fig. [95]5. The F201L variation in IL6 is predicted to be the SNP with the greatest potential, which can lead to the development of OP. The substitution has a free energy change of 1.936 in YASARA FoldX. This shows that the mutated structure is less stable than the wild-type structure is. The results of the molecular docking studies for the IL-6 protein with wild-type and mutated structures carried out in AutoDock 1.5.7 are shown in Table [96]5; Fig. [97]6. Table 4. Results from gnomAD, I-Mutant2.0, Mupro, ConSurf and COACH gnomAD  Hub gene SNP Variation Result Effect   IL6 rs2069849 C-T pass *PLoF   IGF1 rs35767 A-C not pass -   TLR4 rs1057317 C-A not pass -   BGLAP rs1800247 T-C pass with warning -   BMP2 rs2273073 T-G pass *PLoF   TGFB1 rs1800470 G-A pass *PLoF   TP53 rs1042522 G-C pass *PLoF *PLoF - Predicted loss of function  I-Mutant2.0   Hub gene SNP Variation **DDG Effect   IL6 rs2069849 F-L -1.35 Decrease stability   BMP2 rs2273073 S-A 0.01 Increase stability   TGFB1 rs1800470 P-L -0.85 Decrease stability   TP53 rs1042522 P-R 0.38 Increase stability **DDG - delta delta G, DDG < 0: Decrease Stability, DDG > 0: Increase Stability  MUpro   Hub gene SNP Variation **DDG Effect   IL6 rs2069849 F-L -0.58605643 Decrease stability   BMP2 rs2273073 S-A -0.2650236 Decrease stability   TGFB1 rs1800470 P-L -0.050797501 Decrease stability   TP53 rs1042522 P-R 0.16475502 Increase stability **DDG - delta delta G, DDG < 0: Decrease Stability, DDG > 0: Increase Stability  ConSurf   Hub gene SNP Position Region   IL6 rs2069849 201 conserved region   BMP2 rs2273073 37 variable region   TGFB1 rs1800470 10 variable region  COACH   Hub gene SNP Position Prediction Score   IL6 rs2069849 201 consensus binding residue 0.02   BMP2 rs2273073 37 not a binding residue -   TGFB1 rs1800470 10 not a binding residue -  YASARA FoldX   Hub gene SNP Variation Prediction ΔΔG   IL6 rs2069849 F201L Unstable 1.936 ΔΔG < 0: Increase Stability, ΔΔG > 0: Decrease Stability [98]Open in a new tab Fig. 4. [99]Fig. 4 [100]Open in a new tab Results from the evolutionary conservation prediction tool ConSurf server showing the conservation score of the amino acids of the IL-6 protein Fig. 5. [101]Fig. 5 [102]Open in a new tab A ConCavity results of ligand binding site prediction for position 201 of IL-6 from the COACH server. F201 is shown as a yellow circle. B 3D image of the protein structure modeled by the SWISS-MODEL server for the native IL-6 protein Table 5. Results from molecular docking studies of IL-6 protein Protein Ligand Binding energy Ligand efficiency Inhibition constant Intermolecular energy Vdw_hb_ dissolvation energy Electrostatic energy No. of hydrogen bond Wild IL-6 GP 130 -5.05 -0.36 200.28 -5.34 -4.3 -1.04 3 Mutant IL-6 (F→L) -4.6 -0.33 427.57 -4.89 -3.59 -1.31 1 [103]Open in a new tab Fig. 6. [104]Fig. 6 [105]Open in a new tab 3D and 2D interaction of (A & B) mutant IL6 (F→L) with GP130 and (C & D) wild-type IL6 with GP130 Pathway enrichment analysis Pathway enrichment analysis in ClueGO for OP genes revealed 199 significant KEGG pathways and 209 Reactome pathways with p values ≤ 0.05. Signaling by cytokines, various interleukins, FGFRs, estrogen, growth hormone, PI3K/AKT, BMP, insulin, leptin, NTRKs, receptor tyrosine kinases, TGF-beta, IGF1R and WNT were found under the significant signaling pathways. Transcription pathways mediated by FOXO, SMADs, RUNX2, RUNX3 and HIF-1 were found to be significant (Fig. [106]7). In the enrichment analysis of genes associated with PMOP, 180 KEGG pathways and 44 Reactome pathways with p values ≤ 0.05 were obtained. Cytokine signaling, signaling by DAP12, extranuclear estrogen, interleukins, TCF, WNT, PI3K/AKT, platelet, activin, erythropoietin, VEGF and transcriptional pathways involving RUNX2 and FOXO were found to be enriched and all pathways are shown in the supplementary data. Enrichment analysis was also performed in Enrichr for both sets of genes (Fig. [107]8). Fig. 7. [108]Fig. 7 [109]Open in a new tab KEGG pathway enrichment in ClueGO and Enrichr for OP associated genes. A. Visualization of pathway enrichment in ClueGO B. Percentage of terms per group C. Pathway enrichment from Enrichr Fig. 8. [110]Fig. 8 [111]Open in a new tab KEGG pathway enrichment in ClueGO and Enrichr for PMOP-associated genes. A Visualization of pathway enrichment in ClueGO B Percentage of terms per group C. Pathway enrichment from Enrichr Commonality analysis of molecular pathways Dataset intersections of enriched KEGG pathways from gene set analysis of PMOP and OP in ClueGO were detected in a multiple list comparator. Shared and unique pathways of the gene sets associated with PMOP and OP are listed in the supplementary data. The Ras signaling pathway, the VEGF signaling pathway, aldosterone-regulated sodium reabsorption, natural killer cell-mediated cytotoxicity, Alzheimer’s disease, the cAMP signaling pathway, the sphingolipid signaling pathway, axon guidance, the hematopoietic cell lineage, the B cell receptor signaling pathway, the Fc epsilon RI signaling pathway, Fc gamma R-mediated phagocytosis, inflammatory mediator regulation of TRP channels and GnRH secretion are the KEGG pathways unique to PMOP. The Reactome pathways unique to PMOP include the regulation of Insulin-like Growth Factor (IGF) transport and uptake by insulin-like growth factor binding proteins (IGFBPs), ECM proteoglycans, GP1b-IX-V activation signaling, signaling by activin, DAP12 interactions, DAP12 signaling, FCERI mediated MAPK activation, RHO GTPases Activate NADPH Oxidases, GPVI-mediated activation cascade, signaling by erythropoietin and erythropoietin activates phosphoinositide-3-kinase (PI3K). The unique KEGG pathways associated with PMOP shows that it belongs to classes -Immune system, signal transduction, endocrine system, development and regeneration, excretory system, sensory system and neurodegenerative disease. Among the Reactome pathways unique to PMOP, the parent classes included metabolism of proteins, extracellular matrix organization, platelet activation, signaling by TGFB family, innate immune system, signal transduction, signaling by Rho GTPases and Erythropoietin (Fig. [112]9). Most Reactome pathways are immune signaling pathways associated with angiogenesis/endothelial function. A comparative analysis of the top KEGG pathways and unique PMOP pathways is given in Table [113]6. Fig. 9. [114]Fig. 9 [115]Open in a new tab KEGG and Reactome pathways unique to PMOP and their parent classes. A. KEGG class and the corresponding KEGG pathways. B. Reactome event hierarchy and the corresponding Reactome pathways Table 6. Comparative study with previously published data for the top KEGG pathways and pathways unique to PMOP Pathways Role in bone physiology/osteoporosis Role in immune signaling or angiogenesis/endothelial function MAPK signaling pathway Inhibition of MAPK pathway reduces inflammatory bone loss [[116]33] at the same time MAPK signaling is essential for bone formation [[117]34]. MAPK signaling is a mediator through which osteocytes regulates angiogenesis in bone [[118]35]. JAK-STAT signaling pathway Has important role in osteogenesis and in sensing mechanical stress, can regulate mechanically induced bone formation [[119]36]. Activation of JAK-STAT signaling pathway can promote angiogenesis and leptin can induce this activation [[120]37]. Leptin is essential for maintaining bone mass and quality by stimulating osteoblast differentiation [[121]38]. Leukocyte transendothelial migration Found to be higher in OP patients. Enhances osteoclastogenesis [[122]39]. Crucial for inflammation [[123]40]. Parathyroid hormone synthesis, secretion and action Elevated level of Parathyroid hormone stimulates bone resorption and increase bone turn over resulting in loss of BMD [[124]41] whereas intermittent elevation of PTH have anabolic effects in bone [[125]42]. PTH augments angiopoietin-1 release and promote angiogenesis [[126]43, [127]44]. Rheumatoid arthritis Pro-inflammatory cytokines secreted, leads to bone resorption and bone loss [[128]45]. Inflammation and angiogenesis are the key features [[129]46]. Type I diabetes mellitus (T1DM) A risk factor for OP since diabetes contribute to lower peak bone mass [[130]47]. Chronic inflammation in T1DM may promote pathological angiogenesis or inhibit normal angiogenesis. Chronic hyperglycemia can cause endothelial dysfunction [[131]48]. NF-kappa B signaling pathway RANK ligand-induced osteoclastogenesis is mediated by NF-kappa B [[132]49]. Involved in the suppression of angiogenesis [[133]50] PI3K-Akt signaling pathway Involved in the regulation of OP. Promote osteoblast proliferation and differentiation [[134]51]. Expression of angiogenic factors is modulated. Regulate angiogenesis [[135]52]. Wnt signaling pathway Promote osteoblast differentiation. Blocks osteoblasts and osteocyte apoptosis. Regulates osteoclastogenesis [[136]11]. Through it’s interaction with other angiogenic signaling such as Notch Signaling and VEGF signaling, promotes or inhibits angiogenesis. Macrophage-derived Wnt signaling is involved in physiologic as well as pathologic angiogenesis [[137]53]. Apelin signaling pathway Apelin are expressed in osteoblasts and it acts as a mitogenic agent for osteoblasts [[138]54]. Shown to regulate osteogenic differentiation of human mesenchymal stem cells [[139]55]. Apelin is suggested as a therapeutic target for endothelial dysfunction based on the studies in human umbilical cord endothelial cells (EC) and it is also having role in promoting proliferation of EC and angiogenesis [[140]56]. Cytokine-cytokine receptor interaction Production of pro-inflammatory cytokines IL-1β, IL-6, TNF-α were found to be higher in osteoporotic post-menopausal women [[141]57]. Cytokine RANKL stimulates osteoclastogenesis and regulate osteoclastic bone resorption [[142]10]. Cytokines in dysregulated states can lead to osteoporosis [[143]58]. Production of pro-inflammatory cytokines can induce immune signaling that culminate in inflammation [[144]46]. TGF-beta signaling pathway Have important role in the maintenance of bone mass by interacting with other cytokines in the bone tissue. Elevated levels induce inflammatory osteoclastogenesis [[145]59] TGF-beta signaling upregulate proangiogenic factors in cancers with highly vascularized networks [[146]60]. Hippo signaling pathway Involved in turning on osteoclastogenesis, facilitating osteoporosis [[147]61, [148]62]. Involved in modulating angiogenesis [[149]63]. Inflammatory bowel disease Formation of pro-inflammatory cytokines IL-1β, IL-6, TNF-α in inflammatory bowel disease affects bone resorption and contributes to the development of OP [[150]64]. Immune signaling that leads to inflammatory condition is seen. There is pathological angiogenesis where Hippo pathway has a significant role [[151]65]. Rap1 signaling pathway Regulates osteoblast differentiation [[152]66]. Regulate integrins that facilitate cell-ECM adhesion which is integral for endothelial cell movement towards angiogenic signals [[153]67]. Thyroid hormone signaling pathway Have essential role in bone remodeling and energy metabolism in bone [[154]7]. Known to regulate angiogenesis in tumour. Thyroid hormones bind to integrin αvβ3 on the endothelial cell surface and promote endothelial cell proliferation, migration and capillary formation via activating MAPK/ERK signaling pathway. Modulate the expression of pro-angiogenic factors [[155]68, [156]69] HIF-1 signaling pathway Modulates osteoclast formation and bone resorption [[157]70]. Involved in osteogenesis-angiogenesis coupling [[158]71]. Regulation of Insulin-like Growth Factor (IGF) transport and uptake by Insulin-like Growth Factor Binding Proteins (IGFBPs) Insulin-like growth factor type 1 (IGF-1), functions in balancing and coupling bone resorption and bone formation [[159]72]. IGF binding proteins 3 and 4 regulate sprouting angiogenesis. Insulin-like growth factor 2 is found to be essential for sprouting angiogenesis [[160]73]. ECM proteoglycans Modulate osteoblastogenesis from bone marrow stromal cells [[161]74]. Has role in biomineralization and bone formation [[162]75]. Extracellular matrix remodeling in endothelial cells due to physical forces can lead to the induction of angiogenesis [[163]15]. GP1b-IX-V activation signalling GP1b-IX-V, the receptor on platelets is implicated in inflammation [[164]76]. The pro-inflammatory cytokines released can activate osteoclasts leading to bone loss. GPIb-IX-V is involved in inflammatory complications, in addition to its primary role in haemostasis and thrombosis [[165]77] Signaling by Activin Has control over matrix mineralization [[166]78]. Under inflammatory environment activin A reduce vasculogenesis [[167]79] DAP12 signaling DAP12, the signaling adaptor protein is involved in regulating the formation of multinucleated osteoclasts [[168]80]. DAP-12 is involved in osteoclast formation [[169]81]. Known to mediate microglial inflammatory response [[170]82]. FcεRI mediated MAPK activation FcεRI is the IgE receptor expressed on mast cells [[171]83]. Mast cells are seen to undergo clonal proliferation, called systemic mastocytosis and is seen associated with certain cases of osteoporosis [[172]84]. According to database GeneCards FcεRI mediates IgE effector functions in myeloid cells resulting in inflammatory response. Mast cell activation can lead to neovascularization and angiogenesis, where the activation is by the IgE receptor FcεRI [[173]85]. Signaling by Erythropoietin Erythropoietin induces bone formation. Bone loss is also stimulated during erythropoietin stimulated erythropoiesis [[174]86] Mediate physiological as well as pathological angiogenesis [[175]87]. Ras signaling pathway Associated with bone formation [[176]88]. Ras/MAPK signaling pathway is important for bone homeostasis [[177]89]. IL-6 signaling activates Ras pathway [[178]90]. Ras signaling regulates vascular endothelial growth factor and can stimulate endothelial cells [[179]91]. VEGF signaling pathway Promote angiogenesis and bone repair [[180]92]. Involved in the coupling of angiogenesis and osteogenesis [[181]93]. Plays pivotal role in angiogenesis [[182]94]. Aldosterone-regulated sodium reabsorption Have significant impact on calcium metabolism. Aldosterone can act through the mineralocorticoid receptor on osteoclast cells and result in bone resorption to release calcium. Excessive aldosterone can lead to bone loss [[183]95, [184]96]. Inhibit angiogenesis by its action on VEGFR [[185]97]. [186]Open in a new tab Discussion Understanding the etiology of OP is crucial, since the condition affects the most dynamic tissue in the human body. The hub genes associated with OP identified are ACTB, AKT1, ALB, GAPDH, HIF1A, IL1B, IL6, INS, JUN, MAPK3, STAT3, TGFB1, TNF, TP53, CTNNB1, EGFR, ESR1, MYC, PPARG and SRC. Several experimental studies have correlated these genes with OP and many of them were found to ameliorate OP when pathways guided by them are altered. The flavonoid baicalein [[187]98] and the drug tanshinone [[188]99] suppress OP via the AKT signaling pathway. Interleukin-6 (IL-6) is an osteoclastogenic cytokine that mediates the resorption of bone, at altered levels [[189]100, [190]101]. Signal transducer and activator of transcription 3 (STAT3) plays a role in the regulation of bone remodeling via various signaling cytokines [[191]102], and when STAT3 signaling is altered, osteoclast-mediated bone resorption is enhanced in the bone structure. It also plays an important role in osteogenesis–angiogenesis coupling by influencing VEGF production [[192]103, [193]104]. VEGFA stimulates angiogenesis in bone tissue and enhances the formation, differentiation and function of osteoblasts, resulting in an increase in BMD [[194]105]. RUNX2 encodes a nuclear protein that functions as a TF and confers osteogenic effects [[195]106] and it is involved in the vascularization of bone tissue via the induction of angiogenesis along with HIF-1 alpha [[196]107]. The 16 hub genes associated with postmenopausal osteoporosis were ESR1, CTNNB1, GAPDH, IGF1, IL1B, IL6, MAPK3, SPP1, TNF, VEGFA, RUNX2, TLR4, BGLAP, BMP2, HIF1A and SOX9. Among them CTNNB1, ESR1, IGF1, IL1B, IL6, MAPK3, RUNX2, TNF and VEGFA are common in OP and PMOP. GAPDH, SPP1, TLR4, BGLAP, BMP2, HIF1A and SOX9 are unique to the PMOP set. Studies involving stimulation and inhibition of the TLR4 pathway confirm its role in angiogenesis [[197]108–[198]111]. Osteoblasts secrete the bone gamma-carboxyglutamate protein encoded by BGLAP, also called osteocalcin (OC), which has an important role in bone remodeling [[199]112] and angiogenesis [[200]113]. BMP2 has a role in the epigenetic control of OP [[201]114] and its deviant signaling is detected in OP patients [[202]115]. It plays an important role in osteogenesis–angiogenesis coupling [[203]116, [204]117]. Thirteen of the total 28 hub genes identified from the OP and PMOP gene sets were found to play important roles in angiogenesis (Tables [205]2 and [206]3). SNP analysis of single nucleotide variations associated with hub genes viz. predictions of loss of function, structural stability, functional region and ligand binding revealed that the F201L variation in IL-6 as the most deleterious variant. This mutated protein is less stable than the native protein according to the ΔΔG value obtained from YASARA FoldX. There are chances for the mutated structure to be eliminated from the cellular system. This would hamper the ability of IL-6 to maintain its function in a particular tissue. The IL-6 receptor consists of a protein complex consisting of IL-6R and the signal transducer gp130. IL-6 after binding with membrane bound IL-6R subsequently bonds with the signaling receptor protein gp130 [[207]90]. Comparative analysis of the binding energies revealed a significant decrease in the binding affinity of mutant IL-6 for gp130 (-4.6 kcal/mol) compared with that of the wild-type protein (-5.05 kcal/mol). This observation, coupled with a concomitant increase in the inhibition constant (Ki) from 200.28 nM to 427.57 nM, strongly suggests that the F201L mutation adversely affects the protein–ligand interaction. Furthermore, the reduction in hydrogen bond formation from three in the wild-type complex to one in the mutant complex indicates a disruption of the critical hydrogen bonding network at the protein–ligand interface. The changes in van der Waals, hydrogen bonding and dissolvation energy (Vdw_hb_dissolvation energy) suggest alterations in the overall interaction profile between the protein and ligand. While an increase in electrostatic interactions was observed in the mutant complex, this increase was insufficient to compensate for the overall decrease in binding energy. These findings provide computational evidence supporting the hypothesis that the F201L mutation is deleterious and likely contributes to the pathogenesis of the associated diseases. However, further experimental validation is necessary to conclusively establish the functional consequences of this mutation. Various studies have reported an increase in the IL-6 level in OP [[208]101, [209]118, [210]119]. This particular mutation, rs2069849 was reported to increase the risk of osteoporosis in a study in Chinese postmenopausal women [[211]120]. From docking studies, it can be inferred that the less stable interleukin-6 (IL-6) produced may have difficulty in binding to its ligand and consequently different immune responses and biological processes may be adversely affected. If the IL-6 produced is quickly degraded, it can lead to insufficient immune responses and disrupt various signaling pathways, especially the Ras pathway and the JAK/STAT pathway, which are involved in regulating immune responses, inflammation and cell proliferation [[212]90]. Less stable IL-6 would results in reduced activation and regulation of these pathways. IL-6 regulates the production of anti-inflammatory and proinflammatory signals, the balance of which is disrupted if IL-6 is unstable. Its angiogenesis inducing effect may also be compromised. However, many studies have reported an increased plasma levels of IL-6 in OP patients. Increased levels may induce dysregulated angiogenesis and inflammation. It can be inferred that there may be two different conditions of IL-6, prevailing in OP patients. One is, less stable IL-6 leading to reduced signaling events and the other one, increased levels of IL-6 leading to abnormal/dysregulated angiogenesis. In a study it was shown that this polymorphism (SNP) in IL6 is not associated with circulating IL-6 levels [[213]121]. Another possibility is that the less stable IL-6 that is ineffective in binding to the ligand may be circulating in the system which leads to the detection of higher levels of this protein in patients. This SNP of IL6 was also correlated with the highest serum levels of IL-6, was linked with the worst symptoms in COVID-19 and was suggested to be used as a biomarker to identify patients with worst symptoms of COVID-19 [[214]122]. The same SNP is also associated with a chronic reduction in energy in women with breast cancer [[215]123]. The IL6 SNP rs2069849 was observed to have a protective effect against immunoglobulin-A deficiency [[216]124]. Reactome and KEGG pathway enrichment analyses revealed that immune signaling and angiogenesis/endothelial function are disrupted in OP. Signaling by cytokines, various interleukins, FGFRs, estrogen, growth hormone, PI3K/AKT, BMP, insulin, leptin, NTRKs, Receptor Tyrosine Kinases, TGF-beta, IGF1R and WNT are significant signaling pathways associated with OP. The complete sets of Reactome and KEGG pathways are given in the supplementary data. Transcription pathways mediated by the TFs - FOXO, SMADs, RUNX2, RUNX3 and HIF-1 were found to be significant. RUNX2 and HIF-1 were also identified among the hub genes. The enrichment of RUNX genes highlights the importance of the Wnt/beta catenin pathway in the pathophysiology of OP [[217]125]. The Wnt signaling pathway can stimulate the expression of RUNX2 and in this way is involved in maintaining the bone mineral density [[218]126]. The TFs identified in this study are among the several TFs associated with OP in Caucasian women [[219]127]. TFs can be utilized as an effective targets for modulating metabolic disorders. Cytokine signaling, signaling by DAP12, extra-nuclear estrogen, interleukins, TCF, Wnt, PI3K/AKT, platelet, activin, erythropoietin, VEGF and transcriptional pathways involving RUNX2 and FOXO were enriched in the context of PMOP. Inflammatory cytokines such as IL-6, IL-1 and TNF alpha, which are expressed in a greater amounts in COPD patients, are considered as one of the factors leading to development of OP [[220]6] in such patients. These studies reveal the role of inflammatory pathways in OP. The signaling pathways enriched in this study shows that inflammatory/immune environment, angiogenesis, transcriptional control, kinase action and hormonal control of the body, if altered, may lead to the development of OP. This metabolic change can be ameliorated by restoring altered conditions. Recognizing where the alteration occurred, is the key to normalizing the disturbed bone remodeling. Most of the Reactome and KEGG pathways unique to the PMOP gene set (Fig. [221]9) were reported to be invariably related to immune signaling as well as to angiogenesis/endothelial function. Platelet activation plays an important role in angiogenesis and GP1b-IX-V are key platelet receptors [[222]128]. Activin, which is stored in the bone matrix, plays an important role in regulating angiogenesis [[223]129, [224]130]. In inflammatory angiogenesis, Toll-like receptors reportedly have a regulatory function [[225]131]. The VEGF, RHO-GTPases and erythropoietin pathways are essential for angiogenesis [[226]94, [227]132, [228]133]. DAP12 plays an integral role in the differentiation of osteoclasts as well as haematopoiesis [[229]134]. While normal angiogenesis is integral to bone fracture healing as well as maintaining normal BMD [[230]135], pathogenic angiogenesis is detrimental to bone function [[231]104]. The pathways - Natural killer cell mediated cytotoxicity, Hematopoietic cell lineage, B cell receptor signaling pathway, Fc epsilon RI signaling pathway and Fc gamma R-mediated phagocytosis are coming under immune system pathways in KEGG. These findings indicate that PMOP is intricately associated with the immune status of the body. Many recent studies in mice have linked angiogenesis with bone formation. Peng et al., 2016 attributed impaired angiogenesis to decreased bone formation in streptozotocin-induced osteoporotic mice [[232]136]. The role of vascular components in OP has been established in HIF signaling pathway activated ovariectomized mice [[233]137]. Moreover, Notch signaling disruption in mice links angiocrine factors with osteogenesis [[234]138]. Significant hub genes and pathways have revealed the fact that angiogenesis is a crucial factor in the development of OP. Angiogenesis–osteogenesis coupling is an important function for the maintenance of BMD. the continuous resorption and formation of bone require proper angiogenesis to maintain the BMD [[235]139]. The hub gene IL6 can stimulate neovascularization and angiogenesis. In various studies, particularly in cancer, IL-6 has been linked to the induction of angiogenesis and VEGF [[236]140–[237]142]. In rheumatoid arthritis, IL-6 is known to induce angiogenesis where VEGF also plays a major role. The literature suggests that IL-6 induces angiogenesis, as well as osteoclast formation and inflammation, thus leading to bone resorption [[238]143, [239]144]. While angiogenesis is important for osteogenesis in bone remodeling, it can be speculated that a defective angiogenesis may be detrimental, since in most cases of OP, the levels of IL-6 are high [[240]145]. Several studies have associated IL-6 with defective angiogenesis. In glucocorticoid–induced osteoporotic mice, the oversecretion of IL-6 inhibited β-catenin activity and induced defective osteogenesis. Insulin-like growth factor 2 is found to be essential for sprouting angiogenesis and gene expression knockdown studies in human umbilical vein endothelial cells revealed that IGF binding proteins 3 and 4 regulate sprouting angiogenesis [[241]73]. Extracellular matrix (ECM) components play key role in angiogenesis [[242]146]. Proteoglycans in the ECM play important roles in biomineralization and bone formation [[243]75]. Studies in ECM proteoglycan deficient (biglycan and decorin) mice have shown that they are involved modulating osteoblastogenesis from bone marrow stromal cells [[244]74]. Activins play a role in bone formation and have been proposed for therapeutic purposes because of their ability to control matrix mineralization [[245]78]. In an inflammatory environment, activin A has been shown to reduce vasculogenesis, according to studies in endothelial cells and perivascular cells exposed to inflammatory stimuli [[246]79]. The Glycoprotein (GP)Ib-IX-V complex, which is expressed on the surface of platelets are involved in inflammatory complications in addition to its primary role in hemostasis and thrombosis [[247]77]. DAP-12 is involved in osteoclast formation, the overexpression of which can increase osteoclastogenesis [[248]81], and it is known to mediate the microglial inflammatory response resulting in neuroinflammation in the central nervous system [[249]82]. Ras signaling regulates vascular endothelial growth factor and can stimulate endothelial cells and various other pathways resulting in angiogenesis [[250]91]. VEGF signaling plays a pivotal role in angiogenesis [[251]94]. Bone formation, angiogenesis and fracture healing was noted in erythropoietin treated mice with femoral segmental defects [[252]147]. Physiological and pathological angiogenesis can be mediated by erythropoietin according to the signaling molecules present in the tissue, particularly VEGFs and VEGFRs [[253]87]. Aldosterone inhibits angiogenesis through its action on VEGFR [[254]97]. By producing proangiogenic factors, natural killer cells are reported to be involved in physiological and pathological angiogenesis mainly in tumor growth and in uterine tissue [[255]79]. The axon guidance molecules semaphorins and ephrins are integral in establishing functional coupling between bone and the nervous system to maintain bone homeostasis [[256]148]. Several of these pathways, which play crucial roles in bone physiology, are also involved in angiogenesis, although most of the studies that illustrate this angiogenic potential are carried out in tissues other than bone, mostly in tumor growths. It can be inferred by comparing the molecular pathways with existing literature that the disruption of normal angiogenesis, either increased or decreased or in defective form can lead to a decrease in BMD. Various studies have shown that an increase in angiogenesis results in bone formation and protection from OP [[257]137, [258]149]. Several biomolecules viz. IL-6 which induces angiogenesis were reported to be high in osteoporotic patients as stated above. Thus, it can be rounded up that there is a close association between angiogenesis and the incidence of OP through the different molecular pathways and hub genes. Whether it is the increased, decreased or defective angiogenesis that leads to OP needs to be confirmed with more experimental studies. Several genes and signaling pathways control or influence angiogenesis, the balance of which may direct normal angiogenesis. Disruption of angiogenesis can be attributed to changes in hormone levels as well as physical activity in osteoporotic subjects, mainly postmenopausal women and aged people, which is reflected in altered molecular pathways and gene expression. Estrogen deficiency during menopause can stimulate the pathway GnRH secretion, as a result of which FSH and LH increase in the circulation, a condition that subsides with aging [[259]150, [260]151]. BMD can be negatively affected by increasing FSH levels, since it induces osteoclastogenesis and bone resorption [[261]152]. FSH levels in plasma decrease with untrained mechanical loading of the body in women [[262]153]. In some forms of mechanical loading, a decrease in FSH levels has also been noted in males [[263]154]. It seems that untrained physical activity could keep in check the increasing FSH levels in the plasma after menopause. Estrogen deficiency can lead to increased production of sclerostin, which is a bone tissue specific inhibitor of the Wnt-beta catenin pathway [[264]116]. Additionally, there is a difference in the expression of sclerostin from osteocytes according to the loading of bones with mechanical strain [[265]155]. Sclerostin has also been shown to be associated with angiogenesis in in-vitro studies [[266]156]. Continuous excessive exercise, which is exhausting, makes the body prone to an inflammatory state and leads to organ damage [[267]157]. Loading of bones leads to the release of IL-11 from osteoblasts and osteocytes, through which Wnt signaling is upregulated and osteoblastogenesis is promoted, as well as enabling the modulation of adipogenesis [[268]158]. Adipocyte dysfunction can be modulated, or its appearance can be delayed by exercise/physical activity of the body in postmenopausal women [[269]159]. In a study of postmenopausal women, physical activity was found to be the parameter most closely related to hip fracture [[270]151]. Estradiol levels as well as BMD were found to be increased in postmenopausal women who carried out regular anaerobic exercise [[271]160]. An increase in bone mass during adolescence, which influences the peak bone mass attained by an individual, is strongly influenced by the mechanical loading of bone which acts via estrogen receptor alpha [[272]161]. A study of bone defect healing in mice revealed that exercise improves angiogenesis which in turn accelerates bone healing [[273]162]. Ovariectomized mice show improved bone mineral density while angiogenesis is induced along with osteogenesis [[274]137]. Thus, loading of bones via moderate mechanical strain may help prevent the development of PMOP to a certain extent. Immune signaling can promote and regulate angiogenesis and endothelial function. They can also modulate angiogenesis. Fracture healing requires a balanced immune signaling leading to normal angiogenesis [[275]25]. Microfractures that occur in bone need to be repaired to maintain normal BMD [[276]163]. Immune signaling leading to chronic inflammation can result in endothelial dysfunction. Endothelial dysfunction is considered as a predictor of the development of OP in postmenopausal women [[277]164, [278]165]. Improper immune signaling can lead to the release of certain proinflammatory cytokines that can damage endothelial cells. Moreover, dysregulated angiogenesis can lead to endothelial dysfunction and vice versa. Endothelial cells play a pivotal role in angiogenesis by responding to angiogenic signals such as VEGF and angiocrines. Moreover, various immune cell signaling pathways result in the production of angiogenic factors such as VEGF and FGF, and contribute to the formation of blood vessels at the site of tissue injury. Proinflammatory cytokines such as TNF-α, IL-1β and IL-6 can induce the production of angiogenic factors [[279]94, [280]166, [281]167]. Thus immune signaling and angiogenesis/endothelial function are interconnected. The immune pathways are either elicited by inflammatory cytokines or have an inflammatory function in addition to having a role in angiogenesis. The inflammatory environment in the body can adversely affect bone remodeling resulting in osteoclast predominance as evidenced by studies in COVID-19 patients, where there is a proinflammatory cytokine surge. This has also proven with in-vitro studies and in mouse models [[282]168]. Various studies have confirmed the independent role of inflammation in OP [[283]169]. Angiopoietin 2, which is released after endothelial activation can mediate angiogenesis as well as inflammation through various cytokines viz. TNF, IL-1 and VEGF [[284]166]. Glucocorticoids are used to suppress inflammation, and glucocorticoid (GC)-induced osteoporosis (GCOP) is the most common cause of iatrogenic osteoporosis [[285]170]. The long-term use of GCs promotes the differentiation and maturation of osteoclasts and inhibits osteoblasts leading to OP [[286]171]. Therefore, when glucocorticoids are used for the treatment of inflammation associated with COVID-19, they can lead to the development of iatrogenic OP. The findings from this study can lead to a better understanding of the disease progression of OP. These findings can lead to the discovery of novel treatments and diagnostic measures. The visualization and quantification of angiogenesis in bone tissue could be applied as an early diagnostic method utilizing hub genes and molecular pathways. Manipulation of molecules involved in angiogenesis, particularly angiocrines can lead to improved treatment strategies for increasing bone density. Personalized treatment methods according to a patient’s angiogenic profile could be identified through advanced research in this area. This study can have a significant impact on bone regenerative medicine. The identified molecular pathways identified could be utilized to promote vascularization in bone tissue engineering. High-risk individuals could be motivated to promote healthy angiogenesis via lifestyle interventions such as physical activity, exercise and a balanced diet as a preventive strategies. Conclusion The present study identified 28 genes, ACTB, AKT1, ALB, GAPDH, HIF1A, IL1B, IL6, INS, JUN, MAPK3, STAT3, TGFB1, TNF, TP53, CTNNB1, EGFR, ESR1, MYC, PPARG, SRC, IGF1, SPP1, VEGFA, RUNX2, TLR4, BGLAP, BMP2 and SOX9, that play critical roles in the etiology of osteoporosis. Pathways mediated by the transcription factors FOXO, SMADs, RUNX2, RUNX3, HIF-1, and TP53 were enriched. Thirteen hub genes and the KEGG pathways identified as unique to PMOP are directly associated with immune signaling and angiogenesis/endothelial function. Reactome pathway analysis reinforces the role of angiogenesis. The substitution F201L in IL6 was identified as the most deleterious SNP associated with the hub genes of OP. Mechanical loading of bones, which can influence angiogenesis and in turn the immune environment of bone, can be suggested as a preventive strategy in postmenopausal osteoporosis. The role of angiogenesis in OP etiology could be further established by systems analysis and cell culture/animal trials. Inflammation plays a direct role in OP and when glucocorticoids are used to treat inflammation, it can accelerate the occurrence of OP. Hence, better options for modulating inflammatory immune pathways could include antibodies and antagonists such as interleukin-1 receptor antagonists. Supplementary Information [287]Supplementary Material 1.^ (603.5KB, xlsx) Acknowledgements