Identification and validation of ubiquitination-associated genes of senile osteoporosis based on bioinformatics analysis

Background: Senile osteoporosis (SOP) is linked to the ubiquitination process, with dysregulation of ubiquitin-mediated protein turnover disrupting bone remodeling and resulting in decreased bone mineral density (BMD). This study aimed to identify biomarkers related to ubiquitination in SOP and explore their molecular regulatory mechanisms.

Methods: Transcriptomic data of SOP samples (categorized by high and low BMD) were obtained from public databases. Differential expression analysis, protein-protein interaction networks, and the CytoHubba plugin (using Maximum Neighborhood Component and degree algorithms) were utilized, alongside the Least Absolute Shrinkage and Selection Operator, to identify ubiquitination-related genes (URGs) as potential SOP biomarkers. The diagnostic potential of these biomarkers was assessed through a Support Vector Machine model and a nomogram. Their molecular mechanisms were further investigated using enrichment analysis, immune infiltration analysis, and the construction of regulatory networks. Expression levels of the biomarkers were validated in a SOP rat model, with enzyme-linked immunosorbent assay applied to detect relevant indices.

Results: RPS27A and UBE2E1 were significantly underexpressed in low BMD samples and demonstrated a strong ability to differentiate between patients with varying BMDs, making them potential diagnostic biomarkers for SOP. A positive correlation was observed between RPS27A and UBE2E1 (cor = 0.35, P = 0.026). Both genes were involved in neurodegenerative diseases, critical cellular functions, and key intracellular signaling pathways. Additionally, RPS27A showed a positive correlation with macrophages and monocytes, whereas UBE2E1 exhibited a negative correlation with T follicular helper cells (Tfh) and T helper 17 cells (Th17). The transcription factor MAX and miRNA hsa-miR-106b-5p were identified as potential regulators of both biomarkers. Western blot, immunohistochemistry, and reverse transcription quantitative PCR further confirmed significantly lower expression of RPS27A and UBE2E1 in the SOP group compared to the Sham group.

Conclusion: This study successfully identified RPS27A and UBE2E1 as key biomarkers for SOP, demonstrating their diagnostic potential and involvement in important biological pathways and immune responses, thus offering new prospects for therapeutic interventions.

1 Introduction

Osteoporosis (OP) is a chronic metabolic bone disorder that poses a significant risk of fractures and mortality in the elderly. It is typically classified into primary and secondary forms (1). Senile OP (SOP), a subtype of primary OP, is primarily caused by age-related abnormalities in bone metabolism in individuals over the age of 70 (2). The pathophysiology in older women is often linked to postmenopausal estrogen deficiency and the natural aging process, leading to the perception that the development of primary OP, particularly in older women, is inevitable (3). The pathogenesis of OP in elderly women is multifactorial. In addition to aging and estrogen deficiency, factors such as poor diet, physical inactivity, and insufficient sunlight exposure can also contribute to reduced bone mass and mineral density (4). Bone mineral density (BMD) is the most commonly used indicator of bone health, encompassing the structural integrity, metabolic transformation, cumulative bone damage (microfractures), and mineralization of bone tissue (5). BMD measurement provides a means to assess bone health status and predict the degree of OP, which is also critical for assessing fracture risk (6). While BMD remains the gold standard for diagnosing OP, its diagnostic utility is not without limitations (7). In recent years, researchers have sought to explore the mechanisms underlying SOP from perspectives such as immune system dysfunction, mitochondrial impairment, and genetic predisposition, aiming to identify novel diagnostic approaches (8–10). Therefore, the identification of new and effective biomarkers for the early diagnosis and prevention of SOP is of paramount importance.

Ubiquitination is a key post-translational modification of proteins and represents the most significant and effective protein degradation pathway in mammalian cells (11). The ubiquitination process involves three key enzymes: ubiquitin-activating enzyme (E1), ubiquitin-conjugating enzyme (E2), and ubiquitin protein ligase (E3). This process plays a critical role in the degradation of short-lived proteins and the regulation of essential biological functions, including cell proliferation, differentiation, and the growth and development of tissues and organs (12–14). Recent research has increasingly highlighted the significant role of ubiquitination in bone remodeling. OP is associated with the abnormal expression and dysfunction of ubiquitination and de-ubiquitination enzymes (15). For instance, the ubiquitin ligase Smurf1 can ubiquitinate downstream BMP molecules, leading to their degradation and inhibition of osteoblast differentiation (16, 17). Additionally, studies have shown that the ubiquitin-binding enzyme UBE2E3 influences the senescence and osteogenic differentiation of bone marrow stromal cells (BMSCs) (18). Conversely, deubiquitinating enzymes remove ubiquitin from substrate proteins, thereby modulating bone remodeling (19, 20). Hence, ubiquitin plays a pivotal role in the onset and progression of OP in the elderly, warranting further investigation into its specific mechanisms in bone reconstruction.

This study was based on the transcriptome data of “elderly postmenopausal women with different bone densities” from public databases. Through bioinformatics analysis, the ubiquitination-related genes RPS27A and UBE2E1 were selected as potential key biomarkers. The study systematically explored the diagnostic value and molecular regulatory mechanisms of these two genes in postmenopausal osteoporosis, and verified their expression levels and correlations with bone metabolism indicators in a model of elderly postmenopausal rats. To date, no research has adopted a combined strategy of bioinformatics and experimental validation to jointly identify RPS27A and UBE2E1 as biomarkers related to ubiquitination in elderly osteoporosis. Therefore, this study provides new theoretical basis for the early diagnosis, prevention strategies, and targeted treatment development of this disease.

2 Methods

2.1 Data source

