GEO microarray data acquisition

Osteoporosis expression profile datasets GSE7158 and GSE56815, and the BMSC-related dataset GSE105145 were obtained through the GEO database. Among them, the GSE7158 dataset included 14 normal control samples and 12 osteoporosis patient samples, the GSE56815 dataset included 20 normal control samples and 20 osteoporosis patient samples, and the GSE105145 dataset included 3 samples of BMSCs and 3 dental pulp MSC samples.

Datasets GSE7158 and GSE56815 were further merged and the ComBat function of the package sva was used to correct batch effects. All analyses in this paper were performed under R version 4.2.1.

Differential gene expression analysis

Differential analysis was performed through the “limma” package of R language, selecting the differentially expressed genes (DEGs) in each dataset with |logFC| > 1 and p < 0.05 as the screening condition. Volcano maps were drawn using the “ggplot2” package of R language and differential gene expression heat maps were drawn using the “heatmap” package.

WGCNA

To further understand the interrelationship between the merged genes from the GSE7158 and GSE56815 datasets, and to identify the central genes of the gene modules, the co-expression network of osteoporosis samples was constructed with the WGCNA algorithm. The sample cluster analysis was performed using the R language “hclust” package. Gene co-expression networks were constructed using the R language “WGCNA” package. First, the expression of individual genes was transformed into a similarity matrix based on the Pearson correlation coefficient between pairs of genes. A hierarchical clustering dendrogram was constructed using an expression matrix to detect outliers. Then, the scale-free topological fitting exponent was calculated as a function of soft threshold power based on the network topological analysis function. The topological overlap matrix-based dissimilarity measure was then used to calculate the connectivity between genes in the module with similar expression profiles, and hierarchical clustering was utilized to produce a hierarchical clustering dendrogram of genes. Finally, a dynamic mixed cutting method was used to divide these genes into different modules.

Gene selection by GeneCards analysis

Genes associated with “PMO” and “osteoblast differentiation” were searched by the GeneCards database and ranked by the Relevance score from high to low. The selected genes were imported into the STRING database for protein-protein interaction (PPI) analysis, where the “species” was limited to “human”. The PPI network of regulatory relationship was drawn through the Cytoscape software (v3.6.0).

Gene functional enrichment analysis

The analysis results were intersected using the R language “VennDiagram” package, and a Venn diagram of the intersected DEGs was obtained. GO (biological processes, cellular components, and molecular functions) and KEGG enrichment analyses were performed using the R language “ClusterProfiler” package. The cellular functions and pathways mainly affected by potential targets and key targets were analyzed. p < 0.05 was considered statistically significant.

Lasso regression analysis

Regression analysis is a mathematical model that constructs quantitative relationships between variables based on determining the dependent and independent variables. The “glmnet ()” function package in R Studio software was called to load the candidate gene matrix. When α = 1, the appropriate λ value was selected, and the internal validation was implemented by the ten-fold cross-validation method to determine the best model.

M – number of samples, k – number of parameters, λ – regularization parameters.

Receiver operating characteristic (ROC) analysis

The “pROC ()” function in the R Studio software was used to draw the ROC curve for the model. The overall survival (OS) included in the data was selected as the dependent variable. With the specificity as the horizontal axis and the sensitivity as the vertical axis, the ability of the AUC evaluation model to distinguish the risk of PMO prognosis was calculated. In the case of AUC > 0.5, the closer the AUC was to 1, the better the model could predict the PMO prognosis.

Cell culture and treatment

