Identification of Cystatin 3 as a Potential Diagnostic Biomarker for Osteoporosis Using Machine Learning

Introduction

Osteoporosis (OP) is the most common systemic bone metabolic disease, posing a serious threat to the health of middle-aged and elderly people worldwide and leading to an extremely high risk of fractures and heavy medical burdens.¹ OP is characterized by low bone mass, deterioration of bone microarchitecture, and increased bone fragility, which makes it prone to fragility fractures in the hip, spine, wrist and other parts. Common factors leading to OP include aging, estrogen deficiency in postmenopausal women, nutritional imbalance, sedentary lifestyle, chronic diseases, and long-term use of glucocorticoids, etc. At present, the treatment of OP mainly focuses on preventing fractures, with interventions using drugs such as bisphosphonates and calcitonin. Although these drugs can delay bone loss, they cannot completely reverse the already damaged bone structure, and long-term use carries certain risks of adverse reactions, such as the potential to cause osteonecrosis of the jaw and gynecological tumors.² For patients who have already suffered fragility fractures, in addition to drug therapy, rehabilitation interventions and fall prevention measures are also needed to reduce the possibility of re-fracture.³ The current difficulties and pain points in the treatment of osteoporosis are as follows: first, the diagnostic and screening measures are single, and there is a lack of early diagnostic methods; second, the treatment effect is unsatisfactory, which may easily lead to serious side effects. Therefore, there is an urgent clinical need to find new targets for the diagnosis and treatment of osteoporosis to improve the accuracy of diagnosis and the effectiveness of treatment, while reducing the occurrence of side effects.⁴

The pathogenesis of OP involves a complex process of bone metabolic imbalance, including weakened osteoblast function, enhanced osteoclast activity, disrupted bone remodeling, disordered calcium and phosphorus metabolism, and changes in the bone microenvironment. Among these, the imbalance between bone resorption and bone formation constitutes its core pathological mechanism.⁵ Reduced bone mineral density in OP patients is closely linked to the occurrence of fractures, while the deterioration of bone quality indicators such as bone microarchitecture and bone turnover rate further exacerbates the fracture risk. The progression of OP may be synergistically affected by multiple factors, including endocrine disorders, osteoblast apoptosis induced by oxidative stress, interference of immune-inflammatory factors with bone metabolism, and abnormal expression of bone metabolism-related genes regulated by genetic factors.⁶ Notably, the specific immune-related molecular mediators underlying OP progression remain incompletely understood, and few studies have systematically linked immune infiltration to potential biomarkers, leaving a critical research gap.

Additionally, in this series of complex pathological processes, the abnormal expression of Cystatin 3 (CST3) may play an important role. As a secretory protein widely present in human tissues and body fluids, CST3 belongs to the cystatin superfamily. It can participate in various physiological processes such as extracellular matrix metabolism, immune regulation, and bone metabolism by inhibiting the activity of cysteine proteases, all of which are central to the pathogenesis of OP.⁷ Studies have found that in rheumatoid arthritis (RA), CST3 is upregulated at sites of bone destruction but fails to counteract the matrix degradation effect of cathepsin B.⁸ In addition, as a serum marker for renal function testing, CST3 can assess the potential risk of hip fracture by reflecting the state of renal function,^9,10 establishing a link between CST3 and bone fragility. These studies suggest that CST3 is a key molecule regulating bone homeostasis. However, the diagnostic or prognostic value of CST3 in OP and its related roles, including the regulation of osteoblast differentiation, osteoclast activity, and its molecular mechanisms in accelerating the progression of OP by affecting bone metabolic balance, have not been fully elucidated in existing studies. Therefore, it is necessary to investigate the clinical significance and potential molecular mechanisms of CST3 in OP.

To address the aforementioned gaps, this study employs machine learning and immune infiltration analysis. Machine learning can screen specific candidate molecules from high-dimensional omics data by reducing bias and capturing complex patterns. Meanwhile, given the well-established role of immune dysregulation in OP, immune infiltration analysis is crucial, as it helps explore the potential interactions between candidate molecules and immune cells that influence bone metabolism. The main objectives of this study include clarifying the expression pattern and diagnostic potential of CST3 in OP, exploring its functional network and association with immune infiltration, and laying the foundation for CST3 to serve as a novel diagnostic biomarker for OP. This study integrates the perspectives of systems biology, immunology, and translational medicine into bone metabolism research. This strategy enhances the rigor and innovation of the research findings, addressing the limitations of previous studies that lacked systematic biomarker screening or validation.

Materials and Methods

Data Acquisition and Processing

