Identification and Analysis of Neutrophil Extracellular Trap-Related Genes in Osteoarthritis by Bioinformatics and Experimental Verification

Introduction

OA is a degenerative joint disease that primarily affects the elderly population. It is a multifactorial disorder with a complex pathogenesis, involving a variety of joint tissues. In addition to the well-established degradation of articular cartilage, OA encompasses a comprehensive joint pathology, encompassing the synovial membrane, subchondral bone, menisci, and infrapatellar fat pad. These joint components play pivotal roles in preserving joint homeostasis and function. Their alterations profoundly influence biomechanical changes, thus exerting a significant impact on the overall disease progression and its severity.^1–4 The cartilage in the joints continues to suffer wear and tear as it functions as a cushion, which over time may lead to pain, stiffness, and limited range of motion.⁵ The following factors are strongly associated with osteoarthritis: age, obesity, genetic predisposition, and joint injury.⁶ In OA patients, the inflammatory response in the synovial membrane is driven by the release of damage-associated molecular patterns (DAMPs) during the degradation and breakdown of articular cartilage.^7,8 The inflammatory response leads to cellular damage and tissue injury, further accelerating the destruction of joint structures. Older adults with OA experience difficulty with daily activities, decreased independence, and increased risk of falls and fractures. Chronic joint pain could also lead to a decrease in sleep quality and even anxiety, further reducing their physical function.⁹ However, there is a lack of clinically effective treatment options for OA, and the molecular mechanisms and triggering factors of the disease are not well understood.

Upon entry of bacteria or viruses into the body, neutrophils receive signals and activated. Activated neutrophils release Peptidyl arginine deiminase 4 (PAD4) into the nucleus, which induces the release of DNA and granulocyte material from the nucleus to the extracellular to form NETs.¹⁰ NETs have different roles in various diseases and are similar to a double-edged sword. The NETs were initially thought to be a special mechanism for protecting the host from pathogens and destroying them.¹¹ However, NETs have recently been found to increase the inflammatory response in autoimmune diseases such as rheumatoid arthritis and systemic lupus erythematosus.¹² A significant increase in polymorphonuclear neutrophil [PMN]-NETs was observed in RA patients.¹³ Neutrophils in the synovial fluid of patients with RA are accompanied by high levels of expression of peptidylarginine deiminases during disease onset.¹⁴ Recent bioinformatics-based studies have found that genes associated with NETs can be of clinical value as prognostic for patients with pan-cancer and ischemia reperfusion injury, and have an intimate relationship with NETs.^13,15 Yet, comprehensive and in-depth studies on OA and NETs are still lacking.

Our study attempts to comprehensively analyze the subtle relationships between OA, NRGs and NETs-related signaling pathways (Neutrophils pathway, Neutrophil degranulation and Neutrophil granule constituents pathway). Based on OA-related datasets in the GEO database and NRGs collected in the literature, DE-NRGs was obtained after analysis. The GSVA algorithm was used to observe changes in NETs-related signaling pathways in OA patients. The CIBERSORT algorithm to evaluate the percentage of immune cell infiltration in OA patients. Next, a machine learning algorithm was used to determine the strongly correlated DE-NRGs and a Nomogram was constructed to help the clinical diagnosis and treatment of OA. Finally, we examined the previous results through animal experiments (mouse OA model), human OA samples and in vitro cellular experiments. In this study, we systematically analyzed the relationship between OA and NETs to help clinicians to diagnose and treat OA based on NETs as an entry point.

Materials and Methods

Data Collection

GSE55235 (con=10, OA=10), GSE55457 (con=10, OA=10), GSE98918 (con=12, OA=12) and GSE117999 (con=10, OA=10) downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/) for subsequent analysis.To ensure the robustness of the analysis, the four datasets were combined and batch effects were removed before further investigation.NRGs (n=69) were collected through the relevant literature,¹⁵ with specific information in Supplementary Table S1.

