FOXO4 May Be a Biomarker of Postmenopausal Osteoporosis

Introduction

Osteoporosis is a chronic skeletal disorder characterized by decreased bone strength, it should be which increases risk of fractures.¹ The main cause of osteoporosis is an imbalance between bone-forming osteoblasts and bone-resorbing osteoclasts,² which can lead to disability and morbidity, especially in the hip and spine.³ The most frequent type of osteoporosis is postmenopausal osteoporosis (PMOP), which particularly affects the elderly.⁴ The main risk factors for this disease are age, smoking, nutrition and medication,⁵ and as many as one in every two Americans over the age of 50 was affected by osteoporosis in 2020.⁶ The number of global fractures is expected to double in the next 20 years,⁷ generating a significant financial burden and public health challenge. PMOP is of particular concern because it leads to an increased risk of fractures, with subsequent negative impacts on health in older women.⁸

Table 1 The 8 Feature Genes Selected by Least Absolute Shrinkage and Selection Operator (LASSO)

The only source of osteoclast precursors consists of circulating monocytes, which can differentiate into a variety of cells, including macrophages and dendritic cell subtypes.⁹ Peripheral blood monocytes (PBMs) are able to secrete cytokines that affect osteoclast differentiation, activation, and apoptosis.¹⁰⁻¹³ The CD14/CD16 subtype of mature PBMs can migrate to sites of bone resorption, and PBMs can secrete osteoclast factors and play an important role in osteogenesis.^14–16 Moreover, PBMs could differentiate into osteoclasts.¹⁷ Therefore, PBMs have recently received much attention in studies of bone pathophysiology.¹⁸

MicroRNAs (miRNAs) are epigenetic regulators that extensively regulate gene expression, affecting multiple biological processes including bone metabolism.¹⁹ For example, miR133a in circulating monocytes has been identified as a potential marker of PMOP.²⁰ Additional miRNAs have been associated with osteoporosis^21,22 and PMOP.²³

In the present study, we employed a bioinformatics approach to identify key genes associated with PMOP, in particular potential regulators of PMOP pathology, using the Gene Expression Omnibus (GEO) database. Multi-scale embedded gene co-expression network analysis (MEGENA) and least absolute shrinkage and selection operator (LASSO) regression were performed to establish diagnostic gene signatures for PMOP and identify potential therapeutic targets.

Materials and Methods

Data Preprocessing

GSE56814, GSE13850, GSE56815 and GSE74209 were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/).²⁴ The article only collects 134 samples between high and low BMD with postmenopausal women. GSE56814, based on the GPL5175 platform, contains gene expression profiles of blood monocytes from 26 subjects with high bone mineral density (BMD) and 16 subjects with low BMD, all of postmenopausal women. GSE13850, based on the GPL96 platform, contains gene expression profiles of whole blood samples from 40 women with postmenopausal women smokers, of whom 20 had high BMD and 20 low BMD. GSE56815, based on the GPL96 platform, contains gene expression profiles of blood monocytes from 40 postmenopausal caucasian females, 20 with high BMD and 20 with low BMD. GSE74209, based on the GPL20999 platform, contains miRNA expression profiles of fresh femoral neck trabecular bone samples from 6 PMOP and 6 absence of osteoporosis with postmenopausal female subjects (according to BMD and T-score measurements [mean±  SD]: 0.882 ±0.158 and- 0.342 ±1.622, respectively). The data quality in GSE56815 and GSE13850 were evaluated using principal component analysis (PCA). To reduce batch effects from non-biological technical biases involved in the expression profiles of GSE56815 and GSE13850 as a combined dataset, the “ComBat” algorithm²⁵ and the SVA package were applied.²⁶ The workflow of the present study is shown in Figure 1.

Figure 1 Study workflow. The workflow includes four steps. Step 1: differential expression analysis was performed to obtain the differentially expressed genes (DEGs) in combined dataset (GSE56815 and GSE13850) and GSE56814, enrichment analysis was performed using DEGs; Step 2: constructing a MEGENA to obtain module genes, next, module genes with AUC > 0.60 were constructed the LASSO model to obtain feature genes; Step 3: 8 target genes were obtained that a regulated network was constructed using differentially expressed micro-RNAs (DEmiRs) in GSE74209 and feature genes; Step 4: expression profiles of combined dataset were performed ssGSEA, explored the immune cell infiltration and correlation of FOXO4 and biological pathways in PMOP. Abbreviations: AUC, area under receiver operating characteristic curve; DC, dendritic cell; MEGENA, multiscale embedded gene co-expression network analysis; LASSO, least absolute shrinkage and selection operator; logFC, log (fold change); PMOP, postmenopausal osteoporosis; ROC, receiver operating characteristic; ssGSEA, single-sample geneset enrichment analysis.

