Bioinformatics analysis of oxidative phosphorylation-related differentially expressed genes in osteoporosis

Abstract

Background

Osteoporosis (OP) is a common metabolic bone disease characterized by decreased bone mass and increased fracture risk. Recent studies suggest that oxidative phosphorylation (OXPHOS) plays a crucial role in the pathogenesis of OP. This study aims to investigate the differential expression and potential functional roles of OXPHOS-related genes in OP.

Methods

We downloaded gene expression data from two OP-related datasets, GSE56815 and GSE7429, using the GEOquery package. We also collected OXPHOS-related genes from the GeneCards and MsigDB databases. The limma package was used for differential expression analysis of GSE56815, and differentially expressed genes (DEGs) were identified. We intersected these DEGs with OXPHOS-related genes to identify OXPHOS-related differentially expressed genes. Functional enrichment analyses, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG), were conducted using the clusterProfiler package. Additionally, we performed gene set enrichment analysis (GSEA). The Mann–Whitney U test analyzed differences in the expression of OXPHOSRDEGs, and their diagnostic potential was assessed by Receiver Operating Characteristic (ROC) curves. Correlation analysis, Protein–Protein Interaction (PPI) network construction, mRNA-miRNA, mRNA-TF interaction network construction, and immune infiltration analysis using CIBERSORT were also conducted. RAW264.7 cells were induced in vitro for 3 days to differentiate towards osteoblasts, and RT-PCR assay was used to verify the differentiation and detect the differential expression of target genes.

Results

Our results identified 31 DEGs in GSE56815, with 26 upregulated and 5 downregulated genes. Among these, we identified 10 OXPHOSRDEGs: VPS35, TBC1D2, UBQLN2, SH3GLB2, WWP1, NFKBIA, MFSD10, SLC2 A3, RP2, and ZNF91. GO and KEGG enrichment analyses revealed significant involvement of these genes in mechanisms such as the positive regulation of the protein catabolic process and the endocytosis pathway. ROC analysis demonstrated high diagnostic accuracy for VPS35 (AUC = 0.832) and TBC1D2 (AUC = 0.751). Correlation analysis indicated strong relationships between certain OXPHOSRDEGs. The PPI network highlighted 8 hub genes with significant functional similarity among them.

Conclusion

This study systematically elucidates the differential expression and potential mechanisms of OXPHOS-related genes in OP through comprehensive bioinformatics analyses. The identified key genes offer valuable insights into the molecular underpinnings of OP and present potential diagnostic biomarkers and therapeutic targets for further investigation.

Introduction

OP is a systemic disease of bone metabolism characterized by decreased bone mass, increased bone fragility, and increased fracture risk [1]. This condition represents a significant public health issue, threatening millions of people worldwide, particularly menopausal women and the elderly [2]. Primary osteoporosis is the most common type of osteoporosis (OP), accounting for approximately 80–90% of all cases of osteoporosis. It mainly includes postmenopausal osteoporosis (Type I) and senile osteoporosis (Type II). According to the International Osteoporosis Foundation, approximately 200 million people globally suffer from osteoporosis globally, with an estimated 8.9 million fractures annually attributed to the condition. The morbidity and mortality associated with osteoporotic fractures are substantial, with hip fractures alone resulting in a 20–30% increase in mortality within one year. Current diagnostic methods have limitations in predicting fracture risk and do not fully capture the complex pathophysiology of OP [3]. Thus, a new biomarker and therapeutic target are needed to improve the diagnosis and management of osteoporosis.

OXPHOS is a crucial metabolic pathway responsible for ATP production in mitochondria and plays an important regulatory role in neurological diseases and tumors [4]. It has been shown that mitochondrial oxidative stress induces bone formation disorders by causing apoptosis of osteoblasts, leading to bone metabolism disorders, which in turn contribute to the development of osteoporosis [5]. Moreover, some scholars have identified OXPHOS genes associated with bone mineral density (BMD) and fracture risk, further confirming the close relationship between osteoporosis and mitochondrial oxidative stress. However, the mechanisms of OXPHOS-related genes in osteoporosis remain unclear and require further research and exploration.

In this study, we conducted an analysis of gene expression data from the GEO database for primary osteoporosis by using differential expression analysis and functional enrichment analysis methods. Additionally, we constructed the PPI networks of these differentially expressed genes and performed immune infiltration analysis to further explore their possible molecular mechanisms in the development of primary osteoporosis. Our findings deepen the understanding of OXPHOS-related genes in primary osteoporosis and provide potential biomarkers and therapeutic targets for bone metabolism disorders.

Materials and methods

Data download

