In-silico analysis predicts disruption of normal angiogenesis as a causative factor in osteoporosis pathogenesis

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.

Introduction

Osteoporosis (OP) is a serious disorder of bone metabolism. The symptoms of this condition are fragility fractures and low bone mineral density (BMD) [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 [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 [3]. The health-related quality of life of postmenopausal women with OP is affected to variable extents [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 [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% [6].

Normal bone mass is maintained in the body via bone remodeling via signals such as mechanical loading, microfractures, stress, and hormones [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 [8]. Understanding the molecular signals that regulate osteoclasts and osteoblasts is essential in the development of drugs and treatment strategies [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 [10] and the Wnt signaling pathway is integral for bone formation [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 [12]. Type H vessels in the bone marrow are a specialized capillary subtype that is involved in the coupling of angiogenesis and osteogenesis [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 [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 [15]. The impairment of angiogenesis due to changes in the ECM caused by cellular senescence during aging [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 [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 [17]. Teriparatide, a synthetic hormone, human recombinant PTH (1–34) is used to increase bone mineral density in PMOP patients [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 [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 [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 [21]. The intestinal flora and their metabolites can be used to monitor patients with PMOP [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 [9]. The monoclonal antibody denosumab can lead to hypocalcemia and ONJ [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 [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 [9].

The molecular etiology of OP is complex and has not been fully established [24]. Several metabolic, immune, cell cycle and other pathways are perturbed during the progression of OP [25, 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 [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 [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 [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 [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 [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 [30, 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. 1).

Fig. 1

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 (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 (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 (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 (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: 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. 2. When pooled, the total number of hub genes identified for OP and PMOP was 28.

Fig. 2

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 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. 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

Fig. 3

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 2 and 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

Table 3 Key functions of 20 hub genes associated with OP according to GeneCards, NCBI-Gene, UniProt, Enrichr GO molecular function and literature mining

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 4. The position of the SNP rs2069849 was predicted as a conserved region by ConSurf as shown in Fig. 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. 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 5; Fig. 6.

Table 4 Results from gnomAD, I-Mutant2.0, Mupro, ConSurf and COACH

Fig. 4

Results from the evolutionary conservation prediction tool ConSurf server showing the conservation score of the amino acids of the IL-6 protein

Fig. 5

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

Fig. 6

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. 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. 8).

Fig. 7

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

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. 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 6.

Fig. 9

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

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 [98] and the drug tanshinone [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 [100, 101]. Signal transducer and activator of transcription 3 (STAT3) plays a role in the regulation of bone remodeling via various signaling cytokines [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 [103, 104]. VEGFA stimulates angiogenesis in bone tissue and enhances the formation, differentiation and function of osteoblasts, resulting in an increase in BMD [105]. RUNX2 encodes a nuclear protein that functions as a TF and confers osteogenic effects [106] and it is involved in the vascularization of bone tissue via the induction of angiogenesis along with HIF-1 alpha [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 [108,109,110,111]. Osteoblasts secrete the bone gamma-carboxyglutamate protein encoded by BGLAP, also called osteocalcin (OC), which has an important role in bone remodeling [112] and angiogenesis [113]. BMP2 has a role in the epigenetic control of OP [114] and its deviant signaling is detected in OP patients [115]. It plays an important role in osteogenesis–angiogenesis coupling [116, 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 2 and 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 [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 [101, 118, 119]. This particular mutation, rs2069849 was reported to increase the risk of osteoporosis in a study in Chinese postmenopausal women [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 [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 [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 [122]. The same SNP is also associated with a chronic reduction in energy in women with breast cancer [123]. The IL6 SNP rs2069849 was observed to have a protective effect against immunoglobulin-A deficiency [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 [125]. The Wnt signaling pathway can stimulate the expression of RUNX2 and in this way is involved in maintaining the bone mineral density [126].

The TFs identified in this study are among the several TFs associated with OP in Caucasian women [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 [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. 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 [128]. Activin, which is stored in the bone matrix, plays an important role in regulating angiogenesis [129, 130]. In inflammatory angiogenesis, Toll-like receptors reportedly have a regulatory function [131]. The VEGF, RHO-GTPases and erythropoietin pathways are essential for angiogenesis [94, 132, 133]. DAP12 plays an integral role in the differentiation of osteoclasts as well as haematopoiesis [134]. While normal angiogenesis is integral to bone fracture healing as well as maintaining normal BMD [135], pathogenic angiogenesis is detrimental to bone function [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 [136]. The role of vascular components in OP has been established in HIF signaling pathway activated ovariectomized mice [137]. Moreover, Notch signaling disruption in mice links angiocrine factors with osteogenesis [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 [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 [140,141,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 [143, 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 [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 [73]. Extracellular matrix (ECM) components play key role in angiogenesis [146]. Proteoglycans in the ECM play important roles in biomineralization and bone formation [75]. Studies in ECM proteoglycan deficient (biglycan and decorin) mice have shown that they are involved modulating osteoblastogenesis from bone marrow stromal cells [74]. Activins play a role in bone formation and have been proposed for therapeutic purposes because of their ability to control matrix mineralization [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 [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 [77]. DAP-12 is involved in osteoclast formation, the overexpression of which can increase osteoclastogenesis [81], and it is known to mediate the microglial inflammatory response resulting in neuroinflammation in the central nervous system [82]. Ras signaling regulates vascular endothelial growth factor and can stimulate endothelial cells and various other pathways resulting in angiogenesis [91]. VEGF signaling plays a pivotal role in angiogenesis [94]. Bone formation, angiogenesis and fracture healing was noted in erythropoietin treated mice with femoral segmental defects [147]. Physiological and pathological angiogenesis can be mediated by erythropoietin according to the signaling molecules present in the tissue, particularly VEGFs and VEGFRs [87]. Aldosterone inhibits angiogenesis through its action on VEGFR [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 [79]. The axon guidance molecules semaphorins and ephrins are integral in establishing functional coupling between bone and the nervous system to maintain bone homeostasis [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 [137, 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 [150, 151]. BMD can be negatively affected by increasing FSH levels, since it induces osteoclastogenesis and bone resorption [152]. FSH levels in plasma decrease with untrained mechanical loading of the body in women [153]. In some forms of mechanical loading, a decrease in FSH levels has also been noted in males [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 [116]. Additionally, there is a difference in the expression of sclerostin from osteocytes according to the loading of bones with mechanical strain [155]. Sclerostin has also been shown to be associated with angiogenesis in in-vitro studies [156]. Continuous excessive exercise, which is exhausting, makes the body prone to an inflammatory state and leads to organ damage [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 [158]. Adipocyte dysfunction can be modulated, or its appearance can be delayed by exercise/physical activity of the body in postmenopausal women [159]. In a study of postmenopausal women, physical activity was found to be the parameter most closely related to hip fracture [151]. Estradiol levels as well as BMD were found to be increased in postmenopausal women who carried out regular anaerobic exercise [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 [161]. A study of bone defect healing in mice revealed that exercise improves angiogenesis which in turn accelerates bone healing [162]. Ovariectomized mice show improved bone mineral density while angiogenesis is induced along with osteogenesis [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 [25]. Microfractures that occur in bone need to be repaired to maintain normal BMD [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 [164, 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 [94, 166, 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 [168]. Various studies have confirmed the independent role of inflammation in OP [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 [166]. Glucocorticoids are used to suppress inflammation, and glucocorticoid (GC)-induced osteoporosis (GCOP) is the most common cause of iatrogenic osteoporosis [170]. The long-term use of GCs promotes the differentiation and maturation of osteoclasts and inhibits osteoblasts leading to OP [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.