During the search in the Gene Expression Omnibus (GEO) database, we set “osteoporosis” as the core search term, and then further narrowed down the scope to retain only “series” data and “array expression analysis” data derived from Homo sapiens. Through this search process, a total of 67 datasets were obtained. Based on the research needs of this study regarding osteoporosis, we selected four specific expression datasets from these, completed the download, and the datasets in question are GSE35956, GSE35958, GSE56116, GSE57273, and GSE230665 (Table 1). After moving into the data processing stage, we retrieved the annotation information included in the corresponding platform files for each selected dataset. Ultimately, we completed the format conversion from probe matrices to gene expression matrices, which laid a data foundation for subsequent analyses.

Table 1 Details Regarding the GEO Dataset Employed

Data Quality Control and Data Analysis

We used the limma package in R software to perform quantile normalization on the GSE35956, GSE35958, GSE56116, and GSE57273 datasets. By adjusting the expression quantiles of different samples to a common distribution, we achieved alignment of expression distributions between samples, ensuring comparability across datasets. After normalization, log2 transformation was applied to the data. Subsequently, dataset merging was performed: first, we screened for common row names (gene names) across all datasets, retaining only genes present in all four datasets; then, we merged the data column-wise (sample-wise), ultimately forming a combined expression matrix where rows represent genes and columns represent samples from all datasets. Next, the ComBat algorithm from the sva package in R software was used for batch effect correction. This method models batch effects as random variables, adjusts expression values while preserving biological variation, and is a commonly used and effective algorithm for batch effect correction in gene expression data. To ensure data comparability, we performed PCA on the pre-normalization and post-normalization datasets separately to correct for batch effects. The screened differentially expressed genes (DEGs) were visualized, including the generation of volcano plots and heatmaps for DEGs. A statistical significance threshold of p < 0.05 was set to identify DEGs and their expression changes between the OP group and the control group. Gene Ontology (GO) enrichment analysis was conducted using the “clusterProfiler” package in R software, with the generation of bubble plots and bar charts. In addition, the Metascape tool was used to perform enrichment analysis on DEGs in the datasets, as well as interaction analysis and protein-protein interaction (PPI) enrichment analysis on the enrichment terms of DEGs in the datasets, with a p-value < 0.05.

Screening Characteristic Genes

To achieve the prediction of disease status and the screening of key characteristic genes, this study employed two machine learning algorithms for analysis. Among them, the Least Absolute Shrinkage and Selection Operator (LASSO) algorithm falls into the category of regression analysis tools, and its core lies in optimizing the model through regularization technology—it not only reduces interference from redundant features but also effectively improves the accuracy of disease status prediction. The other is Random Forest (RF), a well-established technique in the field of ensemble learning. It performs particularly well in large-scale classification tasks and scenarios with skewed data, and its key advantage is reflected in the decision tree construction process: each node randomly selects a subset of features, which in turn reduces the risk of model overfitting and enhances the robustness of results. In the phase of characteristic gene screening, this study simultaneously applied two algorithms, namely LASSO regression and RF. In specific operations, we first obtained the gene sets screened by each of the two algorithms separately, then performed an intersection operation on these two sets, and identified the finally overlapping genes as the initial characteristic genes. To more intuitively present the set relationship of this screening result, a Venn diagram was subsequently plotted using the “VennDiagram” package in R.

Evaluation of the Characteristic Genes in the Dataset and Diagnostic Accuracy

The “pROC” package was used to generate initial ROC curves for the characteristic genes, with careful evaluation of the area under the ROC curve (AUC). An AUC value above 0.7 indicates high accuracy in disease diagnosis, while a value below 0.7 suggests low diagnostic accuracy. The expression levels of the characteristic genes in the GSE35956, GSE35958, GSE56116, and GSE57273 datasets were validated (training set). Meanwhile, violin plots and ROC curves of the characteristic genes in the datasets were generated, and the area under the curve was calculated. In addition, validation was performed in the new validation set GSE230665 (validation set) to identify the final candidate characteristic genes.

Establishment of OP Rat Model and Collection of Bone Marrow Mesenchymal Stem Cells

All experiments were approved by the Animal Research Ethics Committee of Stomatological Hospital, Southern Medical University, and strictly conducted in accordance with relevant guidelines and regulations, with the ethics number: NYKQ-LAEC-2025-010. A total of 8 female Sprague-Dawley (SD) rats aged 8 weeks (weighing 200–250g) were purchased from Guangdong Laboratory Animal Center. All rats were housed in a specific pathogen-free (SPF) environment and randomly divided into two groups (n=4 in each group): the control group and the OP group. To establish the OP model, rats underwent bilateral ovariectomy. Briefly, the laboratory rats were first adaptively fed for 1 week. Anesthesia was induced by inhaling a mixture of 4–5% isoflurane and oxygen, and then maintained with 1–3% isoflurane in oxygen. After shaving the abdominal area, routine disinfection with iodophor was performed 3 times. An incision of approximately 1.5–2 cm was made along the midline of the abdomen, and tissues were bluntly separated layer by layer to expose the bilateral ovaries. The adipose tissue and oviducts around the ovaries were carefully dissected. The ovarian blood vessels were double ligated with silk threads, followed by resection of the bilateral ovaries. The abdominal incision was then sutured in order. After surgery, the animals were placed in a warm and quiet environment for recovery, with adequate drinking water and regular feed provided. At 13 weeks after surgery, all rats were euthanized by intraperitoneal injection of an excessive dose of pentobarbital sodium (150 mg/kg). Four rats in each group were sacrificed, and bone tissues were collected. Primary bone marrow mesenchymal stem cells (BMSCs) were extracted from the bone tissues and incubated in a 5% carbon dioxide incubator at 37°C. BMSCs at passage 2 were used for subsequent experiments.