Human bone marrow mesenchymal stem cell line UE7T-13 (Biobw, Beijing, China) was cultured in RPMI-1640 medium (11875093, Gibco) containing 10% fetal bovine serum (FBS; 10099141, Gibco), and 1% penicillin/streptomycin (SV30010, Gibco) at 37°C under 5% CO₂. Then cells were transfected with siRNA plasmids (HanBio, Shanghai, China) targeting PIK3CA (PIK3CA-si) and KLF2 (KLF2-si) as well their corresponding NC (PIK3CA-si-NC and KLF2-si-NC) using Lipofectamine 3000 reagents (L3000015, Invitrogen). The normal cells were cultured with normal culture medium while the remaining cells were cultured with the osteogenic induction medium (#HUXMX-90021, Cyagen Bioscience Inc., Suzhou, Jiangsu, China). siRNA sequences are shown in Supplementary Table SI.

Cell differentiation determination

To induce osteoblast differentiation, BMSCs (1 × 10⁵) were cultured in a endothelial cells (ECs) were cultured on 0.1% gelatin-coated coverslips in 24-well culture plates and osteogenic induction medium (HUXMX-90021, Cyagen) for 1 week. Alkaline phosphatase (ALP) staining and alizarin red-S (ARS) staining were carried out using an Alkaline Phosphatase Assay Kit (P0321S, Beyotime, Shanghai, China) and alizarin red-S solution (C0148S, Beyotime), respectively.

RT-qPCR

Total RNA was extracted from cells and tissues via homogenization in TRIzol reagent (10296010, Thermo Fisher Scientific Inc., Waltham, MA; 1 ml per 100 mg tissues), addition of chloroform (200 ml), and centrifugation at 12000 g and 4°C for 10 min. The upper aqueous phase was collected and mixed with 500 ml of isopropanol to precipitate the RNA. RNA was dissolved in RNase-free water (10–30 ml) and quantified on a NanoDrop ND-3300 Spectrophotometer (Thermo Fisher Scientific Inc.). The RNA (1 mg) was reverse transcribed using TaqMan Reverse Transcription Kit (N8080234, Thermo Fisher Scientific Inc.), and qPCR was performed using PowerUp SYBR Green Master Mix (Thermo Fisher Scientific Inc., A25741). GAPDH and U6 were regarded as internal references and the fold changes were calculated by relative quantification (2^-ΔΔCt method). The primer sequences are shown in Supplemnetary Table SII.

Western blot

Total protein was extracted, separated and transferred onto membranes. The membrane was blocked using TBST (5 ml) with 5% skimmed milk powder (P0216-300 g, Beyotime) and underwent overnight incubation at 4°C with primary antibodies to PIK3CA (p110 & p85) (51 kDa, #ab40776, Abcam, Cambridge, UK), KLF2 (139 kDa, #ab236507, Abcam), Akt1 (56 kDa, #ab238477, Abcam), p-Akt1 (56 kDa, #ab81283, Abcam), collagen I (139 kDa, #ab138492, Abcam), ALP (39 kDa, #ab307726, Sigma), RUNX2 (57 kDa, #ab236639, Abcam), OCN (11 kDa, #ab133612, Abcam), and GAPDH (36 kDa, #PA1-987, Thermo Fisher Scientific Inc.). The next day, the membrane was incubated with HRP-labeled secondary antibody goat anti-rabbit (A0208, Beyotime). ECL (P0018FS, Beyotime) was used for chemiluminescence detection and membranes were photographed by the ChemiDoc XRS + system (Bio-Rad). The band intensity was quantified using Image J analysis software, normalized to GAPDH [25].

ChIP assay

BMSCs cultured in the osteogenic induction medium for 1 week were collected and probed with antibodies to KLF2 (139 kDa, #ab236507, Abcam) and IgG (52 kDa, #ab109489, Abcam, serving as a NC). The primers used for this assay are shown in Supplementary Table SIII [26].

Statistical analysis

The Wilcoxon test was used to compare differences between the two sets of non-normally distributed data. The unpaired t-test was used to compare differences between two sets of normally distributed data. The Pearson correlation analysis was performed to compare the correlation between the two data groups and to calculate the correlation coefficients. The R software (4.2.1), SPSS (22.0), and Prism were used for statistical analysis and mapping. p < 0.05 demonstrated that the difference was statistically significant.

GEO-based differential analysis for screening DEGs in PMO

In order to find the molecular mechanisms to alleviate PMO, differential analysis and WGCNA through GSE7158 and GSE56815 datasets were performed. In combination with the GeneCards website, Lasso regression analysis and ROC curve analysis were carried out, and the biological function of core genes was analyzed by functional enrichment. The BMSC-related dataset GSE105145 was used to screen the transcription factors involved in the regulation, and a PPI regulatory network was constructed to provide a theoretical basis for further understanding the mechanism of PMO and finding potential diagnosis and treatment targets. The bioinformatics analysis process is shown in Supplementary Figure S1.

To screen DEGs in PMO, differential analysis was performed using the combined dataset of GSE7158 and GSE56815. A total of 10419 DEGs were selected, where 8308 DEGs were downregulated, and 2112 DEGs were upregulated (Figures 1 A, B).

Figure 1

GEO-based differential analysis. A – Heat map of DEGs in the GSE7158 and GSE56815 datasets. B – Volcano map of DEGs in the GSE62020 dataset. The green dots indicate downregulated genes, the red dots indicate upregulated genes, and the black dots indicate genes with no significant difference

WGCNA-based screening for the green core module

First, 12530 genes in the expression profiles of the combined dataset of GSE7158 and GSE56815 were selected for cluster evaluation, and a total of 66 samples were used for subsequent WGCNA (Figure 2 A). The best value of the soft threshold (power) was set to 14 to build a gene expression network (Figure 2 B).

Figure 2

Data validation and screening for optimal soft thresholds. A – Heterogeneity viewed by sample clustering. B – Non-scale fitting index analysis of soft threshold. C – Analysis of soft threshold average connectivity

Modules with high similarity of module eigengenes were merged by dynamic mixed clipping method, where the minModules size was set to 30, deepSplit to the default value of 2, and cut height to 0.25. Each module was shown in a different color (Figure 3 A). Pearson correlation analysis was performed between each module eigengene and different groups (Figure 3 B). The six co-expression modules with the top absolute value of the correlation coefficient and p value were selected for subsequent analysis, including the green module (2647 genes), the black module (4413 genes), the purple module (99 genes), the magenta module (241 genes), the pink module (824 genes) and the grey module (4281 genes). Among them, the green module showed the most significant correlation (cor = 0.76, p < 1e-200) (Figure 3 C).

Figure 3

Co-expression network constructed based on the combined dataset of GSE7158 and GSE56815. A – Gene dendrogram. B – Pearson correlation analysis between each module eigengene and different groups. C – Correlation coefficient of the green module

Identification of 197 core DEGs related to PMO

GEO-based differential analysis was performed to screen PMO-related DEGs. The WGCNA screened core green module-related DEGs for Venn analysis, and 2189 common DEGs (Figure 4 A) were obtained. Through the GeneCards website, 1188 PMO-related genes were further screened to intersect with 2189 DEGs, and 197 candidate DEGs related to PMO progression were obtained (Figure 4 B).

Figure 4

Screening for core DEGs and regulatory function analysis. A – Venn diagram of the green module-related DEGs and DEGs in GEO database. B – Venn diagram of PMO-related genes and 2189 DEGs. C – KEGG pathway enrichment analysis of candidate DEGs. D–F – Functional analysis of candidate DEGs in biological processes (D), cellular components (E), and molecular functions (F)

The candidate DEGs were subjected to KEGG and GO functional enrichment analyses. The results of the KEGG pathway analysis showed that the candidate DEGs were mainly enriched in the PI3K-Akt pathway, cytokine-cytokine receptor interaction, Jak-STAT pathway, and Th17 cell differentiation pathway (Figure 4 C). The results of GO functional analysis revealed that the candidate DEGs were mainly enriched in ossification, response to peptide, organic hydroxy compound metabolic process and multicellular organismal homeostasis in biological processes. The candidate DEGs were mainly enriched in the collagen-containing extracellular matrix, extracellular matrix, blood microparticle, and secretory granule lumen in cellular components. The candidate DEGs were mainly enriched in receptor regulator activity, receptor ligand activity, growth factor activity, and cytokine receptor binding in molecular function (Figures 4 D–F).

The above results suggested that the candidate DEGs mainly played roles in ossification and cytokine regulation, and were enriched in extracellular matrix and vesicles. The molecular functions of candidate DEGs were mainly related to cytokine-related enzyme activity and receptor activity. According to the results of the KEGG analysis, the candidate DEGs were involved in the regulation of the PI3K-Akt pathway. It has been found that the activation of the PI3K-Akt pathway could regulate the bone balance of PMO, inhibit osteoblast apoptosis and promote cell proliferation [27]. Therefore, we further investigated the molecular regulatory mechanism of the PI3K-Akt pathway in PMO.

PIK3CA as the core gene of the PI3K-Akt pathway affecting PMO

To further investigate the core regulatory genes in candidate DEGs, we performed Lasso regression analysis of the 14 related genes involved in the PI3K-Akt pathway (GHR, ATF2, IL4R, COL2A1, IRS1, FGF18, JAK2, CSF1, HSP90AA1, FLT3, MAPK3, PIK3CA, CCND1, and MTOR) in the candidate DEGs. The Lasso algorithm used the L1 criterion for model construction, and the change of the regression coefficient under the L1 criterion is shown in Figure 5 A. With a parameter α = 1, ten-fold cross-validation was used for model internal validation. The selection process diagram of the cross-validation parameter λ was drawn, where the log (λ) was indicated on the horizontal axis and the root mean square error value (RMSE) on the vertical axis (Figure 5 B). When the model variable value was 4, the RMSE value was minimal; therefore, the corresponding candidate DEGs (PIK3CA, IL4R, JAK2, and MAPK3) were selected for subsequent analysis.

Figure 5

Screening for the core genes of the PI3K-Akt pathway affecting PMO. A – Lasso regression analysis of candidate DEGs and regression coefficient variation diagram under L1 criterion. B – Selection process diagram of cross-validation parameter λ. C – ROC curve diagram of 4 candidate DEGs. D – Boxplot of the expression of 4 candidate DEGs in the combined dataset of GSE7158 and GSE56815. E – Construction of a PPI network of the top 30 genes related to osteoblast differentiation retrieved from the GeneCards website and core genes of the PI3K-Akt pathway. F – Calculation of the degree value of genes in the PPI network

The ROC curves of the four DEGs were drawn, and AUC values were calculated (Figure 5 C). The results showed that the AUC values of PIK3CA, IL4R, JAK2, and MAPK3 were 0.754, 0.787, 0.793, and 0.792, respectively, all of which met AUC > 0.75, suggesting that the model had a good discriminatory ability to predict PMO prognosis. Further calculation of PIK3CA, IL4R, JAK2, and MAPK3 expression in PMO showed that PIK3CA was significantly less expressed in the disease group than in the control group. In contrast, the other three genes showed the opposite trend (Figure 5 D). Therefore, we selected PIK3CA as the core gene to explore the regulatory mechanism of the PI3K-Akt pathway in PMO.

Previous literature has shown that promoting osteoblast differentiation of BMSCs can effectively alleviate PMO [28]. In order to further explore the mechanism, we screened for genes related to osteoblast differentiation through the GeneCards website, and selected the top 30 genes and core genes of the PI3K-Akt pathway to construct a PPI network (Figure 5 E). From the network, we observed that Akt1 and PIK3CA, the core genes of the PI3K-Akt pathway, were closely related to osteoblast differentiation. Further calculation of the Degree value, which is the most direct measure to characterize node centrality in a network analysis, identified Akt1 to be one of the core genes in the PPI network (Figure 5 F). The above results indicate that the PI3K-Akt pathway can participate in the osteoblast differentiation of BMSCs.

PIK3CA promoted osteoblast differentiation of BMSCs

Having identified PIK3CA as the core gene of the PI3K-Akt pathway affecting PMO, we continued to verify the effect of PIK3CA on the osteoblast differentiation of BMSCs. The mRNA and protein expression of PIK3CA was found to be increased in the BMSCs cultured with osteogenic induction medium versus normal BMSCs (Figures 6 A, B), implying that PIK3CA may facilitate osteoblast differentiation of BMSCs. Subsequently, we designed two siRNA sequences targeting PIK3CA and transfected them into BMSCs. The PIK3CA-si2 sequence showed superior silencing efficiency and was used for subsequent experiments (Figures 6 C, D).

Figure 6

PIK3CA facilitates osteoblast differentiation of BMSCs. A – mRNA expression of PIK3CA in the BMSCs cultured with normal medium and osteogenic induction medium detected by RT-qPCR. B – Western blot of PIK3CA protein in the BMSCs cultured with normal medium and osteogenic induction medium. C – Silencing efficiency of PIK3CA-si in the BMSCs detected by RT-qPCR. D – Silencing efficiency of PIK3CA-si in the BMSCs detected by Western blot. E – Osteoblast differentiation of BMSCs transfected with PIK3CA-si detected by ALP staining and ARS staining. F – Western blot of osteoblast differentiation-related proteins collagen I, ALP, RUNX2 and OCN was decreased in BMSCs transfected with PIK3CA-si. G – Phosphorylation level of Akt1 in BMSCs transfected with PIK-3CA-si detected by Western blot. *p < 0.05, ** p < 0.01

As shown in Figure 6 E, a reduction was noted in the osteoblast differentiation of BMSCs transfected with PIK3CA-si. In addition, the expression of osteoblast differentiation-related proteins collagen I, ALP, RUNX2 and OCN was decreased in BMSCs transfected with PIK3CA-si (Figure 6 F). Meanwhile, silencing of PIK3CA inhibited the phosphorylation level of Akt1 (Figure 6 G). Altogether, these results demonstrate that PIK3CA may be a key gene to facilitate osteoblast differentiation of BMSCs.

The transcription factor KLF2 promoted PIK3CA transcription

The transcription factor is a key molecule of gene expression regulation. To explore whether there was a molecular mechanism of the transcription factor to regulate PIK3CA expression in PMO, 1369 human transcription factors were downloaded through the humantfs website to intersect with the top 100 DEGs in PMO, and only the transcription factor KLF2 was obtained (Figure 7 A). Pearson correlation analysis results showed that in the combined dataset of GSE7158 and GSE56815, transcription factor KLF2 was positively correlated with PIK3CA expression (Figure 7 B). Using the JASPAR website, the binding site of the transcription factor KLF2 (Figure 7 C) and that between KLF2 and the PIK3CA promoter region (Supplementary Table SIV) were analyzed.

Figure 7

Screening for transcription factors regulating PIK3CA expression. A – Venn diagram of DEGs in PMO and human transcription factors. B – Pearson correlation analysis of transcription factor KLF2 and PIK3CA. C – JASPAR website was used to predict the binding site of transcription factor KLF2. D – Heat map of DEGs in the BMSC-related GSE105145 dataset. E – Volcano map of DEGs in the BMSC-related GSE105145 dataset. F – Expression of transcription factor KLF2 in GSE105145. *p < 0.05 vs. the control group

BMSCs can self-renew and multi-differentiate, and they are indispensable components in osteoblast differentiation and skeletal remodeling in PMO. The process is controlled by a multistep molecular pathway, and regulated by different transcription factors and signaling proteins. To explore the expression of the transcription factor KLF2 in BMSCs, we performed a differential analysis of the BMSC-related dataset GSE105145 and selected 4510 DEGs, including 2369 upregulated genes and 2141 downregulated genes (Figures 7 D, E). We found that the transcription factor KLF2 was significantly highly expressed in the BMSC samples, compared to that in the control samples (Figure 7 F). It has been well documented that KLF2 promoted the stemness maintenance and self-renewal of BMSCs [29]. Therefore, we could draw the preliminary conclusion that BMSCs might participate in the transcription and development of PMO by regulating PIK3CA transcription through secretion of the transcription factor KLF2.

KLF2 was involved in the osteoblast differentiation of BMSCs by regulating the PI3K-Akt pathway

To verify whether the transcription of PIK3CA is regulated by KLF2, we first detected the expression of KLF2 in BMSCs cultured with osteogenic induction medium by RT-qPCR and Western blot. The results showed an increase of the mRNA and protein expression of PIK3CA in the BMSCs cultured with osteogenic induction medium versus normal BMSCs (Figures 8 A, B). Then, we designed two siRNA sequences targeting KLF2 and transfected them into BMSCs. RT-qPCR and Western blot identified their silencing efficiency, with the KLF2-si1 sequence showing higher silencing efficiency, and it was used for subsequent experiments (Figures 8 C, D).

Figure 8

KLF2 enhances the osteoblast differentiation of BMSCs by regulating the PI3K-Akt pathway. A – mRNA expression of KLF2 in BMSCs cultured with normal medium and osteogenic induction medium detected by RT-qPCR. B – Western blot of KLF2 protein in BMSCs cultured with normal medium and osteogenic induction medium. C – Silencing efficiency of KLF2-si in BMSCs detected by RT-qPCR. D – Silencing efficiency of KLF2-si in BMSCs detected by Western blot. E – Osteoblast differentiation of BMSCs transfected with KLF2-si detected by ALP staining and ARS staining. F – Western blot showing that expression of osteoblast differentiation-related proteins collagen I, ALP, RUNX2 and OCN was decreased in BMSCs transfected with KLF2-si. G – PIK3CA protein expression and phosphorylation level of Akt1 in BMSCs transfected with KLF2-si detected by Western blot. H – Binding of KLF2 to the promoter region of PIK3CA analyzed by ChIP assay. *p < 0.05, **p < 0.01.

In the absence of KLF2, osteoblast differentiation of BMSCs (Figure 8 E) was weakened in addition to decreased expression of osteoblast differentiation-related proteins collagen I, ALP, RUNX2 and OCN (Figure 8 F). Meanwhile, silencing of KLF2 inhibited the phosphorylation level of Akt1 (Figure 8 G). Further, ChIP results demonstrated that KLF2 could bind to the promoter region of PIK3CA (Figure 8 H), indicating that the transcription of PIK3CA is regulated by KLF2. Overall, KLF2 can regulate the PI3K-Akt signaling pathway through PIK3CA, thereby promoting the osteoblast differentiation of BMSCs.