Differential Expression Analysis

Differentially expressed mRNAs (DEmRs) between low-BMD and high-BMD patients were identified in the combined dataset (GSE56815 and GSE13850) and in the dataset GSE56814 using the “limma” package in R.²⁷ Differentially expressed miRNAs (DEmiRs) between PMOP patients and controls in GSE74209 were identified using the “limma” package. DEmRs and DEmiRs associated with p < 0.05 were considered significant and analyzed further.

Enrichment Analysis

We used the “clusterProfiler” package in R²⁸ to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. P < 0.05 was defined as significant enrichment. Gene set variation analysis (GSVA) was visualized using the GSVA R package.²⁹ The “clusterProfiler” package was used to perform gene set enrichment analysis (GSEA),³⁰ and the “fgsea” package in R was used to display the results.

Construction of Co-Expression Networks

We used the “multiscale embedded gene co-expression network analysis (MEGENA)” package in R³¹ to identify a co-expression network. The first step of MEGENA analysis is fast planar filtered network (PFN) construction, and the significant gene pairs of a planar filtered network (PFN) were calculated. The constructed PFNs were input to multiscale clustering analysis (MCA) to achieve the subsequent analyses.

Construction of the LASSO Regression Model and Analysis of Receiver Operating Characteristic (ROC) Curves

Using logistic regression, we selected the optimal features from the module genes in the top three largest modules with an area under the ROC curve (AUC) > 0.6 to construct a LASSO model using the “glmnet” package.³² The most concise LASSO model was obtained, defined as the one using the fewest characteristic genes to predict the group of samples.

To further assess the diagnostic ability of the LASSO model for osteoporosis, the combined dataset was randomly divided into a training set (80%) and a test set (20%), and the model was also validated using samples in GSE56814 dataset. The ability of the LASSO model to diagnose osteoporosis was evaluated in terms of the ROC curve. Feature genes obtained from the LASSO method were extracted, and AUCs were calculated using the pROC package.³³

Prediction of Target Genes

