Screening of Biomarkers Associated with Osteoarthritis Aging Genes and Immune Correlation Studies

Introduction

Osteoarthritis is a common degenerative joint disease of the middle-aged and elderly, with about 3 million new cases diagnosed each year.¹ The disease is a progressive, chronic disease that is characterized by damage to the cartilage of the joints and involves the entire joint tissue, manifesting itself in joint pain, stiffness, hypertrophy, and restricted movement, and is one of the main causes of disability in the middle-aged and elderly, and places a severe psychological, financial and physical burden on patients.^2,3

The pathological features of OA include degeneration of articular cartilage, narrowing of the joint space, inflammation of the synovial membrane, formation of bone fragments and alteration of mechanical properties.^4–6 Osteoarthritis is a comprehensive joint disease encompassing all joint tissues with risk factors including joint injury, poor joint development, genetics, age, obesity and gender.^7–9 Conventional diagnostic methods, including the combination of clinical presentation, physical examination, and imaging techniques like X-ray, are inadequate for the early detection of osteoarthritis.^10–12 Because of the difficulty of early diagnosis, most patients have limited treatment options, leading to poor prognoses. Screening of early biomarkers using bioinformatics has attracted the attention of researchers, developing early biomarkers of OA for disease diagnosis and providing therapeutic targets is of great significance.^13–15

Cell aging is one of the signs of aging and a kind of cell fate.¹⁶ Cellular senescence is characterized by cell cycle arrest and releasing harmful pro-inflammatory molecules into the pericellular microenvironment.¹⁷ As age is the most common risk factor for OA, aging is related to the onset and progression of OA.^18,19 In clinical practice, the efficacy of treating osteoarthritis (OA) is being evaluated through the use of senolytics. Senolytics are a category of medications that specifically trigger the death of aging cells in order to eliminate them. Examples of senolytics include p53 binding protein Mdm-2 inhibitors and B-cell lymphoma 2 (BCL-2) family inhibitors.²⁰ However, the role of cell senescence-related genes in the pathogenesis and progression of OA needs further study.

Therefore, based on the cellular senescence-related genes obtained from the CellAge database, this study identified OA cellular senescence-related biomarkers through a comprehensive bioinformatics approach, providing a theoretical basis for further understanding of disease pathogenesis and new potential targets for clinical diagnosis and treatment.

Materials and Methods

Gene and Dataset Collection

Five OA-related datasets, namely, GSE55235, GSE55457, GSE82107, GSE12021, and GSE51588, were retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) and enrolled in our study (Table 1). Datasets included in this study have to meet the following criteria: The organism of samples is Homo sapiens; The sample type is synovial tissue; dataset including OA and normal (negative control) samples; The dataset contains more than 3 samples per group. A total of 279 cellular senescence-related genes (CSRGs) were excavated from the CellAge database (https://genomics.senescence.info/cells/, accessed on 26 June 2022) (https://genomics.senescence.info/cells/, accessed on 26 June 2022).

Table 1 Sample Size and Information of the OA-Related Datasets

Preprocessing of Transcriptome Sequencing Data

Four OA-related datasets, namely GSE55235, GSE55457, GSE82107, and GSE12021, were merged by eliminating the batch effect through the combat function in R package “sva” (version 1.0.12).^21,22

Excavating the Candidate Genes for Cellular Senescence-Associated Biomarkers in OA

We firstly determined the OA-related genes by WGCNA to develop a co-expression network.²³ After samples clustering processing by the “goodSamplesGenes” function, outliers were recognized and then eliminated. The soft-thresholding power was computed to make the network contented the scale-free distribution. Then, modules (minimum genes number = 100) were developed by the dynamic tree-cutting method. We then mined the “limma” package for differentially expressed genes (DEGs) between OA and control samples. (version 3.4.3) with |log₂ (FC)| > 0.5 and p-value < 0.05.²⁴ Adopting the Venn diagram, we obtained the intersection of DEGs, OA-related genes, and CSRGs, which applied as candidate genes to screen out cellular senescence-associated biomarkers.

Functional Analysis of Genes

The “clusterProfiler” R package (version 4.0.2)²⁵ was used for the enrichment analysis of Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). Molecular function (MF), biological process (BP), and cellular component (CC) were the three main groups into which GO was divided. The criteria for significance was an adjusted p-value less than 0.05.

Certification of Cellular Senescence-Associated Biomarkers in OA

