Identification of aging-related biomarkers and immune infiltration characteristics in osteoarthritis based on bioinformatics analysis and machine learning

Background Osteoarthritis (OA) is a degenerative disease closely related to aging. Nevertheless, the role and mechanisms of aging in osteoarthritis remain unclear. This study aims to identify potential aging-related biomarkers in OA and to explore the role and mechanisms of aging-related genes and the immune microenvironment in OA synovial tissue. Methods Normal and OA synovial gene expression profile microarrays were obtained from the Gene Expression Omnibus (GEO) database and aging-related genes (ARGs) from the Human Aging Genomic Resources database (HAGR). Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), Disease Ontology (DO), and Gene set variation analysis (GSVA) enrichment analysis were used to uncover the underlying mechanisms. To identify Hub ARDEGs with highly correlated OA features (Hub OA-ARDEGs), Weighted Gene Co-expression Network Analysis (WGCNA) and machine learning methods were used. Furthermore, we created diagnostic nomograms and receiver operating characteristic curves (ROC) to assess Hub OA-ARDEGs’ ability to diagnose OA and predict which miRNAs and TFs they might act on. The Single sample gene set enrichment analysis (ssGSEA) algorithm was applied to look at the immune infiltration characteristics of OA and their relationship with Hub OA-ARDEGs. Results We discovered 87 ARDEGs in normal and OA synovium samples. According to functional enrichment, ARDEGs are primarily associated with inflammatory regulation, cellular stress response, cell cycle regulation, and transcriptional regulation. Hub OA-ARDEGs with excellent OA diagnostic ability were identified as MCL1, SIK1, JUND, NFKBIA, and JUN. Wilcox test showed that Hub OA-ARDEGs were all significantly downregulated in OA and were validated in the validation set and by qRT-PCR. Using the ssGSEA algorithm, we discovered that 15 types of immune cell infiltration and six types of immune cell activation were significantly increased in OA synovial samples and well correlated with Hub OA-ARDEGs. Conclusion Synovial aging may promote the progression of OA by inducing immune inflammation. MCL1, SIK1, JUND, NFKBIA, and JUN can be used as novel diagnostic biomolecular markers and potential therapeutic targets for OA.


Introduction
Age is the paramount risk factor for osteoarthritis (OA), which is one of the most common causes of chronic pain and disability in the elderly. The prevalence of OA is rising due to an increase in the number of elderly and obese people, negatively impacting patients' quality of life and imposing a massive burden on families and society (1,2). The pathogenesis of OA is extremely complex, involving mechanical overload, an increase in inflammatory mediators, metabolic changes and cellular senescence, all of which can interact to promote the development of OA (3). As a result, studying the molecular biology of OA and looking for potential biomarkers of OA is critical for early diagnosis and treatment of OA.
Aging is a complicated biological process. Senescent cells continue to accumulate in the human body as we age, leading to a variety of age-related diseases such as osteoarthritis, cardiovascular disease, Alzheimer's disease, and so on (4)(5)(6). Cellular senescence is one of the first signs of aging, characterized by permanent cell cycle arrest and the release of harmful proinflammatory factors into the surrounding microenvironment, a feature known as Senescence-associated secretory phenotype, SASP (7). The persistence of SASP factors causes chronic, low-grade systemic inflammation, stimulates senescence in neighboring cells and accelerates aging progression (8). SASPs such as IL-1, IL-6, MMP-13, VEGF and other pro-inflammatory factors have been found in OA cartilage and synovial fluid (9). Many synovial fibroblasts positive for p16 and SA-b-Gal have been found in the synovium of elderly and OA patients. Senescent synovial fibroblasts can aggravate synovial inflammation and cause cartilage degradation by secreting pro-inflammatory factors and matrix metalloproteinases (10). Intriguingly, injection of senescent fibroblasts into the knee joint cavity of mice resulted in inflammation, cartilage erosion and osteophyte formation (11). In post-traumatic mouse models, targeted removal of senescent cells in mouse knee articular cartilage and synovium specifically marked by P16 can effectively reduce the development of inflammatory response and articular cartilage injury-related pain (12). Aging is clearly involved in multiple pathways of OA pathogenesis, but the precise mechanism remains unknown.
Bioinformatics is a multidisciplinary field. With the development of high-throughput sequencing technology and the application of machine learning in the medical field in recent years, new ideas for studying the molecular mechanisms of various diseases have emerged (13,14). WGCNA and three machine learning algorithms were applied in this study to screen the Hub OA-ARDEGs of OA and an online database to predict their potential miRNAs and TFs. The ssGSEA algorithm was often used to investigate OA's immune infiltration profile and its relationship with Hub OA-ARDEGs. It provides a new direction for the early diagnosis and treatment of OA. The overall workflow of this study is depicted in Figure 1.