To explore potential target genes, we selected the top three DEmiRs showing the highest and lowest changes in differential expression. The TargetScan (Http://www.targetscan.org/vert_72) database³⁴ was used to predict the potential targets of the miRNAs, and the regulatory network was constructed according to the regulation of top 3 miRNAs with the largest |log (fold change)|.

Implementation of Single-Sample Gene Set Enrichment Analysis (ssGSEA)

We used the ssGSEA method and “GSVA” package³⁰ in R to quantify the infiltration levels of immune cell types. Differences in levels of immune cell infiltration were calculated between low-BMD and high-BMD samples using the “limma” package in R. Potential correlations between immune cells or FOXO4 and immune cells were evaluated using Pearson correlation analysis. Results associated with P < 0.05 were considered statistically significant. In addition, Levels of tumor infiltration by immune cells were also evaluated using CIBERSORT (https://cibersort.stanford.edu/).

Results

Biological Processes and Pathways of DEGs in PMOP

To explore differentially expressed genes (DEGs) in PMOP patients, we used the expression profiles of a combined dataset (GSE13850 and GSE56815), hereafter referred to as “the combined dataset”. This led to the identification of 288 up-regulated and 355 down-regulated genes in low-BMD comparing with high-BMD samples (Figure 2A). PCA showed that the expression of the DEGs could be used to clearly distinguish PMOP and healthy samples in the combined dataset (Figure 2B).

Figure 2 Differential gene expression and principal component analysis. (A) Volcano plot of differentially expressed genes by comparing low/high bone density in GSE13850 and GSE56815. The red color indicates up-regulated genes, green indicates down-regulated genes, and gray indicates genes with similar expression between both groups. (B) Principal component analysis of GSE13850 and GSE56815. Abbreviation: PC, principal component.

The intersecting genes from among the DEGs in the combined dataset and from the DEGs in GSE56814 dataset were mainly involved in the KEGG pathways related to oxidative phosphorylation, fatty acid degradation and caffeine metabolism (Figure 3A). Furthermore, the common DEGs were significantly involved in biological processes related to Wnt signaling, T cell differentiation and cell-cell signaling by Wnt (Figure 3B). GSEA showed DEG enrichment in the activation of cytokine-cytokine receptor interactions, neuroactive ligand-receptor interactions, antigen processing and presentation, while DEG enrichment was linked to inhibition of Parkinson’s and Huntington’s diseases (Figure 3C).

Figure 3 Biological processes and KEGG pathway analysis. The DEGs of combined dataset (GSE13850 and GSE56815) and GSE56814 were participated in biological processes as well as pathways. (A) PMOP-associated genes mainly clustered into three signaling pathways. (B) The PMOP-related genes were predicted to be involved in five biological processes. (C) Top 10 pathways significantly activated (red) or inhibited (blue), as shown by GSEA. Abbreviations: DEGs, differentially expressed genes; PMOP, postmenopausal osteoporosis; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, gene set enrichment analysis.

Construction of a MEGENA Network and LASSO Regression of DEGs in PMOP Patients

To identify novel regulatory targets in PMOP, we constructed a MEGENA network using the intersecting genes identified above (Figure 4A). As shown in Figure 4B, we identified 66 modules and 2354 module genes. The largest module C1_4 consisted of 135 genes, followed by C1_7 with 121 genes and module C1_11 with 102 genes.

Figure 4 MEGENA of PMOP-related genes. (A) The common genes of the combined dataset and GSE56814 in the co-expression network. Each node represents a module, with the larger nodes indicating a higher number of genes. The combined dataset includes GSE13850 and GSE56815. (B) The MEGENA network showing the three largest gene modules. Each color represents one module, and triangles represent key genes in the module. Abbreviations: MEGENA, multiscale embedded gene co-expression network analysis; PMOP, postmenopausal osteoporosis.

Module genes were evaluated for their potential diagnostic roles in the combined dataset and GSE56814 based on AUC analysis (Figure 5A). The top nine identified genes (AUC > 0.60) in both datasets are shown in Figure 5B. We subsequently used the LASSO model to screen the eight feature genes with non-zero coefficients, including CARD8, FOXO4, IL1R2, MPHOSPH6, MPRIP, MYOM1, PRR5L and YIPF4 (Figure 5C, Table 1). The model gave an AUC of 0.847 in the training set (Figure 5D) and 0.812 in the test set (Figure 5E). Furthermore, GSE56814 was used to evaluate the diagnostic performance of the model, giving an AUC of 0.675 (Figure 5F).

Figure 5 LASSO model. (A) Evaluation of the diagnostic value of candidate genes involved in PMOP based on AUC. (B) The top nine genes are shown with AUC > 0.60. Red indicates up-regulated genes and blue indicates down-regulated genes. (C) Determination of the number of feature genes in the LASSO analysis. Receiver operating characteristic curves of the training set (D), test set (E) and GSE56814 (F). Abbreviations: AUC, area under receiver operating characteristic curve; LASSO, least absolute shrinkage and selection operator; PMOP, postmenopausal osteoporosis.

Correlation of DEmiRs with Biological Pathways in PMOP Patients

We next set out to explore differentially expressed miRNAs (DEmiRs) in PMOP patients and so performed a differential expression analysis. We identified 535 DEmiRs in GSE74209, with 297 that were up-regulated and 237 that were down-regulated (Figure 6A). Furthermore, the top three DEmiRs with the largest |log (fold change)| were screened, and feature genes were used to construct a regulatory network. GO functional enrichment analysis showed that regulated genes were mainly involved in the VEGFA-VEGFR2 signaling pathway and blood vessel morphogenesis (Figure 6B). Finally, we built a regulatory network involving five miRNAs targeting ACTR3, CARD8, CCDC134, FOXO4, FST, LSM5, MPRIP and SWAP70 (Figure 6C). Importantly, FOXO4, CARD8 and MPRIP were common feature genes in the LASSO model and were also identified as DEmiRNA target genes. In particular, FOXO4 had the largest AUC and therefore defined as a key gene in PMOP. We subsequently identified the predicted binding site between FOXO4 and hsa-miR-4284 using the TargetScan database (Figure 6D).

Figure 6 Regulatory network of miRNAs. (A) Volcano map of differentially expressed miRNAs comparing PMOP and healthy controls in GSE74209. Red indicates up-regulated genes, green indicates down-regulated genes, and gray indicates genes with similar expression between the two groups. (B) GO enrichment analysis of the top three largest gene modules. Each colored dot indicates a different biological process. (C) Sankey plots showing the biological processes involved in the regulatory relationships between miRNAs and mRNAs. (D) Predicted binding site of FOXO4 and hsa-miR-4284, obtained using TargetScan. Abbreviation: GO, gene ontology.

Immune Cell Infiltration in PMOP

Immune cell types were investigated in PMOP samples using the combined dataset and in GSE56814. As shown in Figure 7A, we found that dendritic cells (DCs), interdigitating dendritic cell sarcoma (iDCs) and activated DCs (aDCs) were significantly up-regulated in PMOP compared with healthy controls. Correlation analysis across 24 immune cell types in PMOP showed significant correlations in the levels of DCs, natural killer (NK) cells, T cells and cytotoxic cells (Figure 7B). As a further exploration of the relationship between immune cells and key genes in PMOP, we found a positive correlation between FOXO4 and DCs, indicating that FOXO4 and DCs may together promote the development of PMOP (Figure 7C). Moreover, we found that monocytes constituted the highest proportion of infiltrated immune cells (Figure 7D), suggesting that monocytes may play a key role in PMOP development.

Figure 7 Correlation between FOXO4 expression and immune cell abundance in PMOP. (A) Expression of differentially expressed genes and levels of immune cells in the indicated datasets. (B) Correlation between immune cell types in PMOP samples. Red represents a positive correlation, and green indicates a negative correlation. (C) FOXO4 positively correlated with DCs. (D) Estimated proportions of 20 immune cell types. Abbreviations: PMOP, postmenopausal osteoporosis; DCs, dendritic cells.

Correlation of FOXO4 and Biological Pathways in PMOP

Evaluation of the involvement of FOXO4 in PMOP indicated an AUC of 0.67 in GSE56814, and an AUC of 0.65 in the combined dataset (Figure 8A). As shown in Figure 8B, FOXO4 was up-regulated in PMOP, and we explored the top 10 biological processes and KEGG pathways that correlated with FOXO4 in PMOP. The positively correlated biological processes of FOXO4 were found to be involved in renal filtration and regulation of protein kinase activity, while the negatively correlated pathways of FOXO4 were involved in protein localization to the cytoskeleton, basement membrane disassembly and latent virus replication (Figure 8C). Furthermore, KEGG pathways of ABC transporters, cocaine addiction and nicotine addiction were positively associated with FOXO4, while the VEGF signaling pathway, beta-alanine metabolism and oxidative phosphorylation were negatively associated with FOXO4 in PMOP (Figure 8D).

Figure 8 Enrichment analysis of FOXO4 in PMOP. (A) Receiver operating characteristic curves for the indicated datasets. The combined dataset includes GSE13850 and GSE56815. (B) Expression of FOXO4 in the common genes of the combined dataset and GSE56814. (C) Correlation of FOXO4 and biological processes in immune cells. The red section indicates activation, and the green section indicates inhibition. (D) Correlation of FOXO4 and pathways in immune cells. Red indicates a positive correlation, and green indicates a negative correlation.

Discussion

PMOP is one of the most common chronic diseases, particularly in elderly women, and it has a serious detrimental impact on the body. In this study, we performed a regulated network and LASSO model analysis to screen for potential key genes in PMOP patients, and found that FOXO4 was associated with immune cell infiltration. In addition, we identified pathways related to the regulation of miRNAs in PMOP patients.

The DEGs between PMOP and healthy controls were significantly enriched in biological processes related to the Wnt signaling pathway, T cell differentiation and cell-cell signaling by Wnt. The Wnt signaling pathway is known to be involved in the regulation of bone homeostasis and may increase the risk of osteoporosis in postmenopausal women.³⁵ Precursor T cells can differentiate into Th1, Th17 and other helper T cells.³⁶ Th17 cells are involved in osteoclast differentiation and thus help maintain bone homeostasis,³⁷ Our study revealed that pathways of cytokine-cytokine receptor interactions, neuroactive ligand-receptor interactions and antigen processing and presentation may be involved in PMOP. This finding implies that natural compounds may be useful for targeting cytokine-cytokine receptor interactions in order to inhibit RANKL-induced osteoclasts and reduce the expression of genes involved in bone resorption,³⁸ For their part, neuroactive ligand-receptor interactions³⁹ as well as antigen processing and presentation by osteoclasts⁴⁰ may affect risk of osteoporosis.

MEGENA was performed to show the correlation between genes and genes in modules, and 358 module genes were used for AUC analysis. We constructed the LASSO model using the genes with an AUC > 0.6. Eight feature genes were ultimately obtained: CARD8, FOXO4, IL1R2, MPHOSPH6, MPRIP, MYOM1, PRR5L and YIPF4. CARD8 as an inflammatory sensor and is highly expressed in blood and lymphoid tissues,⁴¹ suggesting that CARD8 may promote inflammation in osteoporosis. FOXO1, FOXO3 and FOXO4 were expressed in bone cells and maintained redox balance in osteocytes, they regulate bone metabolism and participate in the pathological processes of osteoporosis.⁴² IL1R2 is specifically expressed in peripheral blood mononuclear cells of osteoclasts.⁴³ In contrast, the other feature genes that we identified here have not yet been linked to osteoporosis, so their potential role in PMOP needs to be further explored.

Construction of miRNA-regulated gene expression networks implicated genes involved in VEGFA-VEGFR2 signaling and blood vessel morphogenesis in PMOP. VEGFA can stimulate bone repair by promoting angiogenesis and stimulating pre-osteoblast proliferation to differentiate into mature osteoblasts or, intriguingly, promote osteogenesis or transformation.⁴⁴ VEGFA also modulates angiogenesis through the activation of the VEGFR1 and VEGFR2 receptors.⁴⁵ Blood vessel morphogenesis produces mature vasculature and also affects the formation and function of several important tissues.⁴⁶ Such vascularization combines with bone development and growth during osteogenesis.⁴⁷ Osteoblasts can affect bone synthesis and metabolism and promote the development of PMOP.⁴⁸ Altogether, most PMOP-associated pathways involve osteocytes as well as blood vessels, which may play a role PMOP inception.

Out of the 8 predicted target genes of DEmiRs, FOXO4 yielded the largest AUC within the LASSO model and among the target genes of miRNAs, so we hypothesize that it may be a key gene in PMOP. Immune cells are involved in osteoporosis,⁴⁹ and our study found that among the 24 immune cell types, DCs and T cells are highly expressed in PMOP, with DCs positively associated with FOXO4 expression. DCs are antigen-presenting cells that can activate and regulate T cells.^50,51 Moreover, the homeostasis of T cells can cause inflammation, resulting in substantial bone resorption and loss of bone mineral density.⁵² While FOXO4 has also been studied in bone metabolism, it has mostly been reported in the context of cancers, including leukemia, neuroblastoma, as well as cervical, colorectal, pancreatic and lung cancers.^53,54 We propose that FOXO4 may act as a key molecule in the immune environment of bone, and its further experimental exploration may open up new avenues in onco-immunology.

The biological pathways of ABC transporters, cocaine addiction and nicotine addiction were significantly associated with FOXO4 in PMOP. This is not surprising, as ATP-binding cassette g1 (ABCg1) regulates osteogenesis,⁵⁵ while nicotine addiction increases the risk of osteoporosis.^56,57 Additionally, nitrogen metabolism can partially counteract osteoporosis.⁵⁸ However, the role of cocaine addiction in osteoporosis is unclear.

We screened FOXO4 as a key biomarker of PMOP by bioinformatics and explored its role in biological processes and pathways. We found that FOXO4 is involved in immune cell infiltration, suggesting that FOXO4 may be a key molecule in the bone immune environment. However, the relevant mechanisms need to be further explored. Furthermore, our study has some limitations, especially in studies involving larger samples and animal experiments, which could externally verify and extend our purely bioinformatic analysis here.

Conclusion

Using bioinformatics analysis, we identified FOXO4 as a potential novel biomarker in PMOP, which may act as a key molecule in the immune environment of bone, and its further exploration may open up new avenues in immunology.