qRT-PCR

After collecting BMSCs from rats in the control group and OP group, total RNA was directly extracted from the cells using TRIzol reagent (AG Company, China). After determining the concentration and purity of total RNA in the samples, cDNA was synthesized using a reverse transcription kit (ABclonal, China). Finally, amplification was performed on a PCR instrument using the SYBR Green kit (Biosharp, China). For quantitative real-time polymerase chain reaction (qRT-PCR), GAPDH was used as the internal reference, and the 2-ΔΔCt method was applied to standardize the expression levels of characteristic genes. The primer sequences are shown in Table 2.

Table 2 Primer Sequences for PCR

Co-Expression Analysis and PPI Network Analysis of Characteristic DEGs

Based on the analysis, the final characteristic genes were screened out. The “corrplot” package was used for co-expression analysis to identify genes related to the final characteristic genes among the DEGs (Spearman correlation analysis), and correlation boxplots, heatmaps, and network graphs were generated. In addition, 1200 genes with a correlation (measured by the absolute value of the correlation coefficient) greater than 0.5 with the expression of the characteristic genes were screened from the GSE35956, GSE35958, GSE56116, and GSE57273 datasets, and PPI analysis (using the STRING database) was performed on these genes.

Immune Infiltration and Related Functional Analysis

The CIBERSORT (using the LM22 reference signature matrix) algorithm was used to evaluate the infiltration of immune cells in OP, and an immune cell infiltration matrix was obtained at a significance level of p < 0.05. In addition, GO enrichment analysis was performed on immune-related DEGs to assess the functional tendency of these genes.

Statistical Analysis

Bioinformatics analysis was performed using R software (version R-4.4.2). The Shapiro–Wilk test and Levene test were used to test the normality and homogeneity of variance of qRT-PCR data, respectively. The independent t-test for normally distributed variables was applied to compare data between groups, and the results were expressed as the mean ± standard deviation. Statistical significance was set at p < 0.05. All statistical analyses were conducted using IBM SPSS (version 25.0).

Results

Identification of DEGs

The obtained datasets GSE35956, GSE35958, GSE56116, and GSE57273 were subjected to standardization, log transformation, and merging. After removing batch effects from the merged data, PCA analysis was performed, and PCA plots before and after batch correction were obtained respectively (Figure 1A and B), showing the sample clustering of different dataset samples before and after removing the “batch effect”. Subsequently, 178 DEGs were identified based on the criteria that the adjusted P-value (adj.P.Val) was less than 0.05 and the absolute value of the fold change in expression (|logFC|) was greater than 1, including 108 up-regulated genes and 70 down-regulated genes (Supplementary Table 1). These results were visually presented in the form of a volcano plot and a heatmap of DEGs (Figure 1C and D).

Figure 1 Preliminary screening of DEGs. (A and B) Batch pre-normalization, pre-log-transformation, pre-merging, and post-normalization, post-log-transformation, post-merging, followed by post-batch effect removal and post-PCA for the GSE35956, GSE35958, GSE56116, and GSE57273 datasets, yielding PCA plots before and after batch correction, respectively. (C) Volcano plot of DEGs: red represents upregulated genes, green represents downregulated genes, and gray represents genes with no significant difference. (D) Heatmap of DEGs: red marks DEGs with high expression in each group, and blue marks DEGs with low expression in each group of samples.

Enrichment Analysis of DEGs and Interpretation of Core Functional Modules