Gene expression dataset screening and processing
We downloaded the gene expression profile microarrays (GSE55235, GSE55457, GSE12021, GSE1919, GSE89408) for normal and OA synovial samples from the Gene Expression Omnibus (GEO), with the information shown in Table 1. Probe names were converted to gene names using Perl and with the help of platform annotation files. Background correction and normalization of each dataset using the R package limma and integrating three synovial datasets from the same platform using the R package sva to remove batch effects (15). Two-dimensional PCA clustering plots were used to show the before and after differences in removing batch effects from the samples. For subsequent analysis, the merged dataset served as the training set and the GSE1919 and GSE89408 dataset served as the validation set. (16). After merging and deleting duplicate genes, 543 ARGs were obtained for subsequent analysis. The list of aging-related genes is detailed in Table S1 of the Supplementary Material.

Identification of aging-related differentially expressed genes
ARGs expression matrices were extracted from the training set and analyzed for differences using the R package limma, with | logFC|>0.5; FDR<0.05 as the criterion for screening to obtain aging-related differential genes (ARDEGs) in normal and OA synovial samples (17). A heat map of gene expression was created to visualize the top 30 genes with the most significant differences.

Construction of protein-protein interaction networks network
To evaluate the gene interactions among the ARDEGs, a protein-protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING, Flowchart for Comprehensive Analysis of Aging-Related Genes in Osteoarthritis. Zhou et al. 10.3389/fimmu.2023.1168780 Frontiers in Immunology frontiersin.org https://cn.string-db.org/) database (18), with a confidence score of >0.7 as the cut-off criterion.

Functional enrichment analysis
R package Clusterprofiler and DOSE were used to perform Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Disease Ontology (DO) enrichment analysis of ARDEGs, with q-value<0.05 set as a screening criterion to investigate their biological functions, signaling pathway enrichment, and disease similarities.

Gene set variation analysis
Gene set variation analysis (GSVA) is an unsupervised, nonparametric method for assessing transcriptome gene set enrichment (19). We used the gene sets "Hallmark.all.v2022.1.Hs.symbols" and " c2.cp.kegg.symbols " from the Molecular Signature Database (MSigDB) as reference sets. To assess the enrichment of normal and OA synovium samples, the R package GSVA was used to score the HALLMARKS and KEGG pathways. Significant enrichment was defined as FDR<0.05.

WGCNA and screening for ARDEGs with highly correlated OA features
Weighted Gene Co-expression Network Analysis (WGCNA) is an algorithm that finds biologically significant co-expressed gene modules and explores the relationship between gene networks and disease (20). We used the R package WGCNA to build a weighted co-expression network on the training set, analyzing the genes with the highest 50% expression variance in all expression profiles. The goodSamplesGenes function was used to check for missing values in the data, and the "pickSoftThreshold" function was used to filter and validate the best soft threshold. The Pearson correlation coefficient was used to create the adjacency matrix, which was then transformed into a topological overlap matrix (TOM) with appropriate power values and phase anisotropy (1-TOM). TOM classified the genes into different modules. The genes most associated with OA traits in the module were selected as cor.MM>0.7 and cor.GS>0.5. 0.5 is the screening condition for Hub genes in the module, and the intersection of Hub genes and ARDEGs in the module is taken to obtain OA-ARDEGs. The Pearson correlation test assessed the interaction of OA-ARDEGs at the mRNA level.