We implemented the Least absolute shrinkage and selection operator (LASSO)²⁵ logistic regression by the “glmnet” package (version 4.0.2) and support vector machines recursive feature elimination (SVM-RFE)²⁵ analysis by the “e1071” package (version 1.7–4). The recursive feature elimination (RFE) method was deployed to rank each gene based on its importance to procure feature genes. Finally, the overlapping genes were defined as cellular senescence-associated biomarkers.

Establishment of Diagnostic Model for OA

Based on the biomarkers obtained by machine learning, the “lrm” function in R was utilized to establish a diagnostic model and to visualize this model by drawing a nomogram. The corresponding calibration curve and DCA curve were plotted to assess the accuracy of the diagnostic model and patient benefit.

Single-Gene Gene Set Enrichment Analysis (GSEA)

The correlation coefficient between the expression of the target gene and all genes was calculated and used as a ranking criterion for GSEA enrichment analysis. The GO, KEGG, and reactome gene sets were extracted to act as the reference gene set. GSEA was executed using the R software “clusterProfiler” (version 4.0.2) package. Significance thresholds were adjusted p-value < 0.05.

Relevance Analysis of Cellular Senescence-Associated Biomarkers and Immune Cell Infiltration

The infiltration levels of 22 type immune cells in control and OA samples were profiled by the CIBERSORT algorithm.²⁶ A Wilcoxon test was used to evaluate the differences in the infiltration levels of various immune cells between the control and OA samples. The results were then shown using a violin plot. Using the Pearson approach, we evaluated the correlations between the various immune cells and the biomarkers.

Construction of lncRNA-miRNA-mRNA Network