We downloaded the OP-related datasets GSE56815 [6] and GSE7429 [7] from GEO database [8]. The data platform for the GSE56815 dataset is GPL96. This dataset includes 80 samples from Homo sapiens, consisting of 40 OP patients and 40 Control patients, with all samples included in the analysis. The GSE7429 dataset, also from Homo sapiens and on the GPL96 platform, includes 20 samples: 10 OP patients and 10 Control samples (Table 1). The merged datasets (GSE56815 and GSE7429) used were all from the same chip platform and employed the same GPL96, which indicates technical compatibility and facilitates our integrated analysis. OXPHOSRGs were collected from GeneCards [9] database (https://www.genecards.org/) and MsigDB database (http://www.gsea-msigdb.org/gsea/msigdb/human/search.jsp). Using"Oxidative Phosphorylation"as the keyword, we obtained OXPHOSRGs, and then we used “Protein Coding” as the screening criterion. This resulted in 570 OXPHOSRGs from the GeneCards database. Additionally, we obtained, 3511 OXPHOSRGs from the MsigDB database. The oxidative phosphorylation-related genes obtained in the above way were merged and deduplicated to obtain 3807 OXPHOSRGs.

Table 1 GEO dataset information list

Analysis of DEGs

The limma package [10] was used to perform differential expression analysis on the GSE56815 dataset, identifying the DEGs in OP between different groups. The criteria for DEGs were |logFC|> 0.25 and p_adj < 0.05. We intersected the GSE56815 DEGs with the OXPHOSRGs to identify OXPHOS-related DEGs (OXPHOSRDEGs) for further analysis. A volcano plot was created to illustrate the differential analysis results, and a heatmap was generated to show the OXPHOSRDEGs.

GO and KEGG enrichment analysis

GO analysis which includes biological process (BP), molecular function (MF), and cellular component (CC) categories, was conducted to understand the potential role of genes in an organism[11]. KEGG analysis was a bioinformatics method for resolving and visualizing the function of genes or proteins in an organism and the metabolic pathways[12]. The R package clusterProfiler (for enrichment analysis, version: 4.10.0) was used for both GO and KEGG enrichment analyses and calculated the enrichment significance with the hypergeometric test, with statistical significance criteria set at p < 0.05 and FDR < 0.25. Additionally, we performed Benjamini-Hochberg (BH) correction on the p-values to control the false discovery rate (FDR).

GSEA

GSEA assesses the contribution of DEGs to specific phenotypes by analyzing the enrichment of gene sets in a list of genes ordered by phenotypic relatedness [13]. Genes in the GSE56815 dataset were sorted according to logFC, and all genes involved in the differential analysis were enriched using the clusterProfiler package. The GSEA parameters were set as follows: seed of 2022, 5000 permutations, each gene set containing more than 10 genes and fewer than 500 genes. The gene sets"c2.all.v2022.1.Hs.symbols.gmt [all Canonical Pathways] (3050)"from the MSigDB database for Homo sapiens were used, with significant enrichment criteria of p_adj < 0.05 and FDR < 0.25. We chose the"c2.all.v2022.1.Hs.symbols.gmt"gene set because it is derived from the MSigDB database and covers a wide range of biological pathways, which is suitable for the molecular mechanism research of osteoporosis. This gene set is helpful for identifying bone metabolism-related signaling pathways and enhancing the systematicness and biological significance of the analysis. To control for the confounding factors in the pathways, we adopted the GSEA method and combined with the adjusted p-values to screen out significant pathways to reduce random deviations.

ROC analysis and differential expression analysis of OXPHOSRDEGs

Differences in OXPHOSRDEG expression were analyzed using the Mann–Whitney U test. The R package ggplot2 was used to create group comparison plots for the differential analysis. ROC curves were then plotted for OXPHOSRDEGs to evaluate their diagnostic performance. The ROC curve evaluates the performance of binary classifiers by plotting the true positive rate against the false positive rate [14]. ROC curves for OXPHOSRDEGs in diseases were drawn using the R package pROC, and the diagnostic performance was evaluated through the area under the curve (AUC). The closer the AUC is to 1, the better the diagnostic performance.

Correlation analysis of OXPHOSRDEGs

The correlation of the expression of differential genes in the GSE56815 osteoporosis dataset was analyzed using the Spearman algorithm. The top two differentially expressed genes were displayed through a correlation heatmap. The correlation scatter plot was generated using the R package ggplot2. Similarly, the correlation heatmap of differential genes from the GSE7429 dataset was drawn, and the top two differentially expressed genes with the strongest correlation were selected. The correlation scatter plot for these genes was also created using the R package ggplot2. A larger correlation coefficient indicates a stronger correlation.

The CIBERSORT algorithm was used for the immune infiltration analysis

Using the CIBERSORT [15] algorithm in conjunction with the LM22 feature gene matrix, we filtered out immune cells with enrichment scores greater than 0. The specific results of the immune cell infiltration in the osteoporosis dataset GSE56815 were obtained and displayed on a scaled bar graph. The correlation between immune cells was calculated using the Spearman algorithm, and the correlation analysis was shown in a correlation heatmap using the R package heatmap. The Spearman algorithm was also used to calculate the correlations between immune cells and model genes, retaining results with p < 0.05. The correlations between model genes and immune cells were displayed using correlation bubble plots, created with the R package ggplot2.

PPI network and functional similarity analysis

PPIs consist of individual proteins interacting in various life processes. Systematic analysis of PPI in biosystems helps understand how proteins work in biosystems, the biological signaling response machinery, and energy and material metabolism under specific physiological conditions and the function relationships among proteins. The interactions between known and predicted proteins often extend beyond the STRING database [16]. The PPI network related to OXPHOSRDEGs was constructed by the STRING database, with a minimally required interaction score >0.150 [17].

The semantic comparison of GO annotations provides a quantitative method to calculate the similarity between genes and genomes, forming the basis for many bioinformatics analyses. The GO semantic similarity of hub genes was calculated using the GOSemSim package [18]. The results of functional similarity profiling were displayed using the ggplot2 package.

The GeneMANIA [19] database was used to generate hypotheses about gene function and perform functional gene analysis. For a particular query gene, GeneMANIA can identify genes that might share functions based on their interactions. We predicted genes with similar functions to hub genes using the GeneMANIA online tool and downloaded the interaction network. In the resulting figure, the inner circle represents the DEGs from our study, the outer circle represents genes with similar functions, and the color of the connecting lines indicates their related functions.

Construction of mRNA-miRNA and mRNA-TF interaction network

We predicted miRNAs interacting with OXPHOSRDEGs using the ENCORI database [20]. The screening criterion for mRNA-miRNA interaction pairs was pancancerNum > 10, and the mRNA-miRNA interaction network was visualized using Cytoscape software.

Thousands of binding base sequence matrices and their binding sites were identified from DNA-binding protein ChIP-seq data, and transcriptional regulatory relationships between millions of transcription factors (TFs) and genes were predicted by the CHIPBase database (version 3.0) [20] (https://rna.sysu.edu.cn/chipbase/). TFs binding to hub genes were identified from the CHIPBase database. We screened mRNA-TF interaction pairs using the criterion that the sum of samples found upstream and downstream was more than 3 and visualized the results using Cytoscape software.

Quantitative validation of hub genes using RT-PCR

Dulbecco’s Modified Eagle Medium (DMEM) containing 10% Fetal Bovine Serum (FBS) was used to incubate RAW264.7 cells for 1 day. The control group continued to be cultured for 3 days, while the experimental group was induced with DMEM containing 33 ng/mL Macrophage Colony-Stimulating Factor (M-CSF), 10% FBS, and 50 ng/mL receptor activator of nuclear factor-κB ligand (RANKL) for 3 days. The medium was changed daily. Total RNA from both groups was extracted using an RNA rapid extraction kit and reverse transcribed to cDNA using a reverse transcription kit. RT-PCR detection was performed with the SYBR Green PCR Master Mix, and all reactions were set up with 4 replicate wells. Relative quantification was conducted using the 2−ΔΔCq method.

Statistical analysis

All data processing and analysis in this study were conducted using R software (version 4.2.2). Unless otherwise specified, the independent samples T-test (Student’s t-test) was used for the comparison of normally distributed continuous variables between two groups. The Mann–Whitney U-test was used to analyze differences between non-normally distributed variables. The Kruskal–Wallis test was used for multiple group comparisons. Spearman’s correlation analysis was used to calculate the correlation coefficients between different molecules. Unless specified otherwise, all statistical p-values were double-sided.

Results

Technology roadmap

See Fig. 1.

Fig. 1

Comprehensive workflow and key findings in osteoporosis research. OP osteoporosis, DEGs differentially expressed genes, OXPHOSRDEGs oxidative phosphorylation related differentially expressed genes, GSEA gene set enrichment analysis, GO gene ontology, KEGG Kyoto encyclopedia of genes and genomes, ROC receiver operating characteristic curve, PPI Protein–protein interaction network, TF transcription factors

Dataset correction

Firstly, we normalized the datasets GSE7429 and GSE56815 using the R package sva to obtain the normalized datasets GSE7429 and GSE56815. We then compared the datasets before and after standardization using distribution boxplots (Fig. 2A–D). The distribution boxplot results show that the standardized datasets GSE7429 and GSE56815 effectively eliminate deviations or experimental conditions caused by technical differences, resulting in non-biological variation.

Fig. 2

Corrected data set. A Boxplot of the GSE7429 dataset before correction. B Boxplot of the corrected GSE7429 dataset. C Boxplot of the GSE56815 dataset before correction. D Boxplot of the corrected GSE56815 dataset. Light green represents the Control group, and light orange represents the OP group

Analysis of DEGs

The R package limma was used to obtain DEGs of the dataset GSE56815 for osteoporosis disease between the osteoporosis and control groups. The results show that the GSE56815 dataset contains 31 DEGs (Fig. 3A). To obtain OXPHOSRDEGs, all DEGs obtained from the GSE56815 dataset were intersected with OXPHOSRGs, resulting in a total of 10 OXPHOSRDEGs (VPS35, TBC1D2, UBQLN2, SH3GLB2, WWP1, NFKBIA, MFSD10, SLC2 A3, RP2, ZNF91) for our subsequent analysis. A Venn diagram was plotted to visualize this intersection (Fig. 3B). Based on the Venn diagram results, a heatmap of the differential expression of the 10 OXPHOSRDEGs in the dataset GSE56815 was plotted (Fig. 3C), with gene names arranged in descending logFC order.

Fig. 3

Analysis of DEGs. A Volcano plot showing the differential analysis between the OP group and the Control group in the GSE56815 dataset. B Venn diagram showing the intersection of OXPHOSRGs and differential genes in the GSE56815 dataset. C Heatmap of the differential expression of OXPHOSRDEGs in the GSE56815 dataset. In the heatmap, light blue represents the Control group, and light orange represents the OP group

GO and KEGG enrichment analysis

The 10 OXPHOSRDEGs (VPS35, TBC1D2, UBQLN2, SH3GLB2, WWP1, NFKBIA, MFSD10, SLC2 A3, RP2, ZNF91) were analyzed for GO gene function enrichment to understand the BP, MF, and CC involved (Table 2). These 10 OXPHOSRDEGs were mainly enriched in BP such as regulation of protein catabolic process, positive regulation of protein catabolic process, vacuole organization, regulation of macroautophagy, and positive regulation of cellular protein catabolic process. In terms of MF, they were enriched in organic anion transmembrane transporter activity, GTPase activator activity, anion transmembrane transporter activity, cadherin binding, and GTPase regulator activity. KEGG enrichment analyses showed significant enrichment of the 10 OXPHOSRDEGs in the endocytosis KEGG pathway (Table 2). We presented the GO and KEGG enrichment analysis results using bar graphs (Fig. 4A) and bubble plots (Fig. 4B). Additionally, network diagrams for BP (Fig. 4C), MF (Fig. 4D), and KEGG (Fig. 4E) were drawn based on the gene function enrichment analysis. The connecting lines indicate the labeling of corresponding molecules and entries.

Table 2 GO/KEGG enrichment analysis results

Fig. 4

GO and KEGG analysis. A, B Enrichment analysis of OXPHOSRDEGs by GO and KEGG is shown in bar graphs and bubble plots. C, D Enrichment analysis results for OXPHOSRDEGs by GO are shown in network diagrams. E Enrichment analysis results for OXPHOSRDEGs by KEGG shown in network diagrams. In the network diagram, the red dot represents the specific pathway, and the blue dot represents the specific gene

Gene set enrichment analysis

To determine the impact of gene expression levels on disease occurrence, we analyzed the association between gene expression in the GSE56815 dataset and the involved biological processes, affected cellular components, and performed molecular functions using GSEA. Based on the enrichment results, pathways with significant enrichment were screened, and the top 4 pathways with the largest NES were selected to show the results. For the selected results, we created mountain maps (Fig. 5A) and GSEA classical enrichment maps (Fig. 5B–E). In the GSE56815 dataset, genes were significantly enriched in neutrophil degranulation (Fig. 5B), the IL-4 signaling pathway (Fig. 5C), calcineurin-regulated NFAT-dependent transcription in lymphocytes (Fig. 5D), and IL4-mediated signaling events (Fig. 5E) pathways (Table 3).

Fig. 5

GSEA of the GSE56815 dataset. A Mountain map of four main biological functions. B–E Genes were significantly enriched in neutrophil degranulation, the IL-4 signaling pathway, calcineurin-regulated NFAT-dependent transcription in lymphocytes, and IL4-mediated signaling events pathways

Table 3 GSEA results in GSE56815

ROC analysis and differential expression analysis of OXPHOSRDEGs

10 OXPHOSRDEGs (VPS35, TBC1D2, UBQLN2, SH3GLB2, WWP1, NFKBIA, MFSD10, SLC2 A3, RP2, ZNF91) in the GSE56815 dataset were analyzed using the Wilcoxon rank sum test, and the differential expression results are presented as a comparison plot between groups (Fig. 6A). The expression of all DEGs in the GSE56815 dataset was statistically significant between the two groups. We then plotted the ROC curves according to the expression of these 10 OXPHOSRDEGs in the GSE56815 dataset (Fig. 6B–F). The results showed that VPS35 (AUC = 0.832, Fig. 6B), TBC1D2 (AUC = 0.751, Fig. 6B), UBQLN2 (AUC = 0.752, Fig. 6C), SH3GLB2 (AUC = 0.747, Fig. 6C), WWP1 (AUC = 0.739, Fig. 6D), NFKBIA (AUC = 0.758, Fig. 6D), MFSD10 (AUC = 0.740, Fig. 6E), SLC2 A3 (AUC = 0.726, Fig. 6E), and RP2 (AUC = 0.744, Fig. 6F) had moderate diagnostic accuracy for OP, whereas ZNF91 (AUC = 0.688, Fig. 6F) exhibited lower diagnostic accuracy for OP. The specific AUCs are detailed in Table 4.

Fig. 6

ROC analysis and differential expression analysis of OXPHOSRDEGs. A Group comparison plot of OXPHOSRDEGs in OP and Control groups. B–F OXPHOSRDEGs: ROC curves of VPS35 and TBC1D2 (B), UBQLN2 and SH3GLB2 (C), WWP1 and NFKBIA (D), MFSD10 and SLC2 A3 (E), RP2 and ZNF91 (F) between different groups

Table 4 The area under the curve

Correlation analysis of OXPHOSRDEGs

To explore the association among the differential genes, we conducted a correlation analysis of the expression of differential genes in the GSE56815 dataset using the Spearman algorithm, and a heatmap of correlations of differential genes was drawn (Fig. 7A). The top 2 differential genes with the strongest correlation were screened, and the correlation scatterplot was plotted using the R package ggplot2 to show the most highly correlated gene pairs (Fig. 7B, C). The results showed a strong positive correlation between the hub genes SH3GLB2 and MFSD10 (r = 0.816, p < 0.001) (Fig. 7B), and a moderate positive correlation between SH3GLB2 and NFKBIA (r = 0.650, p < 0.001) (Fig. 7C).

Fig. 7

Different genetic correlation analysis. A Correlation heatmap of differential genes of the osteoporosis dataset GSE56815. B, C Scatter plot of the correlation of differential genes among SH3GLB2 and MFSD10 (B), SH3GLB2 and NFKBIA (C) in the osteoporosis dataset GSE56815. D Correlation heatmap of the correlation of differential genes of the osteoporosis dataset GSE7429. E, F Scatter plot of the correlation among DEGs ZNF91 and VPS35 (E), SLC2 A3, and RP2 (F) in the osteoporosis dataset GSE7429

The correlation heatmap of the differential genes was plotted according to the correlation analysis of the expression of differential genes in the osteoporosis dataset GSE7429 (Fig. 7D). Finally, the most highly correlated top 2 differential genes were selected and the correlation scatterplot was created using the R package ggplot2 to show the top 2 gene pairs with the highest correlation (Fig. 7E, F). According to the results, there was a moderate positive correlation between the hub genes ZNF91 and VPS35 (r = 0.620, p = 0.004) (Fig. 7E), and SLC2 A3 and RP2 (r = 0.583, p = 0.008) (Fig. 7F). Correlation coefficients with an absolute value of less than 0.3 are considered to be weakly correlated or uncorrelated, between 0.3 and 0.5 are considered weakly correlated, between 0.5 and 0.8 are moderately correlated, and above 0.8 are strongly correlated.

CIBERSORT algorithm immune analysis

The infiltration abundance of 22 immune cells in the osteoporosis dataset (GSE56815) was calculated using the CIBERSORT algorithm. Based on the results of the immune infiltration analysis, the proportion of immune cells in the osteoporosis dataset (GSE56815) was calculated (Fig. 8A). The results showed that the infiltration abundance of 13 immune cells in the samples of the GSE56815 dataset was not equal to zero. These cells were: Naïve B cell, Effector memory CD4+ T cell, Gamma-delta T cell, Natural killer cell, Central memory CD4+ T cell, M0 Macrophage, M2 Macrophage, Activated dendritic cell, Mast cell, Eosinophil, Regulatory T cell (Treg), Monocyte, and Plasmacytoid dendritic cell. The correlation of abundance results of these 13 immune cell infiltrates was analyzed using correlation heatmaps of immune infiltrates (Fig. 8B). The results showed that Central memory CD4+ T cells and Natural killer cells had the greatest positive correlation (r = 0.48), while Mast cells and Monocytes had the maximum negative correlation (r = −0.84). Finally, the correlation between the 10 OXPHOSRDEGs (VPS35 TBC1D2, UBQLN2, SH3GLB2, WWP1, NFKBIA, MFSD10, SLC2 A3, RP2, ZNF91) in the osteoporosis dataset (GSE56815) and the abundance of immune cell infiltration was shown using correlation bubble plots (Fig. 8C). These plots demonstrated that SH3GLB2, TBC1D2, and MFSD10 were positively associated with T cell regulatory factors among the 10 OXPHOSRDEGs (r > 0.0, p < 0.05). Additionally, SH3GLB2 and MFSD10 were significantly negatively correlated with Gamma-delta T cells (r < 0.0, p < 0.05).

Fig. 8

GSE56815 immune infiltration analysis by CIBERSORT algorithm. A Histogram of immune cells in osteoporosis datasets (GSE56815). B Correlation heatmap of the abundance of immune cell infiltration within the osteoporosis dataset (GSE56815). C Bubble plot depicting the correlation of OXPHOSRDEGs with the abundance of immune cell infiltration within the osteoporosis dataset (GSE56815). Light green represents a negative correlation, light orange represents a positive correlation, and darker colors indicate a stronger correlation

PPI network and functional similarity analysis

We analyzed the 10 OXPHOSRDEGs (VPS35 TBC1D2, UBQLN2, SH3GLB2, WWP1, NFKBIA, MFSD10, SLC2 A3, RP2, ZNF91) for PPI using the STRING database. Only the genes that had connections to other nodes were retained and regarded as hub genes for subsequent analyses. Next, we built the PPI network consisting of 8 hub genes (VPS35 TBC1D2, UBQLN2, SH3GLB2, WWP1, NFKBIA, MFSD10, SLC2 A3) and presented it using Cytoscape software (Fig. 9A). The functional similarity of these 8 hub genes was analyzed with the GOSemSim software package. The R package GOSemSim calculated the semantic similarity among GO terms, GO term sets, gene products, and gene clusters. The results of the functional similarity analysis between the key genes were presented through box plots (Fig. 9B). The results revealed that SH3GLB2 had the highest functional similarity with other genes.

Fig. 9

PPI network and functional similarity analysis. A PPI network. B Box plot of functional similarity of hub genes. C GeneMANIA website predicts the interoperability network of hub genes with similar functions

The GeneMANIA database (http://genemania.org) was used to analyze gene functions by entering a list of genes to predict, identify, and generate gene priorities. The correlation of the 8 hub genes with other genes was analyzed using the GeneMANIA database (Fig. 9C). The results revealed that there were physical interactions and co-expression between the 8 hub genes and other genes.

mRNA-miRNA, mRNA-TF interaction network

Prediction of miRNAs interacting with 8 hub genes (VPS35 TBC1D2, UBQLN2, SH3GLB2, WWP1, NFKBIA, MFSD10, SLC2 A3) was conducted using the ENCORI database, and the mRNA-miRNA interaction network was visualized using Cytoscape software (Fig. 10A). The mRNA-miRNA interaction network consists of 7 mRNAs (VPS35, TBC1D2, UBQLN2, SH3GLB2, WWP1, NFKBIA, SLC2 A3) and 40 miRNAs, forming a total of 43 mRNA-miRNA interaction relationships. The specific mRNA-miRNA interactions are detailed in Table 5.

Fig. 10

mRNA-miRNA, mRNA-TF interaction network. A The hub genes mRNA-miRNA interaction network diagram: orange circles represent mRNAs, and light green squares represent miRNAs. B The hub genes mRNA-TF interaction network diagram: orange circles represent mRNAs, and green diamonds represent transcription factors (TFs)

Table 5 mRNA-miRNA interaction network nodes

To explore the interactions between the 8 hub genes (VPS35, TBC1D2, UBQLN2, SH3GLB2, WWP1, NFKBIA, MFSD10, SLC2 A3) and other molecules, the CHIPBase database was used to identify transcription factors binding to these hub genes. The mRNA-TF interaction network was then visualized using Cytoscape software (Fig. 10B). This network comprises 5 mRNAs (VPS35, WWP1, NFKBIA, MFSD10, SLC2 A3), 13 transcription factor molecules, and a total of 27 mRNA-TF interaction pairs (Table 6).

Table 6 mRNA-TF interaction network nodes

Validation of hub genes by RT-PCR

First, by detecting the differential expression of MMP9 in the two groups of cells, it was confirmed that the intervention group successfully induced osteoclast differentiation of RAW264.7 cells (Fig. 11A). The expression levels of the 5 hub genes (VPS35, SH3GLB2, WWP1, NFKBIA, and MFSD10) in the intervention group were significantly increased compared to the control group, while the expression levels of TBC1D2 and SLC2 A3 were significantly decreased (Fig. 11B–H). There was no statistically significant difference in the expression level of UBQLN2 between the two groups (Fig. 11I). The primer sequences, reference genes and raw Ct values for RT-PCR experiments are detailed in Table 7.

Fig. 11

Differential expression of hub genes. The expression levels of hub genes in the control and intervention groups were determined through RT-PCR, including MMP9 (A), VPS35 (B), TBC1D2 (C), SH3GLB2 (D), WWP1 (E), NFKBIA (F), SLC2 A3 (G), MFSD10 (H), and UBQLN2 (I)

Table 7 The primer sequences, reference genes and raw Ct values for RT-PCR experiments

Discussion

Osteoporosis, as a metabolic bone disease, is characterized by a decrease in bone mass and degradation of bone tissue microstructure, which increases skeletal fragility and the risk of fractures. This condition significantly affects quality of life and increases morbidity and mortality, particularly among the elderly population [21]. The economic burden associated with osteoporotic fractures is substantial, impacting healthcare systems worldwide [22]. Given the serious public health implications of OP, it is urgently necessary to gain a deeper understanding of its underlying molecular mechanisms to develop more effective diagnostic and therapeutic strategies.

In this study, we focused on the differential expression and functional analysis of oxidative phosphorylation-related genes in OP. Oxidative phosphorylation plays a crucial role in energy production and cellular metabolism and is a key metabolic pathway [23]. Previous research has indicated that disruptions in oxidative phosphorylation may contribute to various pathological conditions, including bone diseases [21]. By integrating gene expression data from GEO datasets (GSE56815 and GSE7429) with OXPHOSRGs from the GeneCards and MsigDB databases, we aimed to identify key genes and pathways involved in osteoporosis. Our findings indicate that OXPHOSRDEGs exhibit significant expression differences between osteoporosis patients and healthy controls, participating in essential biological processes and pathways involved in the disease. This study provides new insights into the molecular processes of OP and offers therapeutic targets and potential biomarkers to improve its diagnosis and treatment.

Several genes emerged as pivotal in the context of oxidative phosphorylation and bone metabolism in our investigation. VPS35, integral to the retromer complex, showed significant upregulation in OP patients, suggesting its involvement in mitochondrial function and autophagy [24, 25]. Its diagnostic potential, supported by an AUC of 0.832, underscores its promise as a biomarker and treatment target. TBC1D2, a GTPase-activating protein, also exhibited notable upregulation in OP, potentially influencing bone resorption pathways through endocytic modulation [26]. UBQLN2, known for its role in protein degradation, similarly displayed upregulation, likely compensating for increased protein misfolding in OP pathogenesis [27]. SH3GLB2, involved in endocytosis, demonstrated significant differential expression, highlighting its role in bone cell signaling and nutrient uptake pathways [28]. WWP1, an E3 ubiquitin ligase, was identified with upregulated expression, suggesting it regulates protein turnover in OP, possibly affecting signaling cascades crucial for bone metabolism [29]. NFKBIA, an NF-κB signaling inhibitor, showed increased expression, implicating its role in inflammation-mediated bone resorption [30]. MFSD10, a transporter protein, was upregulated, potentially influencing bone metabolism via organic anion transport [31]. SLC2 A3, a glucose transporter, exhibited elevated expression, suggesting its involvement in regulating cellular energy supply in bone cells [32]. RP2, involved in intracellular trafficking, displayed differential expression, highlighting its potential role in regulating signaling molecules critical for bone cell function [33]. ZNF91, a transcription factor, exhibited altered expression, potentially affecting gene expression crucial for bone metabolism [34]. These findings collectively underscore the complex molecular landscape of OP, suggesting these genes are key players in its pathogenesis and potential targets for therapeutic interventions aimed at mitigating bone metabolism disorders.

The results of this study indicate that the differentially expressed genes (DEGs) related to oxidative phosphorylation show significant differences between osteoporosis (OP) patients and the control group. This finding is not only statistically significant but also of great clinical importance. For instance, we observed that the upregulation of certain genes is associated with the increased metabolic activity of osteocytes, which may suggest that these genes play a crucial role in the occurrence and progression of osteoporosis. The upregulated genes may promote bone resorption or inhibit bone formation, thereby accelerating bone loss.

Furthermore, the changes in gene expression provide potential biomarkers for personalized treatment strategies, which are crucial for early identification of high-risk patients. For example, if these oxidative phosphorylation-related genes (OXPHOSRGs) can be intervened, it may help to slow down or reverse the progression of osteoporosis. Therefore, subsequent studies should focus on exploring the specific functions of these genes and their potential application in treatment to provide more effective solutions for the clinical management of osteoporosis.

The enrichment analysis results from our study revealed significant involvement of the identified OXPHOSRDEGs in several biological processes and pathways, particularly in the regulation of proteolytic metabolic processes and endocytosis. GO analysis showed that these DEGs were primarily concentrated in processes such as regulation of proteolytic processes and positive regulation of proteolytic processes, as well as molecular functions like organic anion transmembrane transporter activity and GTPase activator activity. The KEGG pathway analysis further highlighted the significant enrichment of these genes in the endocytosis pathway. In the future, we plan to further verify the functions of these genes and explore their specific roles in bone cell metabolism.

Endocytosis is a key lytic process involved in the internalization of exocytotic molecules and membranous proteins, essential for maintaining cellular homeostasis and signaling [35]. Dysfunction in endocytosis has been associated with various diseases, including OP, where it might affect bone remodeling and osteoclast function [36]. The enrichment of OXPHOSRDEGs in the endocytosis pathway suggests that these DEGs might play an important role in the pathogenesis of OP by influencing cellular uptake and degradation processes, potentially affecting bone density and strength.

By participating in the catabolism of proteins and the accumulation of abnormal proteins, these genes maintain cellular homeostasis and contribute to the maintenance of bone mass homeostasis [37]. Activation of the oxidative phosphorylation pathway, an important energy-generating process within the cell, not only provides sufficient energy to support various activities of osteoblasts but also promotes healthy bone growth and homeostatic maintenance. Dynamically expressed proteins play crucial roles in biological processes. These genes also play significant roles in organic anion transmembrane transport and GTPase activation, regulating the transmembrane transport of various metabolites and signaling molecules to maintain normal cellular metabolism and function [38]. The enrichment of OXPHOSRDEGs in these molecular functions suggests that they might influence cellular communication and metabolic processes in bone cells, potentially affecting bone remodeling and osteoporosis progression.

The immune infiltration analysis using the CIBERSORT algorithm revealed significant correlations among various immune cell types in OP patients. The dynamic balance among immune cells is of vital importance for bone metabolism. Specifically, naive CD4+ T cells and resting NK cells exhibited a high positive correlation (r = 0.48), whereas resting Monocytes and Mast cells showed the most substantial negative correlation (r = −0.84). Central memory CD4+ T cells activate NK cells by secreting IL-2 and IFN-γ, enhancing their cytotoxicity and immune surveillance functions. Under low inflammation conditions, NK cells inhibit osteoclastogenesis through IFN-γ. In high inflammation states, NK cells secrete TNF-α and IL-6 to promote bone resorption. Meanwhile, NK cells regulate central memory CD4+ T cells through feedback. Aging leads to dysfunction of central memory CD4+ T cells and NK cells, promoting"inflammatory aging", which in turn causes bone loss. In OP patients, the proportion of T cells increases, resulting in the predominance of pro-inflammatory factors (IL-6, TNF-α), accompanied by the decline of NK cell function, forming chronic low-grade inflammation, which inhibits osteoblast function. The interaction between NK cells and M2-type macrophages can exacerbate bone matrix degradation, especially in postmenopausal OP.

A vital component of the innate immune system, natural killer cells are known for targeting and destroying virally infected cells and tumor cells without prior sensitization [39]. It has been shown that Natural killer cells enhance the proliferation of neighboring T cells and themselves through direct cell-to-cell contact, and respond to specific Ag and CD3 cross-linking to enhance the proliferation of CD8+ and CD4+ T cells [40]. National killer cells are positively correlated with central memory CD4+ T cells, suggesting a synergistic immune response in OP patients, where both innate and adaptive immune systems are engaged. Central memory CD4+ T cells play a crucial role in initiating the adaptive immune response, indicating that antigen presentation and subsequent immune activation might be pivotal in OP [41].

Conversely, the negative correlation between resting Mast cells and Monocytes highlights a potential antagonistic relationship in the OP immune milieu. Mast cells are known for their role in allergic responses and tissue homeostasis, while Monocytes are precursors to macrophages and dendritic cells, playing a significant role in inflammation and tissue repair [42]. The observed negative correlation could imply a dysregulated immune response. Where the balance between pro-inflammatory and regulatory mechanisms is disrupted. Chronic inflammation stimulates increased osteoclast activity, accelerating bone resorption, and may also affect osteoblast function, thereby impeding bone formation, and leading to the development of osteoporosis and brittle fractures.

These immune cell interactions provide valuable insights into the complex immunological landscape of OP. The altered immune cell proportions and their interactions suggest that immune dysregulation might be a contributing factor to disease progression. Understanding these relationships could help identify novel therapeutic targets aimed at modulating the immune response to prevent or mitigate bone loss in OP patients.

RAW264.7 cells are monocyte/macrophage-like cells that originate from the Abelson leukemia virus-transformed cell line of BALB/c mice. Initially, RAW264.7 cells were described as a suitable model for macrophages because they can perform pinocytosis and phagocytosis. Later, further research demonstrated that RAW264.7 cells can respond to stimuli in vitro and subsequently generate multinucleated cells with the expected hallmark features of fully differentiated osteoclasts (RAW-OCs). RAW-OCs have been widely used in osteoclastogenesis research for over 20 years. RAW264.7 cells have multiple potential roles in immune regulation, including immune activation, enhanced phagocytic function, cytokine secretion, antioxidant and anti-inflammatory effects, and regulation of mitochondrial function. The 8 hub genes obtained from the PPI network and functional similarity analysis were validated in raw264.7 cells. RT-PCR assays confirmed that 5 hub genes were up-regulated and 2 hub genes were down-regulated in raw264.7 cells during disruption of osteoblast differentiation. But UBQLN2 showed no significant difference. This might be due to the different samples and ethnicities used, as well as the different primers employed. In vitro, the role of immune cells in OP can be further verified using immunohistochemical methods and immunoblotting techniques to verify their centralized distribution sites. Future studies should focus on elucidating the interactions of these immune cells in OP to develop targeted immunotherapy.

In summary, this study systematically explored the differential expression of genes related to oxidative phosphorylation in osteoporosis and the underlying mechanisms. Through various bioinformatics methods, we identified VPS35 and TBC1D2 as potential biomarkers for the diagnosis of osteoporosis, with AUC values of 0.832 and 0.751 respectively, indicating high diagnostic potential (AUC > 0.7). These findings suggest that these genes may serve as novel biomarkers for the early diagnosis of osteoporosis and be integrated into current clinical practice for the early screening and risk stratification of OP patients. Future studies can further validate these findings to see if they can be replicated in larger patient cohorts and explore their roles in individualized treatment strategies. Additionally, functional experiments are needed to elucidate the roles of these genes in the pathogenesis of osteoporosis, thereby providing possible directions for new therapeutic interventions. The enrichment of OXPHOSRDEGs in these biological processes and pathways provides valuable insights into their potential roles in the pathogenesis of OP. These findings suggest that these hub genes and associated pathways may offer novel therapeutic approaches to treat and prevent OP. Further research is needed to fully understand how these genes contribute to bone health and disease.

Although our study demonstrated good predictive performance on the independent validation set and the results were highly reliable and stable, we recognize that relying on whole blood samples may not reflect the bone-specific gene expression, nor the true changes in the local bone microenvironment, and may miss key bone metabolic pathways. Therefore, it is very necessary to supplement the specific data of bone tissue or bone marrow mesenchymal stem cells. Moreover, further experimental validations such as Western blot or immunofluorescence staining will be more comprehensive in confirming the expression of these genes at the protein level. We plan to consider these additional experiments in future studies to enhance the credibility and comprehensiveness of the research. Furthermore, based on the results of this study, we plan to conduct larger-scale cohort studies in the future to evaluate the applicability of these differentially expressed genes under various populations and environmental factors. At the same time, by integrating systems biology approaches to develop new biomarkers and therapeutic targets, it will help to enhance the prevention and treatment efficacy of osteoporosis.

Conclusion

In conclusion, this study comprehensively analyzes the datasets to elucidate the relationships between differential genes and the pathways involved. It describes the differential expression of OXPHOS-associated genes in OP. Furthermore, the expression trends of the hub genes were validated using RT-PCR. Additionally, potential biomarkers for OP diagnosis and treatment are identified.