We performed GO enrichment analysis on the sequences of the screened DEGs using R software to identify the common characteristics of these genes in molecular functions and biological processes. According to the results of GO analysis, the DEGs are mainly related to biological processes and other pathways, and are closely associated with functions such as regulation of extracellular matrix decomposition in bone tissue, regulation of collagen metabolic process, collagen metabolic process, negative regulation of calcium ion transmembrane transport, hyaluronic acid biosynthetic process, and negative regulation of cation transmembrane transport (Figure 2A and B). We also conducted enrichment analysis on the sequences of the screened DEGs using Metascape (http://metascape.org). The results showed that the DEGs are closely related to functions in bone tissue such as circulatory system process, regulation of collagen metabolic process, regulation of hyaluronic acid biosynthetic process, regulation of calcium ion transport, regulation of extracellular matrix decomposition, neurite development, regulation of synaptic vesicle cycle, and muscle structure development, which are basically consistent with the above results. In addition, these DEGs are also related to pathways such as compound binding, extracellular matrix-related gene sets, and TGFβ receptor pathway (Figure 2C). Figure 2D shows the interaction network of highly enriched items (displayed with CLUSTER-ID shadows, where different colors represent different enrichment pathways). Furthermore, PPI enrichment analysis was performed on the enrichment terms of DEGs in the dataset. The results showed that in the complex network, the key biological functions corresponding to the core modules (subnetworks identified by MCODE) include MCODE1 red module and MCODE2 blue module. Among them, the red module composed of EPN1-EGF-AP2A1-VAMP2 may be involved in key pathways such as endocytosis and vesicle transport. The blue module composed of HLA-DRA-CD74-CTSD may be associated with key pathways such as immunity and antigen presentation (Figure 3).

Figure 2 Functional enrichment analysis of DEGs and interaction network of highly enriched terms. (A and B) GO enrichment analysis of DEGs. (C) Metascape-based enrichment analysis of DEGs. (D) Interaction network of highly enriched terms.

Figure 3 PPI network analysis of DEGs-enriched terms and core functional modules (MCODE).

Counting and Screening of DEGs Using Two Combined Machine Learning Algorithms

The LASSO regression algorithm accurately screened out GSTM2, CST3, ITGB2, GDNF, FLJ36848, HIGD1B, MAB21L2, NFIC, CCIN, GLYATL2, MYBPH, MAGEB3, CLDN8, PYDC1, and UBR7 (Figures 4A). Meanwhile, the RF algorithm identified 8 characteristic genes, including FLJ36848, TMEM154, GSTM1, CST3, PLD3, CHCHD10, BRSK1, CARD6, PDZK1, and DYNLRB2 (Figure 4B). A comparison of the genes screened by these two methods revealed that CST3 and FLJ36848 were common overlapping genes, which constituted the initial disease characteristic genes (Figure 4C).

Figure 4 Screening of characteristic genes for osteoporosis based on machine learning algorithms. (A) Characteristic genes screened by LASSO regression algorithm. (B) Characteristic genes screened by RF algorithm. (C) Overlapping genes from the results of the two algorithms.

Intergroup Differences and Diagnostic Accuracy of Two Genes

We presented the differential expression of CST3 and FLJ36848 genes in the datasets GSE35956, GSE35958, GSE56116, and GSE57273 by plotting violin plots (Figure 5A and B), gene expression level line graphs (Figure 5C), and ROC curves (Figure 5D and E). The results showed that CST3 was up-regulated in the OP group, while FLJ36848 was down-regulated (secondary signature gene, to be discussed in the discussion section). The AUC of both genes exceeded 0.9, indicating that these two genes have high accuracy in distinguishing the disease group from the control group, providing a reference for subsequent studies.

Figure 5 Differential expression and diagnostic efficacy analysis of CST3 and FLJ36848 genes based on multiple datasets. (A and B) Violin plots of differential expression of CST3 and FLJ36848 genes in GSE35956, GSE35958, GSE56116, and GSE57273 datasets. (C) Line graphs of expression levels of CST3 and FLJ36848 genes. (D and E) ROC curves of diagnostic efficacy of CST3 and FLJ36848 genes.

CST3 Exhibits Favorable Diagnostic Efficacy in the Validation Set

The diagnostic efficacy of CST3 and FLJ36848 was validated in the validation dataset GSE230665. In the GSE230665 dataset, the ROC curve of CST3 showed an AUC greater than 0.7, while the ROC curve of FLJ36848 showed an AUC less than 0.7, which consolidated the status of CST3 as the finally identified differentially expressed gene (Figure 6A and B).

Figure 6 Diagnostic efficacy of CST3 and FLJ36848 in the validation set GSE230665. (A) ROC curve and AUC of CST3. (B) ROC curve and AUC of FLJ36848.

Expression of CST3 in BMSCs of Rats in the OP Group

Investigation shows that CST3 is used in both human and rat gene nomenclature. Subsequently, we obtained the mRNA and protein sequences of the rat CST3 gene from the NCBI database (https://www.ncbi.nlm.nih.gov/). The results of the qRT-PCR experiment showed that the expression level of CST3 mRNA (Figure 7A) was higher in the OP group than in the control group, while the expression levels of ALP and OCN mRNA (Figure 7B and C) were both lower in the OP group than in the control group, which was consistent with our initial prediction.

Figure 7 Validation of mRNA expression of CST3 and osteogenesis-related genes (ALP, OCN) in BMSCs of OP model rats. (A) mRNA expression level of CST3 in BMSCs from OP group and control group rats. (B) mRNA expression level of ALP in BMSCs from OP group and control group rats. (C) mRNA expression level of OCN in BMSCs from OP group and control group rats. (The Shapiro–Wilk test and Levene test were used to test the normality and homogeneity of variance of qRT-PCR data, respectively. The independent t-test for normally distributed variables was applied to compare data between groups. *p<0.05, **p<0.01, ***p < 0.001, and ****p < 0.0001).

Related Genes of CST3 in OP

After screening the DEGs, we extracted the expression data of CST3 and calculated the Spearman correlation between CST3 and other DEGs. Heatmaps and network graphs were used to reveal the correlation between CST3 and its related DEGs (Figure 8A and B). Subsequently, boxplots describing the expression of CST3-related DEGs in the OP group and control group of the merged dataset were generated (Figure 8C). In addition, 1200 genes with a correlation greater than 0.5 with CST3 expression were screened from the GSE35956, GSE35958, GSE56116, and GSE57273 datasets, and a PPI network was constructed. The resulting network contained 178 nodes and 584 edges, with node colors and legends indicating differences in gene connectivity. CST3 and highly connected genes showed a specific distribution in the network (Figure 9).

Figure 8 Correlation analysis and expression difference verification between CST3 and other DEGs. (A) Correlation heatmap of CST3 and related DEGs. (B) Correlation network diagram of CST3 and related DEGs. (C) Boxplots of expression differences of CST3-related DEGs between the OP group and the control group in the merged dataset.

Figure 9 PPI network of genes related to CST3 expression (1200 genes with a correlation greater than 0.5, based on GSE35956, GSE35958, GSE56116, and GSE57273 datasets).

CST3 is Associated with Immune Cells and Their Functions in OP

Immune infiltration analysis revealed differences in the infiltration proportions of immune cells between the CST3 high-expression group and low-expression group (Figure 10A). The results showed that immune cells such as Macrophages M2, NK cells, and Mast cells had a higher infiltration proportion in the CST3-High group, which may indicate that CST3 is involved in regulating the recruitment of these immune cells, thereby affecting the immune infiltration microenvironment of OP. In contrast, Dendritic cells and T cells follicular helper had a higher infiltration proportion in the CST3-Low group. Meanwhile, we also plotted stacked bar charts and heatmaps of the proportions of immune cell infiltration in the CST3 high-expression group and low-expression group (Figure 10B and C). Subsequently, we screened out immune function-related DEGs from the DEGs and performed GO enrichment analysis on these genes, which was visualized through bubble charts and bar charts (Figure 11A and B). The results showed that DEGs related to immune function were mainly involved in several categories of biological processes, including: ①extracellular matrix-related decomposition regulation, organization, and external encapsulating structure organization; ② regulatory processes and direct metabolism of collagen metabolism; ③ ion transport-related processes such as negative regulation of calcium ion transmembrane transport, amino-acid betaine transport, cation transmembrane transport, and negative regulation of monoatomic ion transmembrane transport; ④ activation of myeloid cells and other cells in immune responses; ⑤ as well as positive regulation of protein localization to the membrane, regulation of synaptic vesicle cycle, etc.

Figure 10 Analysis of differences in immune cell infiltration between CST3 high-expression and low-expression groups. (A) Differences in immune cell infiltration proportions between CST3 high-expression and low-expression groups. (B) Stacked bar charts of immune cell infiltration proportions in CST3 high-expression and low-expression groups. (C) Heatmaps of immune cell infiltration proportions in CST3 high-expression and low-expression groups.

Figure 11 GO enrichment analysis of DEGs related to immune function. (A) Bubble plot of GO enrichment analysis of DEGs related to immune function. (B) Bar Plot of GO enrichment analysis of DEGs related to immune function (The color represents the adjusted P-value (p.adjus). The redder the color, the higher the significance of the enrichment; the horizontal axis is the Count, which represents the number of genes enriched in that biological process.).

Discussion

The pathological mechanism of OP is complex and multifactorial, involving oxidative stress, estrogen deficiency, and the dysregulation of signaling pathways such as Wnt/β-catenin and RANK/RANKL/OPG.^11,12 While previous studies have identified DEGs in OP, they often rely on single-dataset analysis¹³ or focus solely on in vitro mechanisms without systematic biomarker screening.¹⁴ Our study addresses these gaps by integrating multiple datasets and adopting a machine learning strategy to robustly identify candidate genes. Specifically, we employed a machine learning approach for reliable DEG screening from the integrated multi-gene set, with further validation of diagnostic rigor via ROC analysis (AUC > 0.9) to enhance result reliability. FLJ36848 is one of the candidate genes identified through multi-dataset integration and machine learning screening. However, compared with CST3, it exhibits poorer diagnostic efficacy in the validation set with an AUC < 0.7. Additionally, functional annotations of FLJ36848 in the fields of bone metabolism and immune regulation are extremely limited in existing public databases, and no relevant mechanistic studies or experimental evidence have been reported to support its specific role in the pathogenesis of OP. Thus, we focused on investigating the correlation between CST3 and OP in subsequent studies. The CST3 gene stood out as a key candidate gene associated with bone anabolism, immune response, multiple signaling pathways, and immune cells and their functions in OP. These results have been validated in vivo via qRT-PCR in BMSCs of OP rat models and supplemented with CIBERSORT immune infiltration analysis, filling the gap in previous studies that lacked a connection between computational biology predictions in vivo validation, and the exploration of immune-related mechanisms. Additionally, we explored the functional interaction network of the candidate gene and quantified its diagnostic value. Our findings preliminarily indicate that CST3 plays an important role in the pathogenesis of OP and is expected to serve as a biomarker for the early diagnosis of OP and a potential candidate for further therapeutic investigation.

Abnormal expression of CST3 is closely related to the pathogenesis of various diseases, including Alzheimer’s disease,^15,16 breast cancer,^17,18 prostate cancer,^19,20 lung cancer,^21,22 and chronic kidney disease.²³ In addition, existing studies have shown that CST3 is involved in regulating extracellular matrix metabolism^7,24 and immune-inflammatory responses.^25–28 In the field of bone metabolism research,¹⁴ studies have indicated that CST3 may be a coupling factor derived from human osteoclasts, which can inhibit osteoclast differentiation and bone resorption function.^29,30 Recent studies have shown that USP10-mediated deubiquitination of NR3C1 can regulate bone homeostasis by controlling CST3 expression.³¹ Moreover, CST3 can be regarded as a potential marker for monitoring the response to bisphosphonate therapy and predicting osteoprotegerin levels.^32,33 Therefore, CST3 may play an important role in the regulation of bone metabolism and the efficacy detection of bone-related diseases. Although some studies have used gene expression profiles to identify key genes associated with postmenopausal osteoporosis and found through enrichment analysis that CST3 may play a role in defense response-related pathways, the specific association between CST3 and the pathogenesis of OP has not been systematically elucidated.¹³ On the other hand, previous studies have found that CST3 is closely related to immune cell infiltration and functional regulation, playing a role in the regulation of macrophage^31,34,35 behavior and T cell activity.^36–38 In this study, we advance this field by confirming CST3’s correlation with the OP-specific immune microenvironment via immune infiltration analysis, providing the direct evidence linking CST3 to the immune dysregulation characteristic of OP progression, which complements prior reports of CST3’s immune-regulatory functions in other contexts.

In this study, GO and functional enrichment analyses revealed significant differences in biological processes, including the regulation of extracellular matrix decomposition in bone tissue, the regulation of collagen metabolic processes, the negative regulation of calcium ion transmembrane transport, and the hyaluronic acid biosynthetic process. Additionally, DEGs were mainly enriched in the TGFβ receptor pathway, extracellular matrix-related gene sets, and immune-related pathways. The TGFβ signaling pathway plays a crucial role in maintaining bone metabolic balance during the development of OP. Abnormal regulation of this pathway significantly affects osteoblast differentiation, osteoclast activity, and bone matrix synthesis.³⁹ Studies have shown that the TGF-β superfamily signaling transduction module exhibits obvious expression changes in mesenchymal stem cells from patients with age-related osteoporosis.⁴⁰ A previous study demonstrated that the activation of the TGF-β1/Smads signaling pathway was inhibited in an ovariectomized OP rat model, and activating this pathway significantly promoted osteoblast viability and collagen synthesis.⁴¹ This conclusion was also verified in in vitro experiments.⁴² However, the regulatory role of TGFβ-related signaling pathways is not isolated; their association with extracellular matrix homeostasis and the immune microenvironment collectively forms a complex network underlying the occurrence and development of OP. Previous studies have reported that the imbalance between extracellular matrix decomposition and collagen metabolism is a key mechanism for the destruction of bone microstructure in OP.^43–46 Genes related to these processes among DEGs may reduce collagen degradation by inhibiting cathepsin activity, thereby maintaining bone matrix stability. Notably, the TGFβ pathway can interact with immune-related pathways to regulate the progression of OP. A study found that in an OP rat model, the expression of TNF-α and IL-17 was increased, while the expression of TGF-β and IL-10 was decreased in the colon and bone marrow.⁴⁷ Treg cells and Th17 cells originate from the same precursor cells, and the balance of their differentiation highly depends on signaling pathways regulated by TGF-β: Treg cells can maintain bone homeostasis by inhibiting osteoclast differentiation in vivo and in vitro, while Th17 cells promote osteoclast activation by secreting cytokines, exacerbating bone resorption.⁴⁸ Therefore, the downregulated expression of TGF-β in OP models may disrupt the differentiation balance of Treg/Th17 cells, and this process, together with the abnormal increase of pro-inflammatory factors, constitutes an important immunometabolic mechanism for OP progression. In addition, the polarization phenotype of macrophages is closely related to osteoclast-related mechanisms in OP patients.^49–51 Recent studies have shown that M2 macrophages secrete glutamate-containing extracellular vesicles, which alleviate OP by remodeling the fate of osteoclast precursors.⁵² Similarly, Zhang et al⁵³ demonstrated that BMSC-derived exosomes promote the alleviation of osteoporosis through M2 macrophage polarization. This suggests that CST3 may be involved in the occurrence and development of OP by regulating the TGFβ signaling pathway, modulating immune cell behavior, and maintaining the balance of collagen metabolism. Our study found that high expression of CST3 correlates with increased infiltration of M2 macrophages and enrichment of the TGFβ pathway, but this mechanism remains unconfirmed and requires further functional validation to establish causal links.

Our study also found that the PPI network analysis of DEGs between the healthy control group and the OP group revealed two core functional modules, which provide key clues for elucidating the molecular mechanism of OP. Among them, the red module composed of EPN1-EGF-AP2A1-VAMP2 focuses on endocytosis and vesicle transport pathways. The blue module composed of HLA-DRA-CD74-CTSD is associated with immune regulation and antigen presentation functions, directly pointing to the immune pathogenesis of OP. The discovery of these two modules reveals the cross-regulation between abnormal cellular material transport and immune metabolic imbalance in the pathological process of OP, providing potential research directions for targeting endocytic pathways to improve bone cell signal transduction or remodeling bone immune homeostasis through regulating antigen presentation. Combined with the results of this study, CST3, as a key characteristic gene, its high expression is associated with increased infiltration of immune cells such as M2 macrophages, NK cells, and mast cells, and has significant associations with the TGFβ pathway and extracellular matrix regulatory genes. This suggests that the occurrence and development of OP is the result of the combined action of bone cell signal disorders caused by abnormal endocytic transport, bone immune microenvironment disorders caused by imbalanced immune antigen presentation. This provides an integrated research idea for targeting CST3 to synergistically regulate cellular transport and immune balance to improve bone metabolism.

Our study results indicate that CST3 may be involved in the immune pathogenesis of OP. As a cytoplasmic regulatory factor, when upregulated in OP, CST3 may involve signaling pathways such as the TGFβ receptor pathway, regulation of extracellular matrix decomposition, and macrophage polarization. Our immune infiltration analysis results showed that the infiltration proportion of M2-type macrophages in the CST3 high-expression group increased significantly. A large number of studies have confirmed that macrophage polarization is involved in the process of OP.^51,54 Among them, M2-type macrophages can regulate the functional balance between osteoblasts and osteoclasts by secreting anti-inflammatory factors, such as IL-10 and TGF-β, thereby playing a bone-protective role.^50,55 Studies have shown that the presence of M2-type macrophages can reprogram osteoclast precursors to downregulate the expression of osteoclast-specific genes, reduce symptoms of bone loss, and delay the progression of OP.⁵² Some studies have recorded abnormal numbers and functional changes of M2-type macrophages in OP models.⁵⁶ Estrogen deficiency can mediate the differentiation of M2-type macrophages into osteoclasts, leading to changes in the M1/M2 ratio in OVX mice and promoting the progression of OP. Combined with our research results, CST3 in BMSCs of OP rats was significantly upregulated, while osteogenesis-related genes ALP and OCN were significantly downregulated. We speculate that in OP patients, the bone marrow immune-inflammatory microenvironment may stimulate BMSCs to upregulate CST3 expression, and this phenomenon is correlated with M2 macrophage recruitment and inhibition of pro-inflammatory signals. However, with the progression of OP, this correlation may weaken, failing to reverse the bone metabolic imbalance caused by OP. The upregulated expression of CST3 may be an early manifestation in OP patients, and further research is needed to explore the differences in CST3 expression levels between patients with early-stage and late-stage OP.

In this study, we analyzed the genes related to CST3 in OP, providing important clues for deciphering its functional network in the progression of OP. Through Spearman correlation analysis and network construction, we confirmed that CST3 has significant associations with a large number of DEGs. Heatmaps and network graphs intuitively present this complex molecular interaction pattern, while boxplots further verify the expression differences of CST3-related DEGs between the OP group and the control group, suggesting that these genes may synergistically participate in the occurrence and development of OP with CST3. Furthermore, in the PPI network constructed by screening 1200 genes with an expression correlation greater than 0.5 with CST3 from the dataset, the specific distribution of CST3 and highly connected genes indicates that it might be a key hub in the bone metabolism regulatory network. Considering the high expression of CST3 in OP observed in this study, as well as its association with M2 macrophage polarization and the TGFβ signaling pathway, it is hypothesized that these related genes may exert their effects through two pathways. On one hand, they participate in immune regulation together with CST3, influencing the bone immune microenvironment. On the other hand, they form interactive regulatory relationships with CST3 through extracellular matrix metabolism and osteogenesis/osteoclast-related signaling pathways, collectively regulating bone metabolism. These findings not only provide a comprehensive perspective for understanding the molecular mechanism of CST3 in OP but also lay the groundwork for screening potential synergistic regulatory targets and exploring multi-gene combined intervention strategies. In subsequent studies, functional verification can be performed on core associated genes to clarify their specific interaction modes with CST3 and the detailed roles they play in bone metabolism.

This study has several limitations. Firstly, all human data in this study were derived from the GEO database (GSE35956, GSE35958, GSE56116, and GSE57273), focusing mainly on the analysis of DEGs between OP patients and healthy individuals. However, the sample size after merging is still relatively limited. In the future, it is necessary to verify the expression patterns and diagnostic value of CST3 and related genes in larger-scale, multi-center independent datasets. Secondly, some of our results only verified the expression changes of CST3 in a rat OP model through qRT-PCR, and no functional verification has been conducted in clinical samples, such as bone marrow and bone tissue, from OP patients. Therefore, there is a lack of direct clinical evidence to support the practical application value of CST3 in OP diagnosis and treatment. Subsequent clinical transformation studies are needed to confirm its clinical significance. Thirdly, although this study revealed the association between CST3 and OP, as well as its potential connections with M2 macrophage polarization and the TGFβ signaling pathway through bioinformatics analysis, it has not deeply clarified the specific molecular mechanisms by which CST3 regulates bone metabolism, such as upstream and downstream interacting proteins and the phosphorylation status of key signaling molecules. More in vivo and in vitro functional experiments are required to decipher its functional network. Finally, this study did not consider the heterogeneity of OP. OP includes various subtypes such as postmenopausal osteoporosis, senile osteoporosis, and secondary osteoporosis, and the molecular mechanisms of different subtypes may vary. However, subtype stratification was not performed when merging multiple datasets for analysis in this study, which may mask the functional characteristics of CST3 in specific subtypes. In the future, subgroup analyses targeting different subtypes can be carried out to improve the accuracy of the research. To gain a deeper understanding of the role of CST3 in OP, our team plans to conduct more in vitro and in vivo experiments.

Conclusion

In general, through integrated bioinformatics analysis of GEO transcriptome data and animal experiments, this study preliminarily explored the potential of CST3 as a diagnostic biomarker and a candidate for further therapeutic investigation in OP, as well as its associated potential mechanisms. The key findings not only fill important gaps in the current OP research field but also hold significant implications for clinical practice and basic scientific research. Specifically, via data mining and machine learning techniques, we identified CST3, a gene closely associated with OP progression, and verified its expression pattern in BMSCs isolated from an OP rat model, confirming the association between CST3 and the hallmark features of OP. Meanwhile, bioinformatics analysis revealed the multi-dimensional potential mechanisms of CST3 in OP progression: the high expression of CST3 in OP patients is correlated with increased infiltration of M2-type macrophages and the enrichment of signaling pathways such as TGFβ, linking CST3 to the regulation of the bone immune microenvironment and the maintenance of bone metabolic balance.

The significance of this study lies in its dual contribution to addressing long-standing challenges in OP research and clinical diagnosis. In terms of diagnosis, CST3 compensates for the lack of precise detection methods for OP and shows potential as a non-invasive or minimally invasive biomarker, though its clinical application is pending further validation. For instance, CST3 can be detected in BMSCs and may be extended to blood and tissue sample detection in future studies, which is expected to resolve the current clinical dilemma where OP is difficult to identify in the pre-symptomatic stage. In terms of pathogenesis, this study reveals the critical role of CST3 as a “bridge” connecting bone immunology and bone metabolism, providing an integrated framework for understanding the complex, multi-factorial pathogenic nature of OP. It breaks through the previous isolated analysis of single pathways and clearly demonstrates how immune and metabolic processes may achieve coordinated regulation through CST3. This not only advances the development of basic research on OP but also supports CST3 as a promising candidate for further therapeutic exploration. Regulating CST3 expression may potentially enable the modulation of the bone immune microenvironment, the restoration of TGFβ signaling balance, and the stabilization of the extracellular matrix. Compared with current intervention strategies that focus on a single pathway, this study offers insights into the development of a more comprehensive therapeutic approach.

Despite the limitations of this study, these shortcomings do not diminish its research significance. The results provide further evidence for exploring the association between CST3 and the immune microenvironment as well as bone metabolism signaling pathways in OP, and for discovering new regulatory mechanisms.