The miRDB and StarBase database were deployed to predict the target miRNAs of biomarkers. Then the StarBase database and miNet database were employed to predict lncRNAs targeting the miRNAs. The intracellular localization of predicted lncRNAs was analyzed by the lncATLAS database (https://lncatlas.crg.eu/) to obtain lncRNAs localized in the cytoplasm. The final lncRNA-miRNA-mRNA network was constructed using Cytoscape (version 3.8.2).²⁷

Targeted Drug Prediction

With the CellMiner database, drugs targeting biomarkers were predicted by first selecting drugs that have been approved by the FDA or have been used in clinical trials and second retaining drug-gene relationship pairs with absolute values of correlation coefficients greater than 0.4 and p values less than 0.05. The final gene-drug network was mapped by Cytoscape software (version 3.8.2).

In vivo Validation of Key Genes

This study included 11 osteoarthritis patients and 11 healthy controls. The study included participants over 18 with radiologically verified knee osteoarthritis who satisfied American College of Rheumatology criteria. Orthopaedic surgeons diagnosed osteoarthritis. The study excluded patients with severe physical or behavioral problems, rheumatoid joints, or other joint disorders, and those who declined to participate. The healthy control tissues were taken from synovial joints of trauma patients. Trizol reagent (AG21101, China JSC) was used to extract total RNA from synovial tissues of 11 osteoarthritis patients and 11 healthy controls. A commercial reverse transcription kit (AG1706, China JSC) converted isolated RNA into cDNA. Real-time quantitative analysis was performed on the cDNA samples using the ABI7500 Real Time PCR System (7500 software v2.0.4, ABI, USA). The RT-PCR analysis used the SYBR premixed Ex Taq kit (AG11718, AG, China). The reference gene was 3-phosphoglyceraldehyde dehydrogenase. Table 2 lists the primers used to calculate relative mRNA levels using the 2(-Delta CT) technique. The Ethics Committee at Liaocheng Hospital, Shandong First Medical University approved the study (SWYX: NO.2023-005), and all participants gave informed consent.

Table 2 Primers for RT-PCR

Cell Culture and Transfection

The RPKFCL, derived from healthy synovial tissue of experimental rats, was obtained from Saibaikang (Shanghai, China). The cells were cultured in Fibroblast medium supplemented with 10% fetal bovine serum (FBS) from Gibco and 1% penicillin. The culture was maintained in an incubator at 37 degrees Celsius with 5% carbon dioxide. The whole PIM1 gene was put into the Ubi-MCS-3FLAG-SV40-IRES-puromycin Lentiviral Expression Vector and used to introduce it into RPKFCL by transfection. Perform the transfection of the lentiviral vector containing overexpressed PIM1 or overexpressed NC, together with two additional vectors (psPAX2 and pMD2.G), into the HEK293T cell line. Following a cultivation period of 48 hours, the liquid portion containing the desired substances was gathered and condensed for the purpose of introducing it into the RPKFCL cell line. Afterwards, puromycin was introduced into the RPKFCL culture system to acquire the targeted cells, and the effectiveness of PIM1 transfection was confirmed using q-PCR analysis.

Cell Proliferation Assay

A total of 5000 cells were distributed into each well of a 96-well plate and incubated at a constant temperature of 37 °C in a carbon dioxide-enriched environment of 5%. The CCK8 solution (Beyotime, China) was introduced and left to incubate for 40 minutes at the time of inoculation, as well as at 24 and 48 hours. Next, the measurement of absorbance at a wavelength of 450nm was conducted.²⁸

Cell Viability Analysis

The vitality of the RPKFCL cell line was evaluated using a live/dead cell kit (YEASEN, China) following the instructions provided by the manufacturer. The cells were introduced into six-well plates (Sangong, China) at a density of 1 x 10^⁵. The cells were incubated in an incubator at 37 degrees Celsius overnight, and then hydrogen peroxide (H2O2) was added to the culture system to a concentration of 100 nM, and then continue cultured for 2 hours to emulate the oxidative damage environment associated with osteoarthritis.^29,30 Afterward, the cells were stained with calcein and PI for 20 minutes at room temperature, and fluorescence images were captured using a fluorescence microscope manufactured by OLYMPUS in Japan.

Cell Senescence Analysis

The vitality of the RPKFCL cell line was evaluated using a live/dead cell kit (YEASEN, China) following the instructions provided by the manufacturer. Consistent with the steps mentioned above, RPKFCL were treated with hydrogen peroxide to mimic the oxidative damage environment of osteoarthritis.^29,30 Subsequently, the cells were stained with calcein and PI for 20 minutes at room temperature, and fluorescence images were captured using an fluorescence microscope (OLYMPUS, Japan).

Cell Cloning Experiment

The clonogenic ability of the RPKFCL cells was assessed by the execution of a colony formation experiment. A total of 1500 cells were distributed evenly over six-well plates, and the culture media was replaced at regular intervals of 3 days. Following a two-week period, the cells were subjected to fixation using a 4% paraformaldehyde solution obtained from Biosharp (China). Then, staining was performed using crystal violet obtained from Solarbio (China). The colonies were then subjected to three rounds of washing and then photographed.

Cell Migration Assay

The migratory ability of the cells was tested using a transwell chamber (Corning, USA). RPKFCL cells were introduced at a concentration of 5×104 cells per well in the top well, which contained medium without serum. The lower well was filled with 500 μL of complete medium. Following a 48-hour incubation period at 37°C in a temperature-controlled incubator, the non-migrating cells in the upper well were delicately eliminated. The cells in the lower chamber were treated with paraformaldehyde for 10 minutes to immobilize them, followed by staining with a 0.1% crystal violet solution for 20 minutes. Subsequently, the resulting images were examined and recorded using a microscope.

Statistical Analysis

The values are reported as the means ± standard deviation (SD) of a minimum of three separate experiments. The data were subjected to analysis using SPSS 25.0 software (IBM Corp. Released 2017. IBM SPSS Statistics for Windows, Version 25.0. Armonk, NY: IBM Corp), using an unpaired Student’s t-test to assess differences between the two groups. The receiver operating characteristic (ROC) curves were generated using the “pROC” software (version 1.17.0.1).²⁴ The research was conducted using the R programming language (4.0.2), and the data from separate groups were compared using the Wilcoxon test. When explicit instructions were not provided, a p-value less than 0.05 was considered to show statistical significance.

Results

Candidate Genes for Cellular Senescence-Associated Biomarkers in OA

First, we merged four OA-related datasets (GSE55235, GSE55457, GSE82107, and GSE12021) by removing the batch effect with the “sva” package. The expression box plot and PCA result of pre- and post-standardized were exhibited in Supplementary Figure 1A–D. Combined with Q-Q plot (Supplementary Figure 1), these results indicated that cross-platform normalization successfully removed the batch effect. The merged dataset contained 36 control samples and 40 OA samples. WGCNA was then implemented using the combined dataset. First, sample cluster analysis removed no outliers (Figure 1A). Five was selected as the ideal soft threshold (R2=0.85) to maximise scale-free gene interactions (Figure 1B). Next, 10 modules were created using a gene clustering tree and dynamic tree cutting technique (Figure 1C and D). Based on absolute correlation coefficients larger than 0.5, we chose the MEturquoise, MEbrown, and MEyellow modules as important modules for disease OA (Figure 1E). Thus, important module genes were OA-related (Supplementary Figure 2).

Figure 1 Candidate genes for cellular senescence-associated biomarkers in Osteoarthritis. Clustering diagram of merged dataset samples (A). Soft threshold analysis showed that gene associations were most consistent with the scale-free distribution when β = 5 (B). Clustering dendrogram of genes, dissimilarity based on the topological overlap, and assigned merged module colors and original module colors. (C and D). Heatmap of the relevance between gene modules and osteoarthritis (E). The training set on the differentially expressed genes screened under the conditions of | log2FC |>0.5 and p value<0.05. The red dot is the gene up-regulated, and the green dot is the gene down-regulated (F). The testing set on the differentially expressed genes screened under the conditions of | log2FC |>0.5 and p value<0.05. The red dot is the gene up-regulated, and the green dot is the gene down-regulated (G). Test set and training set (up-regulated gene and down-regulated gene) as well as OA-related gene and CSRG are displayed by overlapping Venn diagram. Four up-regulated genes and 12 down-regulated genes were screened (H and I). GO (J) and KEGG (K) enrichment analyses of 16 genes.

Next, the merged dataset was divided into a training set and a validation set according to a 7:3 ratio. The training set contained 29 OA samples and 25 control samples. The validation set included 11 OA samples and 11 control samples. To mine the variably expressed OA-related genes, we then authenticated the differentially expressed genes (DEGs) between OA and control samples in the training set and testing set respectively. According to |log₂FC| > 0.5 and p value < 0.05, a grand total of 1792 DEGs (OA vs control), including 860 up-regulated and 932 down-regulated genes in the OA samples, were identified in the training set (Figure 1F). Meanwhile, 1528 DEGs (OA vs control), including 805 up-regulated and 723 down-regulated genes in the synovial tissue of OA samples, were uncovered in the testing set (Figure 1G). In the following, we overlapped the DEGs, OA-related genes, and CSRGs. Sixteen DE-CSRGs associated with OA, including four intersecting up-regulated genes (Figure 1H) and twelve intersecting down-regulated genes (Figure 1I), namely EPHA3, TLR3, ZMAT3, TNFSF13, BCL6, MCL1, MYC, SLC16A7, VEGFA, ZFP36, PIM1, CEBPB, CDKN1A, BAG3, CYR61, and SNAI1, were identified as the candidate genes for cellular senescence-associated biomarkers in the synovial tissue of OA. Functional enrichment analysis was used to investigate these 16 genes. 770 GO items (714 BP, 10 CC, and 46 MF) and 33 KEGG pathways (Figure 1J) were generated. The bar chart showed the top 10 GO items per categorization (Figure 1K). We noted that these genes involved many cell differentiation-related and immune-related biological processes. For KEGG pathways, these DE-CSRGs from the synovial tissue associated with OA were linked to the “JAK-STAT signaling pathway”, “PI3K-Akt signaling pathway”, “p53 signaling pathway”, “ErbB signaling pathway”, “Rheumatoid arthritis”, “HIF-1 signaling pathway”, “Cell cycle”, “FoxO signaling pathway”, and “Cellular senescence”.

Cellular Senescence-Associated Biomarkers in OA

To further recognize cellular senescence-associated biomarkers in OA, we immediately executed machine learning analysis based on 16 candidate genes in the training set. As shown in Figure 2A and B, five genes (BCL6, MCL1, SLC16A7, PIM1, and EPHA3) were picked out by LASSO logistic regression. Meanwhile, nine genes (MCL1, EPHA3, PIM1, ZFP36, BAG3, SLC16A7, MYC, CEBPB, and BCL6) were selected by SVM-RFE model (Figure 2C and D). Hence, five overlapping genes (BCL6, MCL1, SLC16A7, PIM1, and EPHA3) were obtained by crossing the genes filtered by two machine learning methods. Next, we analyzed the diagnostic ability of the 5 genes as a whole in the training set, test set, and external validation set GSE51588. The ROC curves showed that the AUC values of the 5 genes as a whole were 0.975, 1, and 1 in the training set, test set, and external validation set, respectively, showing a powerful ability to distinguish OA samples from normal samples (Figure 2E–G). ROC curves of each gene were graphed for the training and testing datasets to evaluate the five genes’ diagnostic usefulness. The area under the curve (AUC) values of each gene surpassed 0.8, indicating that their expression could identify osteoarthritis (OA) samples from normal ones (Figure 2H and I). To further confirm the credibility of the above results, we created the corresponding ROC curves in the external validation set GSE51588. The AUC values of the ROC curves were all greater than 0.6, manifesting the tolerable ability to distinguish OA samples from normal samples (Figure 2J). Therefore, we defined these five genes as cellular senescence-associated biomarkers in OA.

Figure 2 Optimization algorithm for screening important Cellular senescence-associated biomarkers in OA. The 16 candidate genes previously screened were calculated using LASSO logistic regression analysis (A and B). SVM-RFE machine learning algorithm to select genes corresponding to the smallest classifier error as candidate key genes (C and D). ROC curves of 5 genes as a whole in the training set (E), test set (F) and external validation set (GSE51588) (G). The ROC curves of the single gene in the training set and testing set (H and I). The ROC curves of the single gene in the external validation set (J).

The OA Diagnostic Model Consisting of Cellular Senescence-Associated Biomarkers

Based on the five identified cellular senescence-associated biomarkers of OA, a nomogram for predicting the probability of disease, ie, a diagnostic model, was created (Figure 3A). The expression of each biomarker corresponds to a score, and the total score of the five biomarkers is used to estimate the risk of the sample developing OA. The higher the score, the greater the risk of developing the disease. The calibration curve revealed that the difference between the predicted and observed probabilities of disease was quite small, indicating that the nomogram performed satisfactorily (Figure 3B). In addition, the DCA curves also indicated that the nomogram model consisting of five genes had a higher benefit rate than individual genes (Figure 3C).

Figure 3 Diagnostic model for OA. Based on the five biomarkers obtained from the above analysis, a line graph of the biomarkers was constructed (A). Calibration curve of 5 genes (B). The decision curves show the predictive power of the line graph prediction model, with different coloured curves representing different genes and black curves representing the nomogram model (C).

Single-Gene GSEA

To explore the molecular mechanisms of cellular senescence-associated biomarkers in OA, we performed the single-gene GSEA based on GO, KEGG, and Reactome gene sets. The top 10 GO items, KEGG pathways, and Reactome pathways enriched by each gene were displayed in Figure 4 and Supplementary Figure 3. All five genes were related to “Oxidative phosphorylation”. Meanwhile, all five biomarkers were associated with some immune-related pathways, including Toll-Like Receptor Cascade, IL-17 signaling pathway, TNF signaling pathway, and NF-kappa B signaling pathway.

Figure 4 Single-gene GSEA based on KEGG, and reactome gene sets. KEGG results for BCL6 (A). Reactome gene results for BCL6 (B). KEGG results for EPHA3 (C). Reactome gene results for EPHA3 (D). KEGG results for MCL1 (E). Reactome gene results for MCL1 (F). KEGG results for PIM1 (G). Reactome gene results for PIM1 (H). KEGG results for SLC16A7 (I). Reactome gene results for SLC16A7 (J).

The Relevance of Cellular Senescence-Associated Biomarkers to Immune Cells in OA

Since all five cellular senescence-associated biomarkers were linked to immune-related pathways, we then probed the relationship between cellular senescence-associated biomarkers and immune cells in OA. We first compared the infiltration level of 22 immune cell types between OA samples and control samples through the CIBERSORT method. The proportion of each immune cell in each sample was shown in the bar stack (Figure 5A). The violin plot manifested that the fraction of T cells CD4 memory resting, NK cells activated, Dendritic cells activated, Mast cells activated, and Neutrophils in OA samples were significantly lower than control samples, and the fraction of B cells memory, Plasma cells, and Mast cells resting in OA samples were significantly higher than control samples (Figure 5B). The Pearson correlation coefficient was then used to determine the relationship between differential immune cells and indicators of cellular senescence. As illustrated in Figure 5C–H, on the basis of |correlation coefficient|>0.3 and p value<0.05, all five biomarkers were notably correlated with T cells CD4 memory resting and Mast cells resting. BCL6 and MCL1 were prominently correlated with NK cells activated. BCL6 was significantly correlated with dendritic cells activated. BCL6, MCL1, SLC16A7, and PIM1 were significantly correlated with Mast cells activated.

Figure 5 The relevance of cellular senescence-associated biomarkers to immune cells in OA. The composition of 22 kinds of immune cells in each sample showed in a histogram (A). The violin plot manifested that significantly different infiltrates of immune cells in OA (B). Analysis of differential immune cell and model gene correlations (C). Relationship between immune cells and 5 genes (BCL6, MCL1, SLC16A7, and PIM1) (D–H). *p<0.05; **p<0.01; ***p<0.001; ****p<0.0001. Abbreviation: ns, not significant.

OA lncRNA-miRNA-mRNA and Gene-Drug Networks Based on Cellular Senescence-Associated Biomarkers

To further probe the upstream regulatory mechanisms of biomarkers, we predicted 15 shared miRNAs targeting five biomarkers from the miRDB and StarBase database (Supplementary Figure 4A). Next, we predicted 323 shared lncRNAs targeting the 15 miRNAs from the Starbase and miNET database (Supplementary Figure 4B). The intracellular localization of 323 lncRNAs was identified through the lncATLAS database, of which 60 lncRNAs located in the cytoplasm were included in the subsequent lncRNA-miRNA-mRNA network construction (Supplementary Figure 5). Finally, a lncRNA-miRNA-mRNA network with 80 nodes (5 mRNAs, 15 miRNAs, and 60 lncRNAs) and 168 edges was generated using Cytoscape (Figure 6A). The relationship pairs in the lncRNA-miRNA-mRNA network were also exhibited in a Sankey diagram (Figure 6B).

Figure 6 The regulating mechanisms of 5 biomarkers. The lncRNA-miRNA-mRNA network was generated using Cytoscape (A). The relationship network between 5 genes and targeted drugs extracted from the CellMiner database (B). Gene-drug network containing 5 nodes and 26 edges (C).

To initially explore potential drugs targeting the five cellular senescence -associated biomarkers, we predicted 12 drugs targeting PIM1, 7 drugs targeting BCL6, 3 drugs targeting MCL1, 2 drugs targeting EPHA3, and 2 drugs targeting SLC16A7 through the CellMiner database. Ultimately, we constructed a gene-drug network containing 5 nodes and 26 edges (Figure 6C). In this network, vorinostat targeted PIM1.

Verification the Expression of Cellular Senescence-Associated Biomarkers in External Dataset

As illustrated in Figure 7A and B, BCL6, MCL1, SLC16A7, and PIM1were down-regulated in OA samples compared to control samples, and EPHA3 was up-regulated in OA samples compared to control samples in the training set and testing set. We further confirmed the expression in the external dataset (GSE51588). Consistent with the results of the training and testing sets, the expression of SLC16A7 and PIM1 were notably decreased in OA samples and EPHA3 was significantly up-regulated in OA samples compared to control samples (Figure 7C). Probably due to sample heterogeneity, the expression differences of BCL6 and MCL1 between OA samples and control samples were not significant, but there was a trend of decreasing expression in OA samples (Figure 7C).

Figure 7 Gene expression in the OA group compared with the standard group. In the training set (A) and test set (B), compared with the normal group, the expression of five genes in the OA group. Verification of critical gene expression in the external validation set (C). *p<0.05; **p<0.01; ***p<0.001; ****p<0.0001. Abbreviation: ns, not significant.

In vivo Verification of Key Genes

In addition, we also provided the RT-PCR results of the five hub genes (Figure 8). Compared to the healthy group, EPHA3 was expressed in the osteoarthritis group at higher levels, while BCL6, MCL1 SLC16A7, and PIM1 were expressed at lower. This result was consistent with the bioinformatics analysis.

Figure 8 In vivo verification of key genes. The PCR experiments were used to analyse the relative expression of key genes in osteoarthritic tissues. Compared to the healthy group, EPHA3 was expressed in the osteoarthritis group at higher levels, while BCL6, MCL1 SLC16A7, and PIM1 were expressed at lower. ns, not significant; ***p<0.001.

In vitro Validation of Key Genes

The relative expression of PIM1 in the oe-PIM1 group was shown to be significantly higher than that in the oe-NC group by q-PCR analysis (Figure 9A). The CCK8 assay indicated that PIM1 enhances the proliferative capacity of RPKFCL cell line (Figure 9B). The Live/Dead Cell kit suggested that PIM1 enhances the Viability of the RPKFCL cell line (Figure 9C). The cellular senescence analysis revealed that PIM1 enhances the proliferative capacity of RPKFCL (Figure 9D). The clone formation experiments demonstrated that PIM1 improves the proliferation of single cells of RPKFCL cell line (Figure 9E). Transwell assay displayed that PIM1 enhances the migration ability of RPKFCL cell line (Figure 9F).

Figure 9 In vitro validation of PIM1. The relative expression of PIM1 was significantly higher in the oe-PIM1 group (A). The CCK8 assay (B). The Cell viability assay (C). The cellular senescence analysis (D). The clone formation experiment (E). Transwell assay (F). *p<0.05; **p<0.01; ***p<0.001.

Discussion

With rapid advances in sequencing and computer technology, bioinformatics analysis has improved our understanding of the underlying mechanisms associated with complex diseases. Approximately 60% elderly and 10–15% young are afflicted by osteoarthritis disease, and aging is closely associated with the onset and progression of OA.³¹ Therefore, this study identifies cellular senescence-related biomarkers of OA and preliminarily explores their mechanisms of action and targeting agents.

First, we identified 16 differentially expressed cellular senescence-associated genes associated with OA by differential expression analysis and WGCNA. Functional enrichment results showed that these genes are associated with many immune-related biological processes, suggesting that the immune response may play an important role in OA.³² In addition, these genes are also enriched in the JAK-STAT signaling pathway, PI3K-Akt signaling pathway, HIF-1 signaling pathway, cell cycle, and cellular senescence KEGG pathways.³³ Studies have shown that the PI3K/AKT pathway is associated with the development of OA and changes in the extracellular matrix and that PI3K/AKT/mTOR signaling plays a vital role in OA cartilage degeneration and synovial inflammation. The critical role of the JAK-STAT signaling pathway in the senescence of cartilage-derived stem/progenitor cells (CSPCs) has been demonstrated, and the senescence of CSPCs may promote the development of OA.³⁴ The articular cartilage microenvironment is in a hypoxic state, and the hypoxia-inducible factor (HIF-1α) signaling pathway also has an important role in the development of OA and is critical to tissue homeostasis.³⁵ Cell senescence is a state of cell cycle arrest induced by cell stress and is related to the occurrence of OA.^19,20 Therefore, we infer that 16 differentially expressed cellular senescence-related genes related to OA mediate the occurrence of OA through the above pathways.

Then, five biomarkers related to cell senescence were finally determined through machine learning, including BCL6, MCL1, SLC16A7, PIM1, and EPHA3. BCL6 (BCL6 transcription inhibitor) is related to the aging of human fibroblasts.³⁶ The findings of our investigation revealed a decrease in BCL6 levels in individuals with osteoarthritis. BCL6, a member of the POZ/BTB-zinc finger protein family, has been demonstrated to enhance the production of chondrocytes and plays a crucial role in controlling the differentiation of osteoblasts and maintaining bone balance. Deficiency of BCL6 results in heightened osteoclast genesis.^37,38 MCL1 was significantly reduced in OA in this study. MCL1 (MCL1 apoptosis regulator, a member of the BCL2 family) is an anti-apoptotic molecule essential for the survival of lymphocytes, macrophages, and synovial fibroblasts.³⁹ The study demonstrated that reducing MCL1 expression hinders the PI3K/Akt signaling pathway, enhances the inflammatory response, and induces death in chondrocytes affected by osteoarthritis.^40,41 PIM1, a member of the pim gene family, is a crucial regulator of apoptosis, differentiation, proliferation, and tumorigenesis. Pim-1 enhances cell proliferation and differentiation by directly phosphorylating Cdc25C to facilitate the G2/M transition of the cell cycle.^42,43 PIM diminishes the extent of airway inflammation in mice and enhances the viability of Human bronchial epithelial BEAS-2B cells.⁴⁴ Various studies indicate that PIM1 is linked to cellular senescence, PIM1 preserves the integrity and structure of mitochondria by increasing the levels of anti-apoptotic proteins. This process may promote cell proliferation, survival, metabolic activity, and counteract senescence.^45–48 However, there is a lack of relevant studies on Pim-1 in osteoarthritis. Our study revealed that PIM1 enhanced the viability, proliferative capacity, and migratory ability of RPKFCL cells and reduced cellular senescence. EPHA3 has been reported to promote a variety of tumours, however osteoarthritis is unclear.^49,50 However, the relationship between these five cellular senescence-associated genes and OA has been poorly reported and further experimental studies are needed.^51–53 We report for the first time that BCL6, MCL1, SLC16A7, PIM1, and EPHA3 are biomarkers associated with OA cellular senescence.

Previous studies have shown that aging is closely related to cellular immunity, and that deficiencies in the immune system accelerate the aging process.^54,55 The single-sample GSEA results showed that these five genes might participate in immune-related pathways in the occurrence and progression of OA. Through the CIBERSORT algorithm, we found that T cell CD4 memory quiescence, NK cell activation, dendritic cell activation, mast cell activation, and neutrophil decrease, and B cell memory, plasma cell, and mast cell quiescence increase may be related to the occurrence of OA. Correlation analysis showed that BCL6 and MCL1 were significantly associated with NK cell activation. BCL6 was significantly associated with dendritic cell activation; BCL6, MCL1, SLC16A7, and PIM1 were significantly associated with mast cell activation. Previous studies have shown that NK cells and neutrophils play an essential role in experimental osteoarthritis and that the CXCL10/CXCR3 axis facilitates the interaction between the two.⁵⁶ Dendritic cells play a vital role in the development and progression of inflammatory arthritis (IA).⁵⁷ Nie et al find that increased Toll-like receptor expression in OA-derived dendritic cells contributes to inflammatory conditions.⁵⁸ Increasingly studies have shown that mast cells play an important role in the pathophysiology of OA, including inflammation and cartilage destruction.^59,60 However, further studies are needed to elucidate the interaction between these five cellular senescence-related biomarkers and the above immune cells.

The results of ROC curves showed that the combination of five biomarkers could clearly distinguish OA samples from normal samples. Therefore, we constructed a column line graph OA diagnostic model based on the expression of 5 genes. The expression amount of five biomarkers can score this nomogram, and this score can indicate the prevalence probability of the sample. The calibration curves indicate that the predictive accuracy of the column line graphs is high and can provide guidance for the diagnostic evaluation of OA.

We also constructed a lncRNA-miRNA-mRNA regulatory network of 5 biomarkers, which provides direction and basis for further study of the regulatory mechanisms of these 5 genes in OA. Current clinical medications for OA treatment, such as oral NSAIDs and glucosamine, are only used for pain relief and do not prevent or intervene in the progression of the disease, and thus treat it.^61,62 We have identified five biomarkers and predicted 26 targeted drugs for potential targeted therapeutics. These drugs include AZD-3965, Imatinib, ABL-001, ARQ-087, BMN-673, AT-9283, Entosplenitib, Vemurafenib, Tandutinib, Brigatinib, Melphalan, Okadaic acid, ABT-199, AMONAFIDE, PD 184352, MPC-3100, CH-5132799, Copanlisib, Bosutinib, 5-aza-T-dCyd, Ulixertinib, Apitolisib, ABT-737, GDC-0084, Vorinostat, and GSK-2126458. Imatinib exhibits anti-inflammatory properties and inhibits synovial fibroblast proliferation.^63,64 PD184352 acts as an anti-inflammatory and antioxidant, and it reduces osteoarthritis symptoms in mice. Vorinostat inhibits IL-6 signaling in chondrocytes, which is relevant to osteoarthritis.⁶⁵ The relationship between the remaining 23 drugs and osteoarthritis is unclear, and further research is needed to explore their potential as targeted therapies for osteoarthritis.

Finally, we conducted RT-PCR to confirm the expression of BCL6, MCL1, SLC16A7, PIM1, and EPHA3 in the synovial tissue of osteoarthritis and healthy controls. Our findings are consistent with the results of the GEO dataset, where the expression levels of EPHA3 were significantly elevated in osteoarthritis, and the levels of BCL6, MCL1, SLC16A7, and PIM1 were significantly downregulated compared to healthy controls. Our results confirm the potential of the above genes as diagnostic markers for osteoarthritis.

Our study focused on bioinformatics analysis, and although five genes were validated by qRT-PCR, the relatively insufficient sample size may lead to bias, and future studies need to be supported and validated by more clinical samples and data. Second, the mechanism of action of the identified cellular senescence-related biomarkers in the onset and progression of OA requires further experimental studies. The clinical application value of the predicted drugs targeting the five biomarkers needs to be validated in further clinical trials.

Conclusion

We identified, for the first time, cellular senescence-related biomarkers in OA patients. We established corresponding diagnostic models, which can provide a reference for the clinical diagnosis and treatment of OA patients, as well as a basis for further understanding of the pathogenesis of OA.