Identification of Hub ARDEGs with highly correlated OA features
LASSO regression, SVM-RFE, and random forest are methods for screening and identifying Hub OA-ARDEGs. The least absolute shrinkage and selection operator (LASSO) regression is a common data mining method (21). R package glmnet was applied to incorporate OA-ARDEGs into the diagnostic model, set the a value of the glmnet function to 1, obtained the best l by ten cross-validations finally obtained the Aging signature genes based on the best l. Support vector machine recursive feature elimination (SVM-RFE) is a common machine learning method based on embedded methods (22). The R package e1071 helped us find the best variables and remove the feature vectors generated by the SVM, thus obtaining the Aging signature genes. Recursive Feature Elimination (RFE) of the Random Forest algorithm is a supervised machine learning method (23). Aging signature genes were identified by relative importance greater than one when the decision tree was set to 500. The intersection of the three machinelearning filtered Aging signature genes was defined as Hub OA-ARDEGs using R package Venn. The GSE1919 and GSE89408 dataset could validate the receiver operating characteristic curve (ROC) analysis of the diagnostic value of Hub OA-ARDEGs for OA.

Patients samples
The synovial tissue from six patients who underwent total knee replacement surgery for osteoarthritis (OA) and normal synovial tissue from six patients with meniscus injuries were obtained for this study. All patients signed informed consent forms, and the collection, processing, and analysis of the samples were conducted under the guidance of the Ethics Committee of the Guangzhou Red Cross Hospital, affiliated with Jinan University (Ethics Approval No. 2023-001-01).

qRT-PCR
The total synovial tissue RNA was extracted using Trizol (Servicebio), followed by reverse transcription of the total RNA into complementary DNA (cDNA) using Takara Prime Script ® RT Master Mix. Quantitative real-time polymerase chain reaction (qRT-PCR) was conducted using 2× SYBR Green qPCR Hub Mix (without ROX) (Service). The primer sequences for the Hub osteoarthritis-associated differentially expressed genes (Hub OA-ARDEGs) can be found in Table S16 of the Supplementary Material. The GAPDH gene was utilized as an internal reference gene. Each biological sample was subjected to three technical replicates.

Construction of Hub OA-ARDEGs risk prediction model
To improve clinical applicability, we use the R package rms to create a nomogram with Hub OA-ARDEGs, where "Points" represents the score of candidate genes and "Total Points" represents the sum of the scores of all the genes listed above. The accuracy of the nomogram model was determined by calibration, clinical decision, and Clinical impact curve.

miRNA-TF-mRNA regulatory network of Hub OA-ARDEGs
To improve prediction accuracy, we used the miRTarBase (24), Starbase (25), and Targetscan (26) databases to predict the miRNAs of Hub OA-ARDEGs. The Enrichr database (http://amp.pharm.mssm.edu/ Enrichr/) was also applied to predict the transcription factors (TF) of Hub OA-ARDEGs, with a p-value of 0.05 as a screening condition. Construction of miRNA-TF-mRNA regulatory networks and visualization of the networks using Cytoscape (3.9.1).

Immunological characteristics of OA
Single sample gene set enrichment analysis (ssGSEA), an extension of the GSEA method, is widely used in bioinformatics studies related to immune infiltration (27). We calculated enrichment scores for normal and OA synovial samples in 28 immune cells and 13 immune functions using the R package GSVA and visualized the results using the R package vioplot. Spearman correlation analysis was used to correlate Hub OA-ARDEGs with immune cells and immune function.

Statistical analysis
All statistical analyzes were performed in R (version 4.2.2). Comparisons among two groups were made using Wilcox test.
Spearman's correlation analysis was used to understand the relationship between the expression levels of Hub OA-ARDEGs and immune cells and immune function. Differences were deemed statistically significant where P < 0.05. Statistical analyses of qRT-PCR were presented as the mean ± standard deviation for at least three individual experiments, and the statistical significance of differences was determined with the unpaired, two-tailed Student t-test. (*P < 0.05; **P < 0.01; ***P < 0.001). P < 0.05 was considered as statistically significant.

GEO data processing
We integrated three synovial datasets, GSE55235, GSE55457, and GSE12021, containing a total of 29 normal synovial and 30 OA synovial samples. As is shown in Figure 1 of the Supplementary Material, the gene expression level and principal component analysis (PCA) of each sample before and after eliminating the batch effect.

Identification and PPI analysis of ARDEGs
We identified 87 ARDEGs using the R package limma and screening criteria of |logFC|>0.5 and FDR<0.05, of which 32 genes were upregulated in OA and 55 genes were downregulated in OA. Table S2 of the Supplementary Material contains a detailed list of differentially expressed Aging-related genes. The heat map and volcano map were used to depict the differences (Figures 2A, B). PPI protein network interaction analysis revealed that ARDEGs interact closely at the protein level ( Figure 2C). The results of the PPI protein network interaction analysis are available in the Supplementary Material: Table S3.

Functional enrichment analysis of ARDEGs
To better understand the potential mechanisms of ARDEGs in OA, we performed GO, KEGG, and DO enrichment analysis on ARDEGs using the R package clusterProfiler. GO enrichment analysis revealed that the first five ARDEG enrichments were primarily involved in the response to extracellular stimulus, neuron death, gland development, response to nutrient levels, and regulation of neuron death. The top 5 enriched items in Cellular Components (CC) and Molecular Functions (MF) are shown in Figure 3A. Furthermore, KEGG pathway analysis revealed that these ARDEGs were enriched in the HIF-1 signaling pathway, the FoxO signaling pathway, Kaposi sarcoma-associated herpesvirus infection, the PI3K-Akt signaling pathway, and the MAPK signaling pathway, and the pathways interacted closely (Figures 3B, C). DO enrichment analysis reveals disease types with similar pathogenic mechanisms of ARDEGs in OA, such as prostate cancer, male reproductive organ cancer, hepatitis B, female reproductive organ cancer, and hepatitis C ( Figure 3D). Tables S4-6 show detailed results for GO, KEGG, and DO enrichment of ARDEGs.

GSVA enrichment analysis
We investigated HALLMARKS and KEGG pathway enrichment in OA through the GSVA method. TNFA signaling via the NFKB pathway, apoptosis, and the P53 pathway were significantly upregulated in OA when compared to the control group, according to HALLMARKS pathway enrichment results. At the same time, DNA repair, oxidative phosphorylation, and bile acid metabolism were all significantly reduced ( Figure 4A). As is shown in the KEGG pathway enrichment results, the first three pathways were significantly upregulated in small cell lung cancer, adipocytokine signaling pathway, and chronic myeloid leukemia in OA. In contrast, leishmania infection, the pentose phosphate pathway, and fatty acid metabolism were all significantly downregulated ( Figure 4B).

WGCNA
We used the R package WGCNA, expression variance in the first 50% of genes as a screening condition to eliminate less volatile genes, and 6538 genes for co-expression network construction. The value 18 was selected as the optimal soft threshold (R 2 = 0.9) to establish a scale-free network ( Figure 5A). Following that, cluster analysis was used to identify highly similar modules, with the minimum module size set at 60. Using dynamic hybridization shearing, eight gene modules were obtained, with one red module (217 genes) having the highest correlation (cor) with OA (cor = 0.86; P = 3e-18) ( Figure 5C). Furthermore, there was a strong correlation between GS and MM within the red module (cor = 0.72; P=5.8e-36) ( Figure 5D). The genes in the red module with cor.MM>0.7 and cor.GS>0.5 were chosen as the Hub genes and intersected with ARDEGs to yield 20 with OA-ARDEGs ( Figure 5E). The results of co-expression modules for datasets can be found in the Supplementary Material: Table S7.

Correlation and enrichment analysis of OA-ARDEGs
We evaluated the correlation between OA-ARDEGs by the Pearson correlation coefficient. MCL1 was found to be highly correlated with ZFP36 (cor=0.87) and BHLHE40 (cor=0.74) ( Figure 6A). Figure 6B depicts the chromosomal location of OA-ARDEGs. Supplementary Material: Table S8. KEGG pathway enrichment analysis revealed that the first five pathways enriched by OA-ARDEGs were primarily involved in Kaposi sarcoma-associated herpesvirus infection, IL-17 signaling pathway, Human Tcell leukemia virus 1 infection, AGER-AGE signaling pathway in diabetic complications, and TNF signaling pathway ( Figure 6C). GO enrichment analysis shows that OA-ARDEGs are enriched in biological processes such as regulation of transcription from RNA polymerase II promoter in response to stress, regulation of DNA-templated transcription in response to stress, myeloid cell ( Figure 6D). The specific results of KEGG and GO enrichment were shown in Supplementary Material: Table S9, 10.

Identification and validation of Hub OA-ARDEGs
To improve the accuracy of Hub OA-ARDEGs diagnostic OA, we used three machine learning algorithms, LASSO ( Figures 7A, B), SVM-RFE ( Figures 7C, D), and Random Forest ( Figures 7E, F), to screen OA-ARDEGs. After combining the results of the three algorithms, a total of five Hub OA-ARDEGs were obtained, namely, MCL1, SIK1, JUND, NFKBIA, and JUN ( Figure 7G). Results of three machine learning identification Hub OA-ARDEGs are shown in the Supplementary Material: Table S11.

Construction of Hub OA-ARDEGs risk prediction model
We developed a diagnostic nomogram for OA based on the expression of Hub OA-ARDEGs in order to obtain a more clinically applicable diagnostic model for OA. By constructing clinical calibration curves ( Figure 8B

Hub OA-ARDEGs expression and diagnostic value
Based on the expression levels of Hub OA-ARDEGs in the training set and validation set, we found that Hub OA-ARDEGs were significantly downregulated in all OA synovial samples ( Figure 9A-C). ROC curve analysis showed that the 5 Hub OA-ARDEGs and nomogram had high diagnostic values for OA in the training set. MCL1 and nomogram had the highest diagnostic value (AUC=1.000), and the other genes had the following diagnostic values: JUN (AUC=0.960), SIK1 (AUC=0.955), NFKBIA (AUC=0.976) and JUND (AUC=0.968) ( Figure 9D). Figures 9E, F displays the ROC analysis results for the external validation set GSE1919 and GSE89408. The AUC for all five Hub OA-ARDEGs and nomogram in the validation set were greater than 0.5. As a result, the five Hub OA-ARDEGs can be used as reliable biomarkers for the diagnosis of OA and have a high diagnostic value.

qRT-PCR
We performed qRT-PCR on total mRNA from synovial membranes of six patients with meniscal injuries and six patients with OA to further validate the mRNA expression levels of Hub OA-ARDEGs. The results showed that all five Hub OA-ARDEGs were significantly downregulated in the OA synovial samples (pvalue less than 0.05), consistent with the expression in the training set ( Figures 10A-E).

Construction of the miRNA-TF-mRNA regulatory network
By predicting miRNA and TF on Hub OA-ARDEGs, we used Cytoscape (3.7.1) to visualize the regulatory network, which contained 80 miRNAs, 12 transcription factors, and five genes, and obtained a total of 196 miRNA-TF-mRNA regulatory relationships ( Figure 11). Details of miRNA-TF-mRNA regulatory networks were shown in Supplementary Material: Table S13.

Immune infiltration analysis
We used the ssGSEA algorithm to discover that the infiltration levels of Activated B cell, Immature dendritic cell, Macrophage, Regulatory T cell, Central memory CD4 T cell, Memory B cell, and Effector memory CD4 T cell were significantly increased in OA samples. In contrast, the infiltration of Eosinophil, Type 2 T helper cell, and Central memory CD8 T cells in OA samples was significantly reduced (Figures 12A, C). The Spearman's correlation was applied to analyze the interaction between immune cells. The results revealed there are significant correlations between most immune cells, e.g. Macrophage and MDSC were significantly positively correlated (r = 0.89) ( Figure 12B). Checkpoint, HLA, parainflammation, T cell coinhibition, T cell co-stimulation, Type I IFN Response, and other immune functions were significantly activated in OA samples ( Figure 12D). HLA gene expression levels were higher in OA samples, including HLA-DMA and HLA-DRA ( Figure 12E), which indicates that the immune response plays an essential role in the development of OA. JUND, JUN, MCL1, NFKNIA, and SIK1 correlated well with a variety of immune cells and immune functions ( Figures 12F, G). For example, JUN was positively correlated with Activated CD4 T cells, Type 2 T helper cells, and aDCs, and negatively correlated with Macrophage, Mast cells, Memory B cells, and Checkpoint. JUND was negatively correlated with Gamma delta T cells and APC co-inhibition. MCL1 was positively correlated with Natural killer T cell, Activated CD4 T cell, and negatively correlated with CD56 bright natural killer cell and APC co-inhibition. NFKBIA was negatively correlated with Activated B cell, Activated CD8 T cell, etc. SIK1 was positively correlated with Central memory CD4 T cell, Central memory CD8 T cell, and Cytolytic activity, and negatively correlated with Neutrophil, Regulatory T cell, Natural killer cell, Macrophage, APC co-inhibition, Type II IFN Response, and HLA. Data of results were shown in the Supplementary Material: Table S14, 15.

Discussion
OA is a chronic disease characterized by pain, cartilage loss, and joint inflammation, with aging playing a significant role in its progression. To date, OA cartilage has been extensively studied in terms of aging. Oxidative stress, for example, causes chondrocyte senescence by increasing p53 and p21 expression and activating the p38 MAPK and PI3K/Akt signaling pathways (28). Mechanical stress causes chondrocyte senescence and accelerates cartilage catabolism by downregulating FBXW7 (29). Sirt6 slows chondrocyte aging by inhibiting IL-15/JAK3/STAT5 signaling (30). However, a growing body of research has recently recognized that synovial aging plays an important role in OA, but the exact role and mechanisms remain unknown (31). As a result, the purpose of this study was to identify potential biomarkers of aging in OA and to investigate the role and mechanisms of Aging-related genes and immune infiltration in OA synovial tissues, thereby providing new directions and ideas for potential mechanisms and early OA diagnosis.
In normal and OA synovial samples, we discovered 87 ARDEGs. According to the GO and KEGG enrichment analysis results, ARDEGs are primarily involved in response to extracellular stimulus, response to nutrient levels, and HIF-1 signaling pathway, FoxO signaling pathway, PI3K-Akt signaling pathway, and MAPK signaling pathway, which is consistent with previous studies (32)(33)(34)(35).
It has been proposed that ARDEGs may be involved in inflammatory regulation, cellular stress response, cell cycle regulation, transcriptional regulation, and other mechanisms that promote OA progression. The results of the GSVA enrichment analysis revealed that TNFA signaling via the NFKB pathway, apoptosis, MAPK signaling pathway, and the P53 pathway were all significantly upregulated in OA. This is consistent with the functional enrichment of ARDEGs in OA. Using WGCNA analysis and three machine learning screenings, we were able to identify five Hub OA-ARDEGs (MCL1, SIK1, JUND, NFKBIA, JUN). Our findings suggested that Hub OA-ARDEGs had an excellent diagnostic ability to predict OA and were significantly downregulated in synovial samples from OA patients. MCL1 (MCL1 Apoptosis Regulator, BCL2 Family Member) is an anti-apoptotic member of the BCL-2 family of proteins, which is involved in the regulation of apoptosis, cellular senescence, and inflammation and is essential for the maintenance of cell survival and viability (36). MCL1 expression is significantly downregulated in OA and senescent chondrocytes, and miR-34a-targeted inhibition of MCL1 can induce chondrocyte apoptosis as well as promote inflammatory response and ECM degradation (37). RHEB overexpression can, surprisingly, upregulate MCL1 to alleviate chondrocyte senescence and oxidative stress, which may be related to MCL1 inhibiting ROS production and P27 expression (38). MCL1 downregulation in OA synovial tissue and a positive correlation with natural killer cells suggest that MCL1 downregulation may be involved in synovial senescence and synovitis. However, more experimental support is required. SIK1 (Salt Inducible Kinase 1) is a member of the AMPK kinase subfamily that regulates cell cycle, gluconeogenesis, and lipogenesis (39). SIK1 has been shown in studies to inhibit TLR4induced NF-B activation and reduce the expression of proinflammatory cytokines (40). However, SIK1 has not been reported in the OA literature. NFKBIA (NFKB Inhibitor Alpha) is an NFKB inhibitor that reduces the inflammatory response by inhibiting the activity of the dimeric NF-kappa-B/REL complex. Multiple studies have shown that abnormal NF-B activation is associated with chondrocyte catabolism, chondrocyte senescence, and synovial inflammation, and NF-B inhibition may be a potential therapeutic target for the treatment of OA (41,42). Interestingly, the number of chondrocytes activated by NF-B signaling was significantly increased in aged articular cartilage, and activation of IKK-NF-B signaling in chondrocytes accelerated the occurrence of age-related joint tissue degeneration (43). Overexpression of NFKBIA in OA synovial fibroblasts effectively inhibits the expression of destructive enzymes Hub OA-ARDEGs expression difference and ROC curve. such as MMP-1, MMP-3, and MMP-13, as well as ADAMTS4, reducing inflammation (44). These findings support our hypothesis that NFKBIA downregulation is important in promoting OA. JUND (JUND Proto-Oncogene, AP-1 Transcription Factor Subunit) is a functional component of the AP1 transcription factor complex, which is involved in cell proliferation, differentiation, and senescence, primarily through the regulation of oxidative stress levels (45,46). JUND has been shown to protect cells from p53-dependent senescence and apoptosis, making it an appealing molecular target for preventing or delaying age-related cardiovascular disease (47,48). In osteoarthritis, JUND transcriptional activation acts on the LncRNA LOXL1-AS1 to promote chondrocyte proliferation and inflammation, leading to osteoarthritis progression (49). The precise role of JUND in OA has yet to be determined. JUN (Jun Proto-Oncogene, AP-1 Transcription Factor Subunit) is a nuclear transcription-activating protein l family member involved in a variety of physiological processes such as cell cycle progression, differentiation, and apoptosis (50). JUN was discovered to slow aging by suppressing p53 gene transcription, whereas JUN knockout mouse embryonic fibroblasts (MEF) exhibit severe proliferation defects and early senescence (51). JUN binds to BATF to form a complex that regulates the expression of anabolic and catabolic genes in chondrocytes, which is critical in the destruction of OA cartilage (52). JUN overexpression has also been linked to reduced SOX9 transcriptional activity and type II collagen expression in chondrocytes (53). We discovered significant JUN downregulation in OA synovium, but it is unclear whether this promotes synovial aging and inflammation. More research into the specific role and mechanisms of Hub OA-ARDEGs in OA will hopefully lead to an emerging target for OA treatment. Long-term low-level chronic inflammation and innate and adaptive immune system activation play critical roles in all aspects of OA pathogenesis (54, 55). Using ssGSEA analysis, we discovered that macrophages, natural killer cells, regulatory T cells, activated CD8 T cells, central memory CD4 T cells, effector memory CD4 T cells, activated B cells, and memory B cells were significantly FIGURE 11 miRNA-TF-mRNA regulatory network. Diagram of miRNA-TF-mRNA regulatory network, where red circles represent genes; blue V-shapes represent predicted miRNAs; green diamonds represent TF. infiltrated in OA. Checkpoint, HLA, parainflammation, and other immune functions were significantly activated in OA samples. Macrophages, the primary innate immune cells in the OA synovium, are present in 76% of OA patients' knee joints and are significantly associated with knee pain, OA radiographic severity, and osteophytes (56). Through the secretion of inflammatory, growth factors and MMPs, activated macrophages cause chondrocyte senescence and apoptosis (57). NK cells, as one of the innate immune cells, play an important role in monitoring and killing senescent cells. In the liver, NK cells use perforin granule exocytosis to target the clearance of senescent hepatic stellate cells (58). The recruitment and activation of NK cells is thought to accelerate the progression of OA in a mouse model of osteoarthritis with cartilage damage (59). Similarly, the adaptive immune system unquestionably plays a role in OA. Numerous studies have found a significant infiltration of T cells in OA synovial and synovial fluid, second only to macrophages in accounting for 25% of inflammatory cells (60). CD8+ T-cell-induced tissue inhibitor of metalloprotein-1 (TIMP1) expression aggravates osteoarthritis (61). The number of synovial CD4+ T lymphocytes correlates with a visual analog pain scale, promotes Th1 cell polarization, and increases the release of immunomodulatory cytokines, accelerating the progression of OA (62,63). B cells and plasma cells have been found in synovial tissue from OA patients. However, no studies have found a direct link between B cell infiltration and OA progression or severity, and further experimental confirmation of their specific role is required (64). We used Spearman correlation analysis to show that the five Hub OA-ARDEGs have a reasonable correlation with immune cells and functions, implying that synovial immune inflammation may collaborate with aging to contribute to the occurrence and progression of OA. However, there are limitations to this study: 1. The transcriptomic data for this study were sourced from publicly available databases, which may limit access to more clinically relevant information. The variability of patient populations and clinical characteristics cannot be ruled out as a possible influence on the study's findings. 2. The limited sample size had a potential impact on the accuracy of the results, emphasizing the need for a larger sample size and a prospective study design to validate and reinforce the model's findings. 3. We conducted an extensive exploration and analysis of databases and supplemented our findings with limited experimental validations. However, it is important to note that our study would benefit from additional experimental support to further strengthen the conclusions.

Conclusion
In conclusion, this study preliminarily investigated the potential mechanisms of aging-related genes in OA synovial tissue, which revealed that synovial aging could be closely linked to immune inflammation. In addition, five Hub OA-ARDEGs have excellent diagnostic capabilities for OA and may be novel targets for the diagnosis and treatment of OA. However, additional experimental studies are required to support our findings.

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 authors.

Ethics statement
The studies involving human participants were reviewed and approved by Guangzhou Red Cross Hospital of Jinan University (Ethics number 2023-001-01). The patients/participants provided their written informed consent to participate in this study.

Author contributions
QM and LW directed the research and revised the manuscript. JZ and JH conceived the idea and wrote the manuscript, they have contributed equally to this work and share first authorship. ZL and QS are responsible for analyzing and processing the data. ZY processed and modified the image. All authors contributed to the article and approved the submitted version.

Funding
The major project of Bijie Bureau of Science and Technology (Grant numbers [2022]-1) to QM.; the Science and Technology of Program of Guangzhou (Grant numbers 202102080344) to QM.