Transcriptome datasets were sourced from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds). The GSE13850 dataset was derived from the GPL96 platform. It contained B-cell samples from the blood of 20 elderly postmenopausal women with osteoporosis, as well as B-cell samples from the blood of 20 age-matched controls. This dataset was used as the training cohort. The GSE56815 dataset (21) was also based on the GPL96 platform. It included monocyte samples from the blood of 20 postmenopausal women with osteoporosis and monocyte samples from the blood of 20 controls. All post-surgical samples were excluded in the study, and the dataset was ultimately used as the validation cohort. The miRNA sequencing dataset GSE64433 was derived from the GPL18402 platform. It contained whole blood samples from 3 patients with osteoporosis (OP) and 3 controls, and was mainly used for the construction of the regulatory network in this study. The URGs utilized in this study were downloaded from the Integrated Ubiquitination Site and Ubiquitin-conjugating Enzyme Database (iUUCD) (http://iuucd.biocuckoo.org/), initially containing 1,394 genes. After removing duplicates and replacing aliases with gene symbols, 1,367 unique genes were retained.

2.2 Determination of differentially expressed URGs

Differentially expressed genes (DEGs) between high and low BMD samples in the GSE13850 dataset were identified using the limma package (v3.58.1) (22). The screening threshold was set as p-value < 0.05 and |log₂ fold change (log₂FC)| > 0, and the selection of this threshold was based on the literature (23). The resulting DEGs were visualized as a volcano plot using the ggplot2 package (v3.3.6) (24) and as a heatmap using the ComplexHeatmap package (v2.4.2) (25). Subsequently, intersections between these DEGs and the 1,367 URGs were determined using the ggVenndiagram package (v1.2.3) (26), generating a set of DE-URGs. Enrichment analyses of the DE-URGs, including Gene Ontology (GO) (P < 0.05) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (P < 0.05), were conducted using the clusterProfiler package (v3.18.1) (27) to explore their biological functions and associated pathways.

2.3 Selection of candidate genes

To investigate potential protein-level interactions among DE-URGs, these genes were imported into the STRING database (http://string-db.org) to build a protein-protein interaction (PPI) network (highest confidence = 0.4). The resulting PPI network was visualized using Cytoscape software (v3.10.0) (28). The CytoHubba plugin within Cytoscape was then used, applying two algorithms—Degree and Maximal Clique Centrality (MCC)—to identify the top 10 influential genes from each algorithm. These top 10 genes were analyzed for intersection points, ultimately screening for candidate genes relevant to this study.

2.4 Recognition of biomarkers and construction of nomogram

The glmnet package (v 4.1-8) (29) was used to construct a least absolute shrinkage and selection operator (LASSO) regression model based on the candidate genes, with specific parameters set as follows: the penalty type was set as L1 regularization, and 5-fold cross-validation was adopted; during each iteration, the data was randomly divided into a training set (4/5) and a validation set (1/5) to select the optimal lambda value (lambda.min) that minimized the prediction error. Based on this optimal lambda value, feature genes corresponding to non-zero coefficients were extracted to further identify candidate genes that made significant contributions to the binary outcome. Afterwards, based on the genes screened by LASSO, the e1071 package (23) was used to construct an Support Vector Machine (SVM) model, with parameters set as follows: the radial basis kernel was selected as the kernel function (‘kernel=“radial”‘); the regularization parameter C was set to 1; and the kernel function parameter gamma was set to ‘scale=TRUE’. The diagnostic performance of the model was evaluated by plotting a receiver operating characteristic (ROC) curve and calculating the area under the curve (AUC). The performance of the classification model was further assessed by constructing a Precision-Recall (PR) curve to examine the relationship between precision and recall across varying thresholds. A decision curve was applied to assess the model’s net benefit. Additionally, the expression of these genes was analyzed within both training and validation cohorts, comparing samples with high and low BMD. Outcomes were visualized using box plots created with the ggplot2 package. Genes exhibiting significant differences and consistent expression patterns across both datasets were selected as biomarkers (P < 0.05).

Subsequently, the rms package (v 6.5-0) (30) was used to integrate the identified biomarkers and develop a nomogram model within the training cohort. Each biomarker was assigned a specific score, with the cumulative score representing the total score. This total score was then used to estimate the probability of developing SOP, with higher scores indicating a greater likelihood of SOP occurrence. A calibration curve was plotted to evaluate the model’s predictive accuracy, and a decision curve was used to assess its clinical utility. Additionally, a clinical impact curve was generated to analyze the clinical response rate of the model.

2.5 Functional and annotation analysis

Following biomarker identification, a comprehensive analysis was conducted to elucidate their intrinsic functions. In the training cohort, correlations between biomarkers were assessed based on their expression levels using the corrplot package (v0.92) (31) with the Spearman method. Gene Set Enrichment Analysis (GSEA) was then performed using the clusterProfiler package, applied to the expression matrix sorted by correlation (P < 0.05). The reference gene set used was c2.cp.kegg_legacy.v2023.2.Hs.symbols.gmt, sourced from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/). To explore the subcellular localizations of the biomarkers, each was queried in the GeneCards database (https://www.genecards.org/), focusing on entries with the highest confidence level (Confidence = 5) to determine the subcellular localization of the proteins encoded by these genes.

2.6 Immune infiltration analysis

To further evaluate the immune microenvironment in SOP patients with varying BMD, immune infiltration assays were conducted. The ImmuCellAI algorithm (32) was applied to quantify the abundance of 25 distinct immune cells in the OP and control samples from the GSE13850 dataset. Spearman correlation analysis was then applied to explore the relationships between biomarkers and immune cell variations across groups (|cor| > 0.3 and P < 0.05). This integrated approach provided insights into the immune interactions associated with SOP.

2.7 Disease association analysis and drug prediction studies

To assess the relationship between biomarkers and SOP, this study further evaluated their associations using the Comparative Toxicogenomics Database (CTD) (http://ctdbase.org), utilizing the inference score to measure the strength of these associations. In the search for potential therapeutic drugs for SOP, DGIdb (https://www.dgidb.org/) was queried for drugs targeting the identified biomarkers, and a network diagram was constructed for visualization.

2.8 Construction of regulatory network

To elucidate the molecular regulatory mechanisms of the biomarkers, regulatory networks were constructed. First, using the miRNet database (https://www.mirnet.ca) with TRRUST parameters, transcription factors (TFs) interacting with the biomarkers were predicted based on the training cohort. Differential analysis was then performed with the limma package on the GSE64433 dataset, which contained samples from three OP and three control subjects, to identify differentially expressed miRNAs (DE-miRNAs) [|log₂Fold Change (FC)| > 1 and P < 0.05]. The multiMiR package (v3.19) (33) was employed to identify miRNAs potentially regulating the biomarkers (‘table = “validated”, ‘org = “hsa”). The intersection of DE-miRNAs and predicted miRNAs with opposite expression to the biomarkers generated the terminal miRNAs, forming a miRNA-mRNA regulatory network. Additionally, lncRNAs regulating these miRNAs were predicted using the starBase database (https://rnasysu.com/encori/), leading to the construction of an lncRNA-miRNA-mRNA network. All networks were visualized using Cytoscape software.

2.9 Animal model construction and sample collection

Six female Sprague-Dawley (SD) rats, aged 8 weeks, were purchased from Beijing Huafukang Biotechnology Co., Ltd. (Production Permit No.: SCXK (Jing) 2024-0003, Usage Permit No.: SYXK (Dian) K2022-0007). After one week of adaptive feeding, the rats were randomly and evenly assigned to the Sham group (n=3) and the SOP group (n=3). The rats were housed under controlled temperature conditions, with free access to food and water. D-galactose was injected into the nape of the rats’ necks for 45 consecutive days to induce senescence. Following this, ovariectomy was performed to establish the SOP model. Body weights of the rats were recorded weekly. After 14 weeks, all rats were anesthetized, and blood was collected from the heart. The blood was allowed to stand at room temperature for 2 hours before centrifugation to obtain serum samples for enzyme-linked immunosorbent assay (ELISA) testing. Subsequently, the rats were euthanized, and their left femurs were collected for Western blot (WB) analysis and reverse transcription quantitative PCR (RT-qPCR) analysis. All procedures were approved by the Animal Ethics Committee of Yunnan Bestime Biological Technology LLC.

2.10 Hematoxylin-eosin staining and immunohistochemicalanalysis

The right femoral tissue was removed and fixed in 4% paraformaldehyde. After trimming and fixing, the tissue was dehydrated, infiltrated with wax, and embedded in paraffin. Sections were cut and mounted on slides for further analysis. H&E staining was performed to observe the pathological changes in the femoral tissue of the SD rats. IHC analysis was conducted to detect the expression of biomarkers in the femoral tissue. Each experimental technique was repeated three times. Images of the sections were captured using an inverted microscope (Nikon, TS2), and the results were analyzed using ImageJ Pro-Plus software. Statistical analysis was performed using GraphPad Prism 5 software, and intergroup differences were compared by t-test. A P-value < 0.05 was considered statistically significant.

2.11 ELISA analysis

To assess the occurrence of OP, an ELISA detection kit was used to measure the levels of osteocalcin (BGP) (Meimian: LY3497-B, China), tartrate-resistant acid phosphatase 5b (TRACP5b) (Meimian: LY3282-B, China), and osteoprotegerin (OPG) (Meimian: LY2991-B, China) in the rats’ serum. The procedure was followed according to the manufacturer’s instructions. Optical density (OD) values were measured at 450 nm using a microplate reader (BioTek: EIx800, USA) within 15 minutes. Experimental technique was repeated three times. Statistical analysis was performed using GraphPad Prism 5 software, and intergroup differences were compared by t-test. A P-value < 0.05 was considered statistically significant.

2.12 WB analysis

Proteins were extracted from the left femoral tissue using RIPA lysis buffer (Servicebio, Servicebio), and protein concentrations were measured using a BCA protein quantification kit (Beyotime, P0009). The appropriate amount of 5× protein loading buffer (Servicebio, G2013-100ML) was then added to the proteins, and the mixture was heated in a boiling water bath for 10 minutes to denature the proteins. SDS-PAGE electrophoresis was performed to separate the protein compounds, which were then transferred onto a PVDF membrane (Millipore, K2MA8350E). The membrane was blocked and incubated with the corresponding primary antibody solution overnight at 4 °C. On the following day, the membrane was incubated with the secondary antibody at room temperature for 60 minutes. The antibody dilution ratios are detailed in Supplementary Table S1. The membrane was exposed using an ECL chemiluminescent substrate (Affinity, KF8001) on a gel imaging system, and the gray values were analyzed using ImageJ software. Experimental technique was repeated three times. Statistical analysis was performed using GraphPad Prism 5 software, and intergroup differences were compared by t-test. A P-value < 0.05 was considered statistically significant.

2.13 RT-qPCR analysis

Left femur tissue samples from 3 Sham and 3 SOP rats were collected. Total RNA was extracted from these samples using TRIzol reagent (Ambion, USA). mRNA was reverse-transcribed into cDNA using a test kit (Yeasen). The cDNA was diluted with RNase/DNase-free reagents. The specific reaction system and primer sequences are listed in Supplementary Table S2. RT-qPCR analysis was then performed according to the amplification conditions shown in Supplementary Table S2. The data were analyzed using the 2^–ΔΔCt method, with GAPDH as the reference gene for normalization (34). Experimental technique was repeated three times. Statistical analysis was performed using GraphPad Prism 5 software, and intergroup differences were compared by t-test. A P-value < 0.05 was considered statistically significant.

2.14 Statistical analysis

All analyses were performed using R software (v4.2.2). Differences between groups were assessed using the Wilcoxon test. For RT-qPCR, WB, ELISA, and IHC, comparisons of 2^–ΔΔCt values were performed using the t-test (P < 0.05).

3 Results

3.1 Identification of 817 DEGs and 76 DE-URGs

The flowchart of the present study was shown in Figure 1 (Figure 1). A total of 817 DEGs were identified from high and low BMD samples in the training cohort, consisting of 384 upregulated DEGs and 433 downregulated DEGs (Figures 2A, B). Upon overlapping these 817 DEGs with 1,367 URGs, 76 DE-URGs were successfully screened (Figure 2C). Enrichment analysis of the DE-URGs revealed 289 GO terms, including “protein polyubiquitination,” “protein deubiquitination,” and “proteasome-mediated ubiquitin-dependent protein catabolic process” (Figure 2D). These processes are part of the cellular mechanisms that regulate protein levels, activity, and function through dynamic modifications involving the addition or removal of ubiquitin or other small protein groups. In KEGG analysis, two significant pathways were identified: such as ubiquitin-mediated proteolysis (Figure 2E). These processes are crucial for protein quality control and homeostasis, with the ubiquitin-mediated proteolysis pathway playing a vital role in degrading improperly folded proteins recognized in the endoplasmic reticulum (ER), as well as those in the cytosol and nucleus. The ER ensures that only properly folded proteins progress through the secretory pathway, while misfolded proteins are either refolded or marked for destruction by the ubiquitin-proteasome system.

Figure 1

Figure 1. Flow chart.

Figure 2

Figure 2. Identification of differentially expressed genes and ubiquitination-related genes in low bone mineral density samples. (A) Volcano map of differentially expressed genes in high and low bone mineral density samples. Divided by reference lines, the genes in the upper right corner are upregulated genes, the genes in the upper left corner are downregulated genes, and the rest are genes with no significant statistical difference. (B) Heat maps of differentially expressed genes in high and low bone mineral density samples. Yellow indicates highly expressed genes, and blue indicates lowly expressed genes. (C) Venn diagram showing the intersection of ubiquitination associated genes (URG) and DEGs. (D, E) GO (D) and KEGG (E) enrichment analysis of intersection difference URG.

3.2 Combining RPS27A, PSMD4, UBE2E1, UCHL3, and RAD23A into a model of high diagnostic value

The 76 DE-URGs were incorporated into a PPI network, consisting of 50 nodes and 137 edges (Figure 3A). The top 10 genes identified by the MCC and Degree algorithms from the CytoHubba plugin were intersected, yielding six candidate genes: RPS27A, UBB, PSMD4, UBE2E1, UCHL3, and RAD23A (Figures 3B, C). These genes were further refined using a LASSO model, which minimized at lambda = 0.11806, leaving five genes: RPS27A, PSMD4, UBE2E1, UCHL3, and RAD23A (Figures 3D, E). An SVM model was subsequently constructed. The confusion matrix indicated a high number of correct predictions, with the AUC of 0.945, and the PR curve yielded an AUC of 0.943 (Figures 3F–H). These results demonstrate that the model is highly effective in distinguishing between high and low BMD samples. Decision curve analysis (DCA) further confirmed that the model provided a positive net benefit, highlighting its potential value in practical applications (Figure 3I).

Figure 3

Figure 3. Construction of SVM model for candidate genes. (A) The PPI network map was constructed for the intersection difference URG genes. (B) The Venn diagram shows the intersection of MCC and Degree methods to obtain candidate genes. (C) PPI networks of candidate genes. (D, E) TMRGs related differential genes were screened by LASSO. (D) LASSO logic coefficient penalty diagram. The x-axis represents the logarithm of lambdas, and the y-axis represents the model error. (E) LASSO coefficient spectrum. The x-axis represents the logarithm of lambdas, and the y-axis represents the variable coefficients. (F–H). (F) The confusion matrix. The X-axis represents actual values and the Y-axis represents predicted values. The top-left corner of the figure denotes True Positive (TP), the bottom-left corner denotes False Negative (FN), the top-right corner denotes False Positive (FP), and the bottom-right corner denotes True Negative (TN). Values on the diagonal indicate the count/proportion of correct predictions, while off-diagonal elements represent misclassifications. (G) ROC curve. The X-axis is the false positive rate (FPR) and the Y-axis is the true positive rate (TPR). AUC refers to the area under the curve, with a larger value indicating higher accuracy. (H) PR curve. The Y-axis is Precision and the X-axis is Recall, which is obtained by adjusting the threshold for distinguishing between Positive and Negative classes. (I) DCA. The X-axis represents Threshold Probability, and the Y-axis represents the net benefit rate after subtracting harms from benefits.

3.3 Recognition of RPS27A and UBE2E1 as biomarkers

The expression levels of the five filtered genes were evaluated in both the training and validation cohorts. Notably, RPS27A and UBE2E1 exhibited consistent expression patterns across both datasets and showed significant differences (P < 0.05), with lower expression levels observed in samples with low BMD (Figures 4A, B). These findings led to the selection of RPS27A and UBE2E1 as biomarkers for this study, and they were subsequently used to construct a nomogram for predicting the risk of individual patients developing SOP and distinguishing the pathological characteristics of patient groups with varying BMD (Figure 4C). The calibration curve demonstrated an almost perfect alignment between the model’s estimated probabilities and the actual outcomes, with a concordance index (C-index) of 0.713 (Figure 4D). DCA indicated that the model’s net benefit decreased as the threshold probability increased; however, the model performed well across the range of threshold probabilities, except for very low threshold values (Figure 4E). The clinical impact curve highlighted the clinical utility of this predictive model (Figure 4F). By analyzing the expression of URGs in SOP patients, the nomogram model not only provided a deeper understanding of the disease’s molecular mechanisms but also showed promise for developing innovative therapeutic strategies.

Figure 4

Figure 4. Identification of biomarkers (A, B) Expression levels of biomarkers in training sets (A) and validation sets (B) *p < 0.05, **p < 0.01, ***p < 0.001 were considered statistically significant (C) Nomogram of the SOP diagnostic model on the basis of 2 biomarkers. Each gene in the figure corresponds to a line segment marked with scale marks; the length of the segment reflects the magnitude of the factor’s contribution to the outcome event. The sum of individual scores gives the Total Point, and a higher Total Point indicates a better model. (D) The calibration was used to evaluate the diagnostic efficacy of the SOP diagnostic nomogram. (E) DCA illustrating the net binefit assessing the outcome. (F) The clinical impact curve showed the clinical effectiveness of the nomogram model. The red curve (Number of high risk) denotes the number of individuals classified as positive (high-risk) by the simple model at various threshold probabilities; the blue curve (Number of high risk with outcome) denotes the number of true positives at each threshold probability.

3.4 Correlations between biomarkers and biological functions: insights into signaling pathways and health implications

A moderate positive correlation was observed between the two biomarkers, RPS27A and UBE2E1 (cor = 0.35, P = 0.026) (Figure 5A). GSEA revealed that both RPS27A and UBE2E1 were enriched in pathways related to “Oxidative phosphorylation,” “Proteasome,” “Ribosome,” and “Insulin signaling pathway” (Figures 5B, C). These pathways are involved in neurodegenerative diseases, essential cellular functions (such as energy production, protein synthesis, and degradation), and critical intracellular signaling mechanisms, suggesting that these biomarkers could have significant mechanistic implications for SOP patients across varying BMD levels.

Figure 5

Figure 5. Correlation and biological functions of RPS27A and UBE2E1 (A) Correlation analysis showed a moderate positive correlation between RPS27A and UBE2E1 (B, C) GSEA shows the related pathways and functions enriched by RPS27A (B) and UBE2E1 (CD. Subcellular localization of RPS27A and UBE2E1. Yellow indicates content related to RPS27A and UBE2E1, and green indicates subcellular localization.

Additionally, a subcellular-biomarker network constructed using GeneCards showed that RPS27A is localized to the cytosol and extracellular space (Figure 5D). This localization suggests that RPS27A may play a role in fundamental cytoplasmic functions, such as protein synthesis through its interaction with ribosomes, and may also be involved in extracellular processes, such as protein release during cell damage. UBE2E1, on the other hand, was found in both the cytosol and the nucleus, indicating its potential involvement in intracellular transport, nucleo-cytoplasmic exchange, and its possible direct influence on gene expression and the maintenance of genomic stability within the nucleus.

3.5 Biomarkers as representatives of immune microenvironment status in SOP patients

To explore the correlation between biomarkers and the immune microenvironment in SOP patients with different BMD levels, an immune infiltration analysis was performed. The assessment of immune cell infiltration across 25 immune cell types in patients with high and low BMD revealed significant differences in five cell types (P < 0.05) (Figures 6A, B). Macrophages, monocytes, and T regulatory 1 cells (Tr1) exhibited higher abundance in high BMD samples (P < 0.05), while T helper 17 cells (Th17) and T follicular helper cells (Tfh) were more prevalent in low BMD samples (P < 0.05), potentially contributing to the pathogenesis of OP. Correlation analysis demonstrated that RPS27A exhibited moderate correlations with all immune cells except Th17 (|cor| > 0.3), while UBE2E1 showed moderate correlations with Tfh, Th17, and Tr1 (|cor| > 0.3) (Figure 6C). These results suggest that RPS27A and UBE2E1 may regulate specific immune cell types in the microenvironment. In particular, the strongest positive correlation was observed between Tfh and UBE2E1 (cor = 0.43, P < 0.01), suggesting that UBE2E1 may influence autoimmune diseases and inflammatory responses. Additionally, monocytes and RPS27A showed the strongest positive correlation (cor = 0.38, P < 0.05), indicating that RPS27A might modulate or respond to immune processes driven by monocytes.

Figure 6

Figure 6. Immunoinfiltration analysis of RPS27A and UBE2E1. (A, B) Stack (A) and box (B) showed that there were five types of immune cells with significant differences between the disease and control groups. The X-axis represents different samples, and the Y-axis represents the proportion of immune cells. *p < 0.05 was considered statistically significant. (C) Correlation analysis of RPS27A and UBE2E1 with five kinds of immune cells. The X-axis represents biomarkers, and the Y-axis represents differentially expressed immune cells.

3.6 Speculating the potential molecular mechanisms of biomarkers

The association between biomarkers and SOP was further explored using the CTD, which revealed a strong correlation between UBE2E1 and SOP, with an inference score greater than 4. In contrast, RPS27A exhibited a moderate correlation with SOP, with an inference score over 2 (Figure 7A). Additionally, potential drugs targeting these biomarkers were identified using the DGIdb database, with RPS27A interacting with drugs such as EXALUREN, CYCLOHEXIMIDE, ATALUREN, MT-3724, and DORLIMOMAB ARITOX (Figure 7B).

Figure 7

Figure 7. Potential mechanism prediction and drug prediction analysis of RPS27A and UBE2E1 in SOP occurrence (A) The CTD database evaluates the relationship between RPS27A and SOP and the relationship between UBE2E1 and SOP. (B) Potential therapeutic agents for RPS27A and UBE2E1. (C) Network diagram of RPS27A and UBE2E1 interacting with transcription factors. Orange dots represent mRNAs, and purple dots represent transcription factors. (D) Volcanic maps of miRNA differential expression in disease and control groups. Divided by reference lines, the genes in the upper right corner are upregulated miRNA, the genes in the upper left corner are downregulated miRNA, and the rest are miRNA with no significant statistical difference. (E) Heat maps of miRNA differential expression in disease and control groups.Yellow indicates highly expressed miRNA, and blue indicates lowly expressed miRNA. (F) The Venn diagram indicated that the predicted miRNA of RPS27A and UBE2E1 were intersected with significantly differentially expressed mRNA. (G) The miRNA network diagram of RPS27A and UBE2E1 was adjusted jointly. Red dots represent mRNAs, and yellow dots represent miRNAs. (H) Construction of ceRNA network diagram for RPS27A and UBE2E1. Blue dots represent lncRNAs, yellow dots represent miRNAs, and red dots represent mRNAs.

To identify potential regulatory TFs for the biomarkers, a TF-mRNA network was constructed using the miRNet database, resulting in a network comprising 23 edges and 24 nodes. MAX was identified as a regulator of both genes (Figure 7C). Differential analysis of miRNA expression between SOP patients and control subjects revealed 44 DE-miRNAs, with 35 upregulated and 9 downregulated (Figures 7D, E). Using the multiMiR package, 195 miRNAs were predicted to target RPS27A and 51 miRNAs to target UBE2E1. By intersecting these miRNAs with those showing inverse expression to the mRNAs, a set of 12 miRNAs was identified, leading to the construction of a miRNA-mRNA regulatory network (Figure 7F). Among these, three miRNAs, including hsa-miR-106b-5p, were found to co-regulate both RPS27A and UBE2E1 (Figure 7G). Further investigation in the starBase database identified corresponding lncRNAs for these 12 miRNAs, resulting in the construction of a complex lncRNA-miRNA-mRNA network with 88 nodes and 279 edges, which included two mRNAs, 12 miRNAs, and 64 lncRNAs (Figure 7H). This network revealed multiple interactions, such as the regulation of UBE2E1 by lncRNAs like AC005538.2, XIST, NEAT1, MALAT1, and AC130462.1 via hsa-miR-185-5p.

3.7 Successful establishment of the SOP model and verification of biomarkers

Vaginal smear analysis revealed that rats in the Sham group remained in the estrus stage, while those in the SOP group predominantly exhibited white blood cells, indicating successful induction of senescence (Supplementary Figure 1A). No significant differences in the gross appearance of femoral tissues were observed between the Sham and SOP groups (Supplementary Figure 1B). Additionally, the body weight of SD rats increased over time in both groups (Supplementary Figure 1C). H&E staining showed that the bone tissue structure in the Sham group was dense and regular, with tightly arranged and continuous trabeculae. In contrast, the SOP group displayed sparse trabecular structure, increased gaps between trabeculae, and numerous irregular cavity-like structures, confirming successful model establishment (Figure 8A). ELISA results revealed that, compared to the Sham group, BGP and OPG levels were significantly reduced in the SOP group (P < 0.05), while TRACP5b levels were elevated (P < 0.01) (Figure 8B). IHC analysis indicated a decrease in the expression of RPS27A and UBE2E1 in the SOP group compared to the Sham group (Figures 8C, D). WB analysis confirmed that expression levels of both biomarkers were lower in the SOP group than in the Sham group (P < 0.05) (Figure 8E). Notably, RT-qPCR results were consistent with those of WB, showing significantly lower expression of the biomarkers in the SOP group compared to the Sham group (P < 0.01) (Figure 8F).

Figure 8

Figure 8. Validation of biomarkers in a rat model of SOP (A) HE staining of femur tissues in control group and disease group.10× Scale bars, 1mm,200× Scale bars, 50μm. (B) Serum concentrations of BGP, OPG and TRACP5b in control group and disease group were detected by ELISA Representative images from immunohistochemical staining of RPS27A (C) and UBE2E1 (D) in the Femur tissues. 10× Scale bars, 1mm,200× Scale bars, 50μm. (E) Detection of RPS27A and UBE2E1 in femur tissues from control group and disease group by western blot. Total GAPDH was used as a loading control. Data are mean ± s.d. (n=3),*p<0.05, by Student’s t-test. (F) RPS27A and UBE2E1 mRNA in femur tissue of rats in control group and disease group were measured by RT-qPCR. Data are mean ± s.d. (n=3). *p<0.05, ** p < 0.01, and *** p < 0.001 were considered statistically significant by Student’s t-test.

4 Discussion

SOP is a prevalent metabolic bone disease, particularly in postmenopausal elderly women (35). Age-related bone loss is primarily caused by an imbalance in bone remodeling (36). It has been reported that the ubiquitin-proteasome pathway plays a pivotal role in bone remodeling (15). However, the biomarkers and molecular mechanisms of this pathway in elderly patients with OP remain unclear. This study mined the GEO database to analyze DEGs and constructed a PPI network for DEGs related to ubiquitination. By combining the MCC algorithm and Degree algorithm in the Cytoscape plug-in with LASSO analysis, five key genes were identified: RPS27A, PSMD4, UBE2E1, UCHL3, and RAD23A. Subsequently, an SVM model was constructed based on these five genes to assess the diagnostic accuracy of the model in distinguishing between the disease and control groups. Further analysis revealed that RPS27A and UBE2E1 exhibited significant differences and consistent expression trends in both the training and validation sets. This multi-dimensional screening approach evaluates gene importance from both local and global perspectives, reducing the number of candidate biomarkers and improving the accuracy and reliability of the screening process (37–40). This study provides a theoretical foundation and data support for the prevention and treatment of OP in elderly women.

RPS27A (Ribosomal Protein S27a) and UBE2E1 (Ubiquitin Conjugating Enzyme E2E1) are two biomarkers involved in cellular function and various biological processes. The RPS27A gene encodes a component of the 40S ribosomal subunit, playing a pivotal role in protein synthesis (41). The encoded protein is also part of the ubiquitin-proteasome pathway, functioning as a fusion protein with an N-terminal ubiquitin tag (42) that facilitates the degradation of target proteins (42). Previous studies have shown that RPS27A expression is significantly down-regulated in postmenopausal women with OP (43), and recent research has also found down-regulation of RPS27A expression in elderly patients with OP (44), consistent with our findings. Compared with previous studies, this research incorporates two key influencing factors of aging and menopause, focusing on postmenopausal elderly women as the specific population. Through multiple-step screening, RPS27A is identified as a biomarker with diagnostic value for osteoporosis in this specific group. Moreover, the support vector machine disease prediction model constructed further validates this conclusion. Additionally, this research also reveals the association between RPS27A and the immune microenvironment, enriching the understanding of the role of this gene in the pathogenesis of osteoporosis.The UBE2E1 gene encodes a member of the ubiquitin ligase family, playing a critical role in the ubiquitination pathway by attaching ubiquitin tags to proteins, marking them for proteasomal degradation (45). Some studies have suggested that UBE2E1 may be involved in neurodegenerative diseases such as Alzheimer’s disease and Parkinson’s disease (46), and elevated levels of UBE2E1 may also be associated with inflammatory and autoimmune diseases (47). To date, there have been no literature reports indicating a correlation between it and osteoporosis. However, this study is the first to discover the correlation between the two and link it to the disease risk. Some researchers have found that this gene can regulate the classical WNT pathway (48), and the Wnt/β-catenin signaling pathway is considered a key pathway for osteoblast differentiation and has become an important target for the treatment of osteoporosis (49). UBE2E1 may be involved in the occurrence and development of senile osteoporosis through the classical Wnt pathway. This finding has not been reported in the literature yet and requires further verification. Overall, RPS27A and UBE2E1 are two newly identified key ubiquitination genes that hold promise not only as biomarkers for OP but also as potential therapeutic targets. Notably, this study developed a nomogram that integrates these two diagnostic markers with high AUC values and excellent calibration. This nomogram demonstrated outstanding accuracy and reliability in diagnosing OP and is expected to be clinically applied to facilitate early diagnosis of SOP.

GSEA was performed to further explore the potential functions and pathways associated with these two biomarkers in OP. The analysis revealed co-enrichment in seven pathways: “Alzheimer’s disease,” “Huntington’s disease,” “Oxidative phosphorylation,” “Parkinson’s disease,” “Proteasome,” “Ribosome,” and “Insulin signaling pathway.” The insulin signaling pathway is a critical, versatile pathway with multiple anabolic functions beyond glucose homeostasis (50). It regulates essential aspects of organismal health, such as growth, longevity, metabolism, and reproduction (50). Activation of the insulin signaling pathway can enhance osteoblast differentiation and promote bone formation (51, 52), while dysfunction in this pathway leads to bone formation disorders and OP (53). Furthermore, several studies have suggested a link between the development of OP in the elderly and ribosome-related genes and pathways (44, 54). The proteasome pathway recognizes proteins tagged with ubiquitin and degrades them through hydrolysis. Our study further confirms that the ubiquitin-proteasome pathway plays a significant role in OP, particularly in elderly patients. The ubiquitination of proteins generally occurs in the cytoplasm. The subcellular localization analysis of the biomarkers revealed that both RPS27A and UBE2E1 are expressed in the cytoplasm, further supporting the idea that key steps in regulating OP development may occur in the cytoplasm.

During the onset and progression of OP, the immune system plays a pivotal role in bone metabolism. Chronic low-grade inflammation is now recognized as a significant contributor to OP, particularly in postmenopausal women and the elderly (55). In the present study, immune infiltration analysis revealed significant differences in the infiltration levels of macrophages, monocytes, Tfh, Th17, and Tr1 cells. Th17 cells, a subset of T cells, promote inflammation by directly expressing RANKL, which facilitates the interaction between receptor activator of nuclear factor-κB ligand (RANKL)and receptor activator of nuclear factor-κB (RANK). These cells also secrete inflammatory cytokines such as IL-17, TNF-α, and IL-6, which stimulate osteoclasts (OCs) and enhance inflammatory infiltration. This process increases the expression of NF-κB, which further upregulates RANKL expression, promoting OC maturation and differentiation (56). However, the role of Th17 cells in osteoblasts remains controversial, with studies showing both supportive and inhibitory effects on bone metabolism (8). Increased IL-17 levels in postmenopausal OP patients correlate with a higher frequency of Th17 cells (57). Aging also leads to an increased Th17/Treg ratio, which may induce a low-activation state of the immune system, promote the production of inflammatory factors, and contribute to bone loss (58). Our study found that Th17 cells have a high infiltration abundance in elderly postmenopausal OP patients, and UBE2E1 is significantly negatively correlated with Th17, suggesting that UBE2E1 may influence the internal environment of OP patients through Th17 cells.

Tfh, a subset of CD4+ Th cells, play a critical role in B cell activation, proliferation, and differentiation (59). These processes predominantly occur in B cell follicles within secondary lymphoid organs and are essential for mounting effective antibody responses (59). Tfh are also involved in various conditions, including cancer, allergies, autoimmune diseases, and infections (59). Although Tfh have widespread effects, their association with OP has not been well documented. However, Tfh may be involved in bone metabolism disorders caused by immune-related diseases such as rheumatoid arthritis (60), systemic lupus erythematosus (61), and schistosomiasis (62). In our study, both RPS27A and UBE2E1 were significantly negatively correlated with Tfh, and their expression levels were high in OP patients, indicating that RPS27A may also serve as an immune factor in the pathogenesis of SOP. These findings offer new insights for future research, guiding researchers to explore the specific roles of these genes in immune responses and related diseases. Additionally, they may contribute to the development of novel therapeutic strategies for OP.

Through our research combined with previous literature reports, we discovered that UBE2E1 can work in synergy with the E3 ligase RNF34 to ubiquitinate ciliary protein MKS1. This process can precisely regulate the level of the core transcription factor β-catenin in the classical Wnt signaling pathway (48), and the abnormality of β-catenin can activate Th17 differentiation and inhibit Treg cells (15), thereby disrupting the bone immune balance and leading to increased bone resorption. Similarly, in the study of RPS27A, it was found that RPS27A can promote M2 polarization of macrophages. Although this study mainly focused on cancer rather than bone diseases, it revealed the regulatory role of RPS27A on macrophage polarization (63), and macrophage polarization plays a key role in bone remodeling (64, 65). Therefore, we speculate that as a ribosomal protein and ubiquitin precursor, RPS27A may participate in the pathogenesis of osteoporosis by influencing macrophage polarization through releasing certain signals. We believe that the ubiquitination system is the core regulator of immune cell activation and inflammatory factor signal transduction, and its participation in the chronic inflammatory response in the osteoporosis microenvironment may lead to the disruption of bone homeostasis and is an important link in the pathogenesis of osteoporosis.

This study predicted the upstream regulators of RPS27A and UBE2E1 through molecular regulatory network analysis, identifying hsa-miR-106b-5p, hsa-miR-17-5p, and hsa-miR-30c-5p as key regulators of these genes. Notably, hsa-miR-106b-5p has been shown to influence bone formation and OP through the regulation of TCF4 (66). It is hypothesized that RPS27A and UBE2E1 are regulated by the upstream miRNA hsa-miR-106b-5p, contributing to the promotion of OP. This insight could offer new perspectives and potential drug targets for OP research. Finally, we established the SOP model and verified the expression levels of RPS27A and UBE2E1 at both the protein and mRNA levels. It was found that the protein expressions of RPS27A and UBE2E1 in the SOP group were lower than those in the control group. The immunohistochemical results were also in line with expectations, and this was consistent with the expression at the transcriptional level. The above results indicate that RPS27A and UBE2E1 may be involved in the pathological process of SOP through the ubiquitin-mediated pathway. This study provides new insights into the role of the ubiquitin-mediated pathway and supports the potential of these two biomarkers as biomarkers for diagnosing and treating SOP.

In conclusion, by analyzing elderly postmenopausal OP samples and healthy samples from the GEO high-throughput database and applying bioinformatics and machine learning techniques, this study identified RPS27A and UBE2E1 as biomarkers in the ubiquitin-proteasome pathway involved in the development of OP in the elderly. The findings offer valuable insights into the molecular mechanisms of ubiquitination in SOP. However, these conclusions are based on health data analysis, and additional clinical sample data are required for further validation. Future research will continue to explore the role of these mechanisms and conduct experimental studies to confirm their scientific accuracy.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Ethics statement

Ethical approval was not required for the study involving humans in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements. The animal study was approved by Yunnan Bestime Biological Technology LLC. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

XC: Writing – original draft, Writing – review & editing, Data curation, Formal analysis, Investigation, Methodology. JL: Data curation, Investigation, Software, Writing – review & editing. YG: Data curation, Investigation, Writing – original draft. BJ: Formal analysis, Investigation, Writing – original draft. JZ: Investigation, Software, Writing – original draft. YC: Formal Analysis, Writing – original draft.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This research was supported by Hebei Natural Science Foundation (H2022206383).

Acknowledgments

Thanks to the Hebei Natural Science Foundation, Hebei Medical University.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1629276/full#supplementary-material

Abbreviations

SOP, Senile osteoporosis; BMD, Bone mineral density; PPI, Protein–protein interaction; MNC, Maximum Neighborhood Component; LASSO, Least Absolute Shrinkage and Selection Operator; URGs, Ubiquitination-related genes; SVM, Support Vector Machine; ELISA, Enzyme-linked immunosorbent assay; Tfh, T follicular helper cells; Th17, T helper 17 cells; TF, Transcription factor; OP, Osteoporosis; E1, Ubiquitin activating enzyme; E2, Ubiquitin conjugating enzyme; E3, Ubiquitin protein ligase; BMSCs, Bone marrow stromal cells; DEGs, Differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MCC, Maximal Clique Centrality; ROC, Receiver operating characteristic; AUC, Area under the curve; GSEA, Gene Set Enrichment Analysis; WB, Western blot; H&E, Hematoxylin–eosin; IHC, Immunohistochemical; OPG, Osteocalcin; TRACP5b, Tartrate-resistant acid phosphatase 5b; ER, Endoplasmic reticulum; Tr1, T regulatory 1 cells; RPS27A, Ribosomal protein S27a; UBE2E1, Ubiquitin conjugating enzyme E2E1; UBE2E3, Ubiquitin conjugating enzyme E2E3; OC, Osteoclasts; RANKL, Receptor Activator of Nuclear Factor-κB Ligand; RANK, Receptor Activator of Nuclear Factor-κB.

References

1. Amarnath SS, Kumar V, and Das SL. Classification of osteoporosis. Indian J orthopaedics. (2023) 57:49–54. doi: 10.1007/s43465-023-01058-3

PubMed Abstract | Crossref Full Text | Google Scholar

2. Aibar-Almazán A, Voltes-Martínez A, and Castellote-Caballero Y. Current status of the diagnosis and management of osteoporosis. Int J Mol Sci. (2022) 23:9465. doi: 10.3390/ijms23169465

PubMed Abstract | Crossref Full Text | Google Scholar

3. Wu J, Niu L, Yang K, Xu J, Zhang D, Ling J, et al. The role and mechanism of RNA-binding proteins in bone metabolism and osteoporosis. Ageing Res Rev. (2024) 96:102234. doi: 10.1016/j.arr.2024.102234

PubMed Abstract | Crossref Full Text | Google Scholar

4. Hu Y, Xu H, Ji W, Yang J, Li H, Li K, et al. Prevalence of frailty in senile osteoporosis: A systematic review and meta-analysis. Arch gerontology geriatrics. (2025) 130:105718. doi: 10.1016/j.archger.2024.105718

PubMed Abstract | Crossref Full Text | Google Scholar

5. Pizza IC, Bongiorno A, Pedullà M, Albano D, Sconfienza LM, and Messina C. DXA: new concepts and tools beyond bone mineral density. Semin musculoskeletal radiology. (2024) 28:528–38. doi: 10.1055/s-0044-1788579

PubMed Abstract | Crossref Full Text | Google Scholar

6. Zhang C, Yu H, Miao Y, and Wei B. Causal relationship between osteoporosis, bone mineral density, and osteonecrosis: a bidirectional two-sample Mendelian randomization study. J Trans Med. (2025) 23:226. doi: 10.1186/s12967-024-06030-9

PubMed Abstract | Crossref Full Text | Google Scholar

7. Tabib S, Alizadeh SD, Andishgar A, Pezeshki B, Keshavarzian O, and Tabrizi R. Diagnosis osteoporosis risk: using machine learning algorithms among fasa adults cohort study (FACS). Endocrinol Diabetes & Metab. (2025) 8:e70023. doi: 10.1002/edm2.70023

PubMed Abstract | Crossref Full Text | Google Scholar

8. Zhang W, Gao R, Rong X, Zhu S, Cui Y, Liu H, et al. Immunoporosis: Role of immune system in the pathophysiology of different types of osteoporosis. Front endocrinology. (2022) 13:965258. doi: 10.3389/fendo.2022.965258

PubMed Abstract | Crossref Full Text | Google Scholar

9. Liu J, Gao Z, and Liu X. Mitochondrial dysfunction and therapeutic perspectives in osteoporosis. Front endocrinology. (2024) 15:1325317. doi: 10.3389/fendo.2024.1325317

PubMed Abstract | Crossref Full Text | Google Scholar

10. Liu H, Zhao H, Lin H, Li Z, Xue H, Zhang Y, et al. Relationship of COL9A1 and SOX9 genes with genetic susceptibility of postmenopausal osteoporosis. Calcified Tissue Int. (2020) 106:248–55. doi: 10.1007/s00223-019-00629-7

PubMed Abstract | Crossref Full Text | Google Scholar

11. Popovic D, Vucic D, and Dikic I. Ubiquitination in disease pathogenesis and treatment. Nat Med. (2014) 20:1242–53. doi: 10.1038/nm.3739

PubMed Abstract | Crossref Full Text | Google Scholar

12. Hershko A and Ciechanover A. The ubiquitin system. Annu Rev Biochem. (1998) 67:425–79. doi: 10.1146/annurev.biochem.67.1.425

PubMed Abstract | Crossref Full Text | Google Scholar

13. Saritas-Yildirim B and Silva EM. The role of targeted protein degradation in early neural development. Genesis (New York NY: 2000). (2014) 52:287–99. doi: 10.1002/dvg.22771

PubMed Abstract | Crossref Full Text | Google Scholar

14. Varshavsky A. The ubiquitin system, an immense realm. Annu Rev Biochem. (2012) 81:167–76. doi: 10.1146/annurev-biochem-051910-094049

PubMed Abstract | Crossref Full Text | Google Scholar

15. Fan X, Zhang R, Xu G, Fan P, Luo W, Cai C, et al. Role of ubiquitination in the occurrence and development of osteoporosis (Review). Int J Mol Med. (2024) 54:68. doi: 10.3892/ijmm.2024.5392

PubMed Abstract | Crossref Full Text | Google Scholar

16. Li Y and Reverter D. Molecular mechanisms of DUBs regulation in signaling and disease. Int J Mol Sci. (2021) 22:986. doi: 10.3390/ijms22030986

PubMed Abstract | Crossref Full Text | Google Scholar

17. Liang C, Peng S, Li J, Lu J, Guan D, Jiang F, et al. Inhibition of osteoblastic Smurf1 promotes bone formation in mouse models of distinctive age-related osteoporosis. Nat Commun. (2018) 9:3428. doi: 10.1038/s41467-018-05974-z

PubMed Abstract | Crossref Full Text | Google Scholar

18. Liu Y, Cai G, Chen P, Jiang T, and Xia Z. UBE2E3 regulates cellular senescence and osteogenic differentiation of BMSCs during aging. PeerJ. (2021) 9:e12253. doi: 10.7717/peerj.12253

PubMed Abstract | Crossref Full Text | Google Scholar

19. Hariri H, Kose O, and Bezdjian A. USP53 regulates bone homeostasis by controlling rankl expression in osteoblasts and bone marrow adipocytes. J Bone mineral Res. (2023) 38:578–96. doi: 10.1002/jbmr.4778

PubMed Abstract | Crossref Full Text | Google Scholar

20. Lin YC, Zheng G, Liu HT, Wang P, Yuan WQ, Zhang YH, et al. USP7 promotes the osteoclast differentiation of CD14+ human peripheral blood monocytes in osteoporosis via HMGB1 deubiquitination. J orthopaedic translation. (2023) 40:80–91. doi: 10.1016/j.jot.2023.05.007

PubMed Abstract | Crossref Full Text | Google Scholar

21. Zhou Y, Gao Y, and Xu C. A novel approach for correction of crosstalk effects in pathway analysis and its application in osteoporosis research. Sci Rep. (2018) 8:668. doi: 10.1038/s41598-018-19196-2

PubMed Abstract | Crossref Full Text | Google Scholar

22. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | Crossref Full Text | Google Scholar

23. Yang L, Pan X, Zhang Y, Zhao D, Wang L, Yuan G, et al. Bioinformatics analysis to screen for genes related to myocardial infarction. Front Genet. (2022) 13:990888. doi: 10.3389/fgene.2022.990888

PubMed Abstract | Crossref Full Text | Google Scholar

24. Gustavsson EK and Zhang D. ggtranscript: an R package for the visualization and interpretation of transcript isoforms using ggplot2. Bioinformatics(Oxford England). (2022) 38:3844–6. doi: 10.1093/bioinformatics/btac409

PubMed Abstract | Crossref Full Text | Google Scholar

25. Gu Z and Hübschmann D. Make interactive complex heatmaps in R. Bioinf (Oxford England). (2022) 38:1460–2. doi: 10.1093/bioinformatics/btab806

PubMed Abstract | Crossref Full Text | Google Scholar

26. Chen H and Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinf. (2011) 12:35. doi: 10.1186/1471-2105-12-35

PubMed Abstract | Crossref Full Text | Google Scholar

27. Yu G, Wang LG, Han Y, and He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics: J Integr Biol. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | Crossref Full Text | Google Scholar

28. Liu P and Xu H. Potential molecular mechanisms of plantain in the treatment of gout and hyperuricemia based on network pharmacology. Evidence-Based complementary Altern medicine:eCAM. (2020) 2020:3023127. doi: 10.1155/2020/3023127

PubMed Abstract | Crossref Full Text | Google Scholar

29. Li Y, Lu F, and Yin Y. Applying logistic LASSO regression for the diagnosis of atypical Crohn's disease. Sci Rep. (2022) 12:11340. doi: 10.1038/s41598-022-15609-5

PubMed Abstract | Crossref Full Text | Google Scholar

30. Sachs MC. plotROC: A tool for plotting ROC curves. J Stat software. (2017) 79:1–19. doi: 10.18637/jss.v079.c02

PubMed Abstract | Crossref Full Text | Google Scholar

31. Wang L, Wang D, Yang L, Zeng X, Zhang Q, Liu G, et al. Cuproptosis related genes associated with Jab1 shapes tumor microenvironment and pharmacological profile in nasopharyngeal carcinoma. Front Immunol. (2022) 13:989286. doi: 10.3389/fimmu.2022.989286

PubMed Abstract | Crossref Full Text | Google Scholar

32. Miao YR, Zhang Q, Lei Q, Luo M, Xie GY, Wang H, et al. ImmuCellAI: A unique method for comprehensive T-cell subsets abundance prediction and its application in cancer immunotherapy. Advanced Sci (Weinheim Baden-Wurttemberg Germany). (2020) 7:1902880. doi: 10.1002/advs.201902880

PubMed Abstract | Crossref Full Text | Google Scholar

33. Ru Y, Kechris KJ, Tabakoff B, Hoffman P, Radcliffe RA, Bowler R, et al. The multiMiR R package and database: integration of microRNA-target interactions along with their disease and drug associations. Nucleic Acids Res. (2014) 42:e133. doi: 10.1093/nar/gku631

PubMed Abstract | Crossref Full Text | Google Scholar

34. Livak KJ and Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods (San Diego Calif). (2001) 25:402–8. doi: 10.1006/meth.2001.1262

PubMed Abstract | Crossref Full Text | Google Scholar

35. Hejazi K and Askari R. Effects of physical exercise on bone mineral density in older postmenopausal women: a systematic review and meta-analysis of randomized controlled trials. Arch osteoporosis. (2022) 17:102. doi: 10.1007/s11657-022-01140-7

PubMed Abstract | Crossref Full Text | Google Scholar

36. Rendina-Ruedy E and Rosen CJ. Lipids in the bone marrow: an evolving perspective. Cell Metab. (2020) 31:219–31. doi: 10.1016/j.cmet.2019.09.015

PubMed Abstract | Crossref Full Text | Google Scholar

37. Yoo TK, Kim SK, Kim DW, Choi JY, Lee WH, Oh E, et al. Osteoporosis risk prediction for bone mineral density assessment of postmenopausal women using machine learning. Yonsei Med J. (2013) 54:1321–30. doi: 10.3349/ymj.2013.54.6.1321

PubMed Abstract | Crossref Full Text | Google Scholar

38. Wang X, Pei Z, Hao T, Ariben J, Li S, He W, et al. Prognostic analysis and validation of diagnostic marker genes in patients with osteoporosis. Front Immunol. (2022) 13:987937. doi: 10.3389/fimmu.2022.987937

PubMed Abstract | Crossref Full Text | Google Scholar

39. Wu LL, Zhou JX, Jia YM, and Leng H. Screening and bioinformatics analysis of senile osteoporosis genes based on GEO database. Eur Rev Med Pharmacol Sci. (2023) 27:4857–64. doi: 10.26355/eurrev_202306_32602

PubMed Abstract | Crossref Full Text | Google Scholar

40. Deng YX, He WG, Cai HJ, Jiang JH, Yang YY, Dan YR, et al. Analysis and validation of hub genes in blood monocytes of postmenopausal osteoporosis patients. Front endocrinology. (2021) 12:815245. doi: 10.3389/fendo.2021.815245

PubMed Abstract | Crossref Full Text | Google Scholar

41. Eastham MJ, Pelava A, and Wells GR. RPS27a and RPL40, which are produced as ubiquitin fusion proteins, are not essential for p53 signalling. Biomolecules. (2023) 13:898. doi: 10.3390/biom13060898

PubMed Abstract | Crossref Full Text | Google Scholar

42. Luo J, Zhao H, and Chen L. Multifaceted functions of RPS27a: An unconventional ribosomal protein. J Cell Physiol. (2023) 238:485–97. doi: 10.1002/jcp.30941

PubMed Abstract | Crossref Full Text | Google Scholar

43. Yu T, You X, Zhou H, He W, Li Z, Li B, et al. MiR-16-5p regulates postmenopausal osteoporosis by directly targeting VEGFA. Aging. (2020) 12:9500–14. doi: 10.18632/aging.103223

PubMed Abstract | Crossref Full Text | Google Scholar

44. Wang C, Jiang X, Li Q, Zhang YZ, Tao JF, and Wu CA. Identification of core pathogenic genes and pathways in elderly osteoporosis based on bioinformatics analysis. Zhonghua yu fang yi xue za zhi [Chinese J Prev medicine]. (2023) 57:1040–6. doi: 10.3760/cma.j.cn112150-20230221-00140

PubMed Abstract | Crossref Full Text | Google Scholar

45. Wu X, Du Y, and Liang LJ. Structure-guided engineering enables E3 ligase-free and versatile protein ubiquitination via UBE2E1. Nat Commun. (2024) 15:1266. doi: 10.1038/s41467-024-45635-y

PubMed Abstract | Crossref Full Text | Google Scholar

46. Larabee CM, Georgescu C, Wren JD, and Plafker SM. Expression profiling of the ubiquitin conjugating enzyme UbcM2 in murine brain reveals modest age-dependent decreases in specific neurons. BMC Neurosci. (2015) 16:76. doi: 10.1186/s12868-015-0194-y

PubMed Abstract | Crossref Full Text | Google Scholar

47. Wang T, Zeng F, Li X, Wei Y, Wang D, Zhang W, et al. Identification of key genes and pathways associated with sex differences in rheumatoid arthritis based on bioinformatics analysis. Clin Rheumatol. (2023) 42:399–406. doi: 10.1007/s10067-022-06387-6

PubMed Abstract | Crossref Full Text | Google Scholar

48. Szymanska K, Boldt K, Logan CV, Adams M, Robinson PA, and Ueffing M. Regulation of canonical Wnt signalling by the ciliopathy protein MKS1 and the E2 ubiquitin-conjugating enzyme UBE2E1. Elife. (2022) 11:e57593. doi: 10.7554/eLife.57593

PubMed Abstract | Crossref Full Text | Google Scholar

49. Wang X, Qu Z, Zhao S, Luo L, and Yan L. Wnt/β-catenin signaling pathway: proteins' roles in osteoporosis and cancer diseases and the regulatory effects of natural compounds on osteoporosis. Mol Med. (2024) 30:193. doi: 10.1186/s10020-024-00957-x

PubMed Abstract | Crossref Full Text | Google Scholar

50. Paloma Álvarez-Rendón J, Manuel Murillo-Maldonado J, and Rafael Riesgo-Escovar J. The insulin signaling pathway a century after its discovery: Sexual dimorphism in insulin signaling. Gen Comp endocrinology. (2023) 330:114146. doi: 10.1016/j.ygcen.2022.114146

PubMed Abstract | Crossref Full Text | Google Scholar

51. Jiang F, Zong Y, Ma X, Jiang C, Shan H, Lin Y, et al. miR-26a attenuated bone-specific insulin resistance and bone quality in diabetic mice. Mol Ther Nucleic Acids. (2020) 20:459–67. doi: 10.1016/j.omtn.2020.03.010

PubMed Abstract | Crossref Full Text | Google Scholar

52. Ma H, Ma JX, Xue P, Gao Y, and Li YK. Osteoblast proliferation is enhanced upon the insulin receptor substrate 1 overexpression via PI3K signaling leading to down-regulation of NFκB and BAX pathway. Exp Clin Endocrinol diabetes: Off journal German Soc Endocrinol [and] German Diabetes Assoc. (2015) 123:126–31. doi: 10.1055/s-0034-1390422

PubMed Abstract | Crossref Full Text | Google Scholar

53. Wong SK and Mohamad NV. A review on the crosstalk between insulin and wnt/β-catenin signalling for bone health. Int J Mol Sci. (2023) 24:12441. doi: 10.3390/ijms241512441

PubMed Abstract | Crossref Full Text | Google Scholar

54. Su Y, Yu G, Li D, Lu Y, Ren C, Xu Y, et al. Identification of mitophagy-related biomarkers in human osteoporosis based on a machine learning model. Front Physiol. (2023) 14:1289976. doi: 10.3389/fphys.2023.1289976

PubMed Abstract | Crossref Full Text | Google Scholar

55. Sapra L and Srivastava RK. Immunotherapy in the management of inflammatory bone loss in osteoporosis. Adv Protein Chem Struct Biol. (2025) 144:461–91. doi: 10.1016/bs.apcsb.2024.10.013

PubMed Abstract | Crossref Full Text | Google Scholar

56. Wang X, Sun B, Wang Y, Gao P, Song J, Chang W, et al. Research progress of targeted therapy regulating Th17/Treg balance in bone immune diseases. Front Immunol. (2024) 15:1333993. doi: 10.3389/fimmu.2024.1333993

PubMed Abstract | Crossref Full Text | Google Scholar

57. Bhadricha H, Patel V, Singh AK, Savardekar L, Patil A, Surve S, et al. Increased frequency of Th17 cells and IL-17 levels are associated with low bone mineral density in postmenopausal women. Sci Rep. (2021) 11:16155. doi: 10.1038/s41598-021-95640-0

PubMed Abstract | Crossref Full Text | Google Scholar

58. Schmitt V, Rink L, and Uciechowski P. The Th17/Treg balance is disturbed during aging. Exp gerontology. (2013) 48:1379–86. doi: 10.1016/j.exger.2013.09.003

PubMed Abstract | Crossref Full Text | Google Scholar

59. Lei H, Hu J, Zhu J, Li R, Zhao Y, Zhao Y, et al. Global research prospects and trends in TFH cells and tumors: a bibliometric analysis. Front Oncol. (2025) 15:1443890. doi: 10.3389/fonc.2025.1443890

PubMed Abstract | Crossref Full Text | Google Scholar

60. He J and Chu Y. Intestinal butyrate-metabolizing species contribute to autoantibody production and bone erosion in rheumatoid arthritis. Sci Adv. (2022) 8:eabm1511. doi: 10.1126/sciadv.abm1511

PubMed Abstract | Crossref Full Text | Google Scholar

61. Chen F, Wu Y, Ren G, and Wen S. Impact of T helper cells on bone metabolism in systemic lupus erythematosus. Hum Immunol. (2023) 84:327–36. doi: 10.1016/j.humimm.2023.04.003

PubMed Abstract | Crossref Full Text | Google Scholar

62. Li W, Wei C, Xu L, Yu B, Chen Y, Lu D, et al. Schistosome infection promotes osteoclast-mediated bone loss. PloS pathogens. (2021) 17:e1009462. doi: 10.1371/journal.ppat.1009462

PubMed Abstract | Crossref Full Text | Google Scholar

63. Kuai X, Wei C, He X, Wang F, Wang C, and Ji J. The potential value of RPS27A in prognosis and immunotherapy: from pan-cancer analysis to hepatocellular carcinoma validation. ImmunoTargets Ther. (2024) 13:673–90. doi: 10.2147/itt.S493217

PubMed Abstract | Crossref Full Text | Google Scholar

64. Wang S, Xiao L, Prasadam I, Crawford R, Zhou Y, and Xiao Y. Inflammatory macrophages interrupt osteocyte maturation and mineralization via regulating the Notch signaling pathway. Mol Med. (2022) 28:102. doi: 10.1186/s10020-022-00530-4

PubMed Abstract | Crossref Full Text | Google Scholar

65. Li M, Li D, Jiang Y, He P, Li Y, Wu Y, et al. The genetic background determines material-induced bone formation through the macrophage-osteoclast axis. Biomaterials. (2023) 302:122356. doi: 10.1016/j.biomaterials.2023.122356

PubMed Abstract | Crossref Full Text | Google Scholar

66. Chen H, Li S, Wang J, Ma Y, and Yin H. Screening of key biomarkers in osteoporosis: Evidence from bioinformatic analysis. Int J rheumatic diseases. (2023) 26:69–79. doi: 10.1111/1756-185x.14450

PubMed Abstract | Crossref Full Text | Google Scholar