Identification of DE-NRGs in Patients with OA

Combine GSE55235, GSE55457, GSE98918 and GSE117999 and remove the batch effect. The “limma” package was used to identify differentially expressed genes (DEGs) in OA patients, with Adjusted p-value< 0.05 and |logFC| > 0.585 as screening criteria.¹⁶ After intersecting DEGs with NRGs, we successfully obtained the DE-NRGs.

Construction of Nomogram Model Based on DE-NRGs

Construction of a Nomogram DE-NRGs using the “rms” R Package. We evaluate each gene based on the I nfluence of DE-NRGs on the outcome variable, and then add up the evaluations to obtain the overall score. Finally, the total score was followed to provide assistance in the clinical diagnosis of OA patients.

Receiver Operating Characteristics (ROC) Analysis of DE-NRGs

To verify the usefulness of the Nomogram constructed by DE-NRGs for clinical diagnosis, we used the “pROC” package to plot ROC curves to evaluate the accuracy and usefulness of the Nomogram model.¹⁷

CIBERSORT Assessment of Immune Cell Infiltration in OA

The CIBERSORT algorithm was applied to determine the difference in immune cell counts between the control and OA groups, and the results were represented using a violin plot.¹⁸

Gene Set Variation Analysis (GSVA)

GSVA algorithm can evaluate the biological pathways or genomic activity changes in the disease, which is important for the diagnosis and prognosis of the disease.¹⁹ We applied this technique to analyze changes in biological pathways associated with NETs in OA patients.

Identification of Core DE-NRGs

Five core DE-NRGs were filtered by the results of the LASSO, SVM-RFE and RF algorithms.²⁰ The Lasso method involves a penalty applied to the absolute values of the coefficients in order to regulate their magnitude, leading to improved control against overfitting and heightened interpretability of the model.²¹ The SVM-RFE algorithm is used to identify the most important signature genes related to a disease by eliminating noisy or redundant features, reducing the dimensionality of the data and improving the accuracy of the model.²² The RF algorithm for regression is a machine learning technique that builds decision trees using random subsets of data and features. This reduces the correlation between trees and enhances the overall performance.

PPI Network (Protein-Protein Interaction) and Functional Enrichment Analysis

GeneMANIA (http://www.genemania.org) is a database that integrates information on genetic interactions, physical interactions, co-expression, co-localization, and gene enrichment analysis.²³ GeneMANIA can be utilized for construct PPI network and to generate gene function predictions and identify genes with comparable effects. The analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was conducted using the R package “clusterprofiler”. The Gene Set Enrichment Analysis (GSEA) algorithm was used to explore the biological changes that were significantly different between OA patients (median TLR7 expression was used to divide OA patients into two groups).²⁴

Collection of Human Tissue Samples

Human cartilage was obtained from patients undergoing total knee arthroplasty at the First Affiliated Hospital of Chongqing Medical University (Chongqing, China). Human OA cartilage was obtained from osteoarthritis patients undergoing total knee arthroplasty at the First Affiliated Hospital of Chongqing Medical University (Chongqing, China) (n=6). The medial part of the tibial plateau was used as the OA group, while the lateral tibial plateau with less cartilage destruction was used as the normal group (n=6). Collected human biological materials were approved by the Ethics Committee of the First Affiliated Hospital of Chongqing Medical University (2023.4.03 No: K2023-087). Patients provided informed consent following the principles of the Declaration of Helsinki.

Experiment OA in Mice

Male C57BL/6J mice Aged 8–10 weeks (25–30g) were purchased from the Animal Center of Chongqing Medical University. Animals were provided with sufficient food and water, housed in a stable environment at 22–26 degrees Celsius, and maintained on a regular circadian rhythm. The animal experiments covered in the article were in strict accordance with the guidelines of the Chongqing Medical University Animal Care and Use Committee (2022 12.06 No: 2022–620). The animal experiments were approved by the Animal Care and Use Committee of Chongqing Medical University. We randomly divided the mice into two groups (n=5), the Sham group and the DMM group. Destabilization of the medial meniscus (DMM) surgery was used to induce OA in mice. The sham group was performed on the knee with no ligament transection. The mice were fed the same diet for 8 weeks after surgery until they were sacrificed and the joints were collected for histology.

Histological Assessment

Knee cartilage destruction was assessed by safranin-O staining and scored using the OARSI grading system.²⁵ In brief, knee joints from mice and cartilage from human were fixed in 4% paraformaldehyde (PFA), decalcified in 14% EDTA, embedded in paraffin, and sectioned at 5-um thickness. The sections were deparaffinized in xylene, hydrated with graded ethanol, and finally stained with safranin O- fast green (Solarbio, Beijing, China).

Immunofluorescence

Joints from mice were fixed in 4% PFA for 14 h and decalcified in 14% EDTA for 3 d. Furthermore, the cartilage from human was also fixed in 4% PFA for 16 h and decalcified in 14% EDTA for 5 d. Afterward, the joints and cartilage were preserved in an optimum cutting temperature (OCT) compound, cut to a thickness of 10μm, and subjected to immunofluorescence using antibodies recognizing Toll-like receptor 7 (Boster, China, 1:100).The immunofluorescence protocol can be found in our previously published article.²⁶

Cell Culture

The mouse chondrogenic cell line ATDC5 cells were seeded into Dulbecco's modified Eagle’s medium (DMEM, Gibco, USA) which was supplemented with 10% fetal bovine serum (Gibco, USA), penicillin (100 units/mL), and streptomycin (100 μg/mL). Following reaching confluency, the ATDC5 cells were pretreated for 24 hours with tumor necrosis factor ⍺ (TNF-⍺, Abcam, USA) (20 ng/mL) or phosphate buffer saline (PBS).

Quantitative Reverse Transcription Polymerase Chain Reaction (RT-qPCR)

The process of mRNA and cDNA was conducted following the previously described method.²⁶ Primers of interest genes were designed using NCBI Primer-Blast software. Primer sequences in Supplementary Table S2.

Statistics

In this study, all expression matrix data were analysed using R software v4.1.2. The results of the experimental and control groups were analyzed by two-tailed Student’s t-test, and one-way ANOVA was performed for multiple comparisons (p < 0.05 = *p < 0.01 = **p < 0.001 = ***p < 0.0001 = ****p >0.05 = ns).

Results

DE-NRGs Identification in OA Patients

The four OA-related data collected in the GEO database were combined and batch effects were corrected to minimize errors in the subsequent analysis. P-adjust <0.05 and |logFC|>0.585 were used as filtering criteria for the difference analysis of Con and OA groups, and finally 133 down-regulated differentially expressed genes and 120 up-regulated differentially expressed genes (DEGs) were obtained (Figure 1A and Supplementary Table S3). As shown in the heatmap, the top 50 genes whose expression has been significantly up- or down-regulated are highlighted (Figure 1B). Seven DE-NRGs (TLR7, CRISPLD2, SLC25A37, IL1B, IL6, MMP9 and TLR8) implicated in the OA process were identified by intersecting 255 DEGs with 69 NRGs (Figure 1C). The GSVA algorithm analysis of the NETs lineage related pathways (Neutrophils pathway, Neutrophil degranulation and Neutrophil granule constituents.) showed that the pathway activity was significantly elevated in the OA group to the Con group (Figure 1D).

Figure 1 Identification of OA-related DEGs. (A and B) Volcano and heat maps show differentially expressed genes in OA patients. (C) Venn diagram presenting DE-NRGs in OA patients. (D) Heat map of changes in the activity of NETs-related pathways in the Con and OA groups.

Differences in the Immune Characteristics of OA Patients

The CIBERSORT algorithm was used to explore the detail of immune cell infiltration abundance between the OA and Con groups. OA patients had significantly higher abundance of T-cells CD4 memory activated, Macrophages M0 and Mast cells resting than in the Con group. In contrast, Plasma cells, T cells CD4 memory resting, Monocytes, Mast cells activated had a higher percentage of infiltration in the Con group (Figure 2A and B).

Figure 2 Comparison of the immune infiltration landscape between the Con and OA groups. (A) Stacked histogram of changes in the proportion of immune cells. (B) Infiltration of 22 immune cells in the Con group and OA group. p < 0.05= *, p < 0.01= **, p < 0.001= ***, p >0.05= ns.

Identification of Core DE-NRGs and Analysis of Interactions with Trait-Associated Genes

We integrated the results of three machine learning algorithms, LASSO, SVM- RFE, and RF, and determined the most critical DE-NRGs for the OA process. Based on the results of linear regression, LASSO establishes five important DE-NRGs (Figure 3A and B). The results of SVM- RFE and RF indicated that these five genes were ranked as the highest (Figure 3C–F and Supplementary Table S4). Using correlation analysis of these five DE-NRGs we found that TLR7 and MMP9 were significantly positively correlated, while SLC25A37 CRISPLD2 and IL1B were significantly negatively correlated (Figure 4A). Next, we constructed a PPI co-expression network associated with DE-NRG by the GeneMANIA tool and performed KEGG/GO enrichment analysis (Figure 4B–D). The KEGG results focused on MAPK signaling pathway and cytokine-cytokine receptor interactions. Biological pathways, molecular functional outcomes involving cytokine-mediated signaling pathways and binding of cytokine receptors. The results of cellular components suggest that collagen-containing extracellular matrix and secretory granular lumen may be involved.

Figure 3 Machine learning algorithms screen key DE-NRGs. (A and B) Obtaining 5 candidate DE-NRGs through LASSO regression with 10-fold cross-validation. (C and D) Characteristic DE-NRGs (n=5) were established after identification by the SVM-RFE algorithm. (E) The plot depicts a random forest tree where the x-axis represents the number of trees and the y-axis represents the error rate. The dots in red, green, and black represent OA samples, non-OA samples, and all samples, respectively. (F) The horizontal axis represents gene importance and the vertical axis represents DE-NRGs.

Figure 4 DE-NRGs interaction relationship and correlation analysis. (A) Correlation between DE-NRGs. (B) The co-expression network of DE-NRGs. (C) KEGG analysis of co–expressed genes. (D) GO analysis co–expressed genes.

NETs-Related Nomogram Mode Construction and Validation

We developed a NETs related Nomogram model related based on 5 DE-NRGs (TLR7, CRISPLD2, SLC25A37, IL1B, and MMP9) to provide physicians with a streamlined and reliable diagnostic tool (Figure 5A). The calibration curve analysis showed that the accuracy of the Nomogram model for DE-NRGs was very similar to the actual positivity rate (Figure 5B). The results of both the decision curve analysis and the clinical impact curve analysis demonstrated that our Nomogram model constructed based on DE-NRGs can offer a significant level of support for identifying OA, as evidenced in Figures 5C and D. The ROC curve in Figures 5E–I also displays that the classification performance of CRISPLD2, SLC25A37, IL1B, MMP9, and TLR7 is excellent, as indicated by the high Area Under the Curve (AUC) values. Finally, we also verified the expression of CRISPLD2, SLC25A37, IL1B, MMP9, and TLR7 in GSE55235 and showed that these 5 DE-NRGs were also significantly different (Figure S1).

Figure 5 Construction of DE-NRGs related Nomogram to assess clinical value. (A) Nomogram demonstrates the prognostic value of CRISPLD2, IL1B, SLC25A37, MMP9, and TLR7 for OA patients. (B) Calibration curves to assess the degree of similarity between the predicted and true results of DE-NRGs related Nomogram. (C) Decision Curve Analysis to evaluate the sensitivity and specificity of the DE-NRGs related Nomogram. (D) Clinical impact curve to assess the clinical impact of the DE-NRGs related Nomogram at different thresholds. (E–I) ROC analysis results on CRISPLD2, IL1B, SLC25A37, MMP9 and TLR7.

Predicting DE-NRGs Targeted Drugs

We delved deeper into the potential drugs that can effectively target DE-NRGs through analysis of the DGIdb database and visualized the results with Cytoscape software (Figure 6). Regrettably, our analysis did not yield a prediction for a drug targeting CRISPLD2. We ultimately uncovered 34 drugs that could potentially impact DE-NRGs, among which were 4 drugs targeting TLR7, 1 drug targeting CRISPLD2, 28 drugs targeting IL1B, and 1 drug targeting MMP9. The findings suggest that HYDROXYCHLOROQUINE SULFATE can be used as an antagonist of TLR7, while IMIQUIMOD can be utilized as an agonist of TLR7. RILONACEPT and CANAKINUMAB can be utilized as inhibitors of IL1B. Specific information can be found in Supplementary Table S5.

Figure 6 Prediction of DE-NRGs targeted drugs. Network diagram showing that the DGIdb database predicts DE-NRGs related drugs.

ssGSEA and GSEA to Explore the Formation of NETs

To determine whether there are significant differences in the activities of Neutrophils pathway, Neutrophil degranulation and Neutrophil granule constituents pathway, we performed validation in the GSE55235 dataset using the ssGSEA algorithm. The activity of Neutrophils pathway, Neutrophil degranulation and Neutrophil granule constituents pathway in OA group was significantly higher than that in Con group (Figure 7A). Results indicated a significant positive correlation between TLR7 and the Neutrophil degranulation and Neutrophil granule constituents pathway. Additionally, upregulation of MMP9 may boost the activity of Neutrophils pathway, Neutrophil degranulation and Neutrophil granule constituents pathway. CRISPLD2 and IL1B showed a significant negative correlation with the Neutrophil degranulation and Neutrophil granule constituents pathway. SLC25A37 was only negatively correlated with the Neutrophil granule constituents pathway (Figure 7B). Subsequently, we performed gene enrichment analysis on TLR7 in the GSE55235 and GSE55457 datasets. The results showed that TLR7 was enriched in the Neutrophil extracellular trap formation pathway in both GSE55235 and GSE55457 datasets (Figure 7C and D).

Figure 7 Correlation between OA and NETs. (A) Comparison of ssGSEA scores of Neutrophil extracellular trap-related pathways in OA and Con groups. (B) Correlation between DE-NRGs and Neutrophil extracellular trap-related pathways. (C and D) GESA results of TLR7 at GSE55235 and GSE55457. p < 0.05= *, p < 0.01= **, p < 0.001= ***.

Validation of Bioinformatics Results by qPCR and Immunofluorescence

To elucidate the contribution of TLR7, CRISPLD2, SLC25A37, IL1B, and MMP9 in OA, we started experiments to verify the authenticity of human OA samples and established mouse OA models by staining with Safranin O-Fast Green (Figure 8A–C). This step was essential to ensure the reliability of the samples used in our investigation. Subsequently, the mouse OA model underwent medial tibial plateau (MTP) and medial femoral condyle (MFC) scoring to further reinforce the validity of the model (Figure 8D and E). The qPCR results showed that TLR7, CRISPLD2, IL1B, and MMP9 were significantly different in mRNA expression levels and were consistent with the bioinformatics analysis (Figure 8F and G). However, the difference in SLC25A37 mRNA expression was only observed in the mouse OA model (Figure 8F and G). As both the SVM-RFE and RF algorithms demonstrated that TLR7 was the most critical factor, we subsequently focused on analyzing TLR7. In the ATDC5 cell experiment, TNF-α stimulation led to a substantial increase in mRNA expression levels of TLR7, MMP9, and IL6, while CRISPLD2, IL1B, and SLC25A37 were notably downregulated (Figure 8H). Finally, we found that TLR7 expression was elevated in both human OA samples and mouse OA samples by immunofluorescence assay (Figure 8I–L).

Figure 8 Experimental validation of TLR7. (A) Representative images of Safranin O–fast green–stained sections of articular cartilage from human with or without OA (n=6). (B) OARSI scoring system evaluation of human articular cartilage. (C) Representative images of Safranin O–fast green–stained sections of the keen joint from mice after sham operation or DMM surgery (n=5). (D and E) Severity of articular cartilage damage in mice after sham operation or DMM surgery, evaluated using the OARSI scoring system. MFC: medial femoral condyle; MTP: medial tibial plateau. (F and G) Changes of mRNA expression levels of CRISPLD2, IL1B, SLC25A37, MMP9, and TLR7 in human and mice samples (n=3). (H) After TNFα-induced, the mRNA expression levels of Tlr7, Mmp9, Il6, CRISPLD2, IL1B, and SLC25A37 were significantly higher in the ATDC5 cell line (n=3). (I and K) Immunofluorescence staining for TLR7 was conducted on articular cartilage samples obtained from both healthy individuals and OA patients, followed by quantitative fluorescence analysis (n=3). TLR7 stains green. Dapi stains the nucleus. (J and L) Immunofluorescence staining of Tlr7 was performed in articular cartilage from mice at sham or DMM-surgery, and quantitative fluorescence statistics were analyzed. Tlr7 stains green. Dapi stains the nucleus (n=3). p < 0.05= *, p < 0.01= **, p < 0.001= ***, p < 0.0001= ****, p >0.05 = ns.

Discussion

OA is a widespread degenerative joint condition that results in the deterioration of articular cartilage and bone. The progression of OA is intricately associated with aging, joint injury, and muscle wasting around the joints, alongside inflammation of the synovial membrane, degeneration of the infrapatellar fat pad, and meniscal degeneration.²⁷ Approximately 60% of the elderly population is estimated to be affected by OA, causing significant degradation in their quality of life. The OA process has been suggested to be impacted by various cellular death mechanisms such as apoptosis, necroptosis, and cuproptosis,²⁸ but the contribution of NETs to OA remains unclear.

NETs are composed of decondensed chromatin fibers and granule proteins such as histones and neutrophil elastase, which are released by activated neutrophils during degranulation.²⁹ NETs can trap and kill pathogens, but they can also contribute to tissue damage in certain pathologies, such as autoimmune diseases and cancer. Neutrophils, which contain specialized organelles called neutrophil granules.³⁰ These granules harbor a diverse range of molecules, including enzymes, antimicrobial peptides, and various proteins, enabling them to play a significant role in the body’s immune response. Upon encountering a pathogen, neutrophils activate the release of the granules they contain, thereby discharging their contents into the surrounding area and effectively combating the pathogen.

At the outset, we conducted a comprehensive analysis by merging four OA-related datasets, effectively correcting for batch effects, and subsequently identifying 253 DEGs with regulatory relevance, while comparing the control and OA groups. Our analysis resulted in the identification of 7 DE-NRGs including TLR7, SLC25A37, IL11B, CRISPLD2, MMP9, TLR8, and IL6. GSVA evaluates the effect of specific genomic alterations on disease by assigning scores based on the activity of gene sets involved in particular biological processes or functions.³¹ Using GSVA, we found that the gene sets for Neutrophils pathway, Neutrophil degranulation and Neutrophil granule constituents pathway were significantly more activated in OA group compared to the control group. Next, we used the CIBERSORT algorithm to analyze immune cell infiltration in OA patients. We discovered that plasma cells, monocytes, and activated mast cells were significantly more infiltrated in the control group compared to the OA group, while T cells CD4 memory resting, T cells CD4 memory activated, Macrophages M0, and resting mast cells were significantly more infiltrated in the OA group compared to the control group.

Machine learning algorithms offer various advantages, including improved model accuracy, reduced overfitting, enhanced interpretability, and effective handling of high-dimensional data.³² In the field of biomedical research, it is possible to more accurately predict core molecules in disease processes and identify potential clinical features.

The elevated expression of TLR7 in the blood of RA patients makes the synovial tissue more susceptible to reacting to citrullinated proteins and nucleic acids that are released into the extracellular environment through the formation of NETs.³³ The research findings demonstrate a sustained down-regulation of SCL25A37 in individuals with Major Depressive Disorder (MDD), suggesting its potential utility as a biomarker for MDD diagnosis and a promising target for future therapeutic and diagnostic interventions.³⁴ CRISPLD2 is a gene that is differentially expressed in both polycystic ovary syndrome and ovarian cancer and has significant clinical implications for the prognosis of ovarian cancer.³⁵ The levels of serum MMP9 expression in OA patients increase significantly after total joint arthroplasty treatment, reaching a stable level by the fifth day. Our findings align with this observation, indicating the potential utility of detecting MMP9 expression in both blood and joint samples. Such assessments could offer valuable guidance to surgeons, enabling them to refine treatment strategies and mitigate surgical risks, ultimately leading to improved patient outcomes.³⁶ Based on the expression of five DE-NGRs, we created NETs-related nomograms for the purpose of facilitating clinical risk assessment and treatment planning in OA patients. The follow-up five DE-NGR highly related genes were subjected to KEGG and GO analysis.

The KEGG results were most closely related to interactions between cytokine receptors and cytokines, osteoclast differentiation, and MAPK signaling. The biological processes, cellular components, and molecular functional analysis results are related to cytokine-mediated signaling pathways, the extracellular matrix containing collagen, and binding of cytokine receptors, respectively. During the progression of OA, cytokine-cytokine receptor interactions may lead to activation of the MAPK signaling pathway, which then promotes the differentiation of osteoclasts. This increased osteoclast activity can promote destruction of articular cartilage and bone, leading to the characteristic symptoms of OA, such as pain, stiffness, and decreased mobility.³⁷ Additionally, the MAPK signaling pathway plays a role in regulating cytokine-mediated signaling pathways, including those that drive inflammation by producing cytokines.³⁸ Because both SVM-RFE, and RF analyses showed that TLR7 was the top ranked DE-NRGs, we performed a single gene GSEA on TLR7. In GSE55235 and GSE554757, TLR7 was significantly enriched in the Neutrophil extracellular trap formation. Not only that, TLR7 was also significantly correlated with Neutrophil degranulation and Neutrophil granule constituents progression. Thus, we speculate that TLR7 may influence the OA process by activating the Neutrophil extracellular trap formation pathway and promoting the formation of NETs. We confirmed the effect of TLR7 on OA by constructing mouse OA models, human OA samples, and cellular experiments. The results consistently showed that TLR7 expression was significantly increased in OA.

However, our study does have certain limitations. Firstly, the relatively small sample size of the dataset may introduce some bias in our experimental results. To address this, future studies with larger sample sizes and advancements in sequencing technology will allow us to build more accurate prognostic models. Secondly, we acknowledge the need for further investigations into the specific mechanisms through which TLR7 influences the formation of NETs. In conclusion, our research offers valuable insights into the diagnosis and treatment of OA through targeted NETs, while recognizing the potential for future research to expand and deepen our understanding in this field.

Conclusions

In summary, our study comprehensively explored the association between the NETs process and OA and finally identified the hub genes. Next, the hub genes were combined with immune cells, NETs formation and gene-targeted drugs for an integrated analysis. Finally, TLR7 was found to potentially induce the formation of NETs to drive the progression of OA.