The expression profile datasets of IDD tissue samples, GSE34095, GSE124272, and GSE147383, were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The removeBatchEffect function from the limma package was used to merge these three datasets and eliminate batch effects. The merged dataset included 15 disease tissue samples and 15 normal control tissue samples. Additionally, the expression profile datasets GSE150408 (17 normal and disease samples) and GSE70362 (16 normal and 32 disease samples) were also downloaded as external validation.

The limma package was used to analyze the differential gene expression between the disease group and the normal group in both the integrated dataset and GSE150408. The screening criteria were set as |log2FC| > 0.5 and p < 0.05. Based on previously published literature, a set of 4,463 PDGs was identified (Finan Kruger et al., 2017). The VennDiagram package was used to generate a Venn diagram to identify druggable genes that were differentially expressed in both datasets.

To explore the biological functions of the identified druggable genes, GO and KEGG pathway enrichment analyses were conducted using the clusterProfiler package. The significance threshold was set at p < 0.05. The results were visualized using the ggplot2 package.

The CIBERSORT algorithm was employed to predict the distribution of immune cells in the samples, and a t-test was conducted to compare the proportions of immune cells between the disease and normal groups, aiming to investigate immune cell infiltration changes between the two groups. Furthermore, to ensure the reliability of the results, the MCPcounter and TIMER algorithms were used for additional validation.

To further explore the diagnostic role of the differentially expressed druggable genes in IDD, Lasso regression analysis was performed on the integrated dataset using the glmnet package. A 5-fold cross-validation was conducted, with the lambda parameter set to 0.05784, to obtain the optimal binary classification diagnostic model. Additionally, to prevent overfitting, the classification performance of the model was also validated in GSE150408. The pROC package was used to generate ROC curves to quantify the predictive performance of the model.

The GWAS data for the exposure factors (diagnostic model characteristic genes) and outcome variables were obtained from the IEU Open GWAS Project (https://gwas.mrcieu.ac.uk/), and the details are shown in Supplementary Table S1. The extract_instruments function in the TwoSampleMR package was used to select SNPs that were strongly associated with the exposure factors from the GWAS data, using a threshold of P-value < 5e-08. The clump_data function was then applied to remove SNPs with linkage disequilibrium to ensure the independence of instrumental variables, with the parameters set as clump_r2 = 0.001 and clump_kb = 10,000. MR analysis was performed using the inverse variance weighted (IVW) method, weighted median (WM) estimation, and MR-Egger regression, implemented in the TwoSampleMR package. Additionally, IVW and MR-Egger tests were conducted to assess the presence of heterogeneity and pleiotropy in the results.

Human NP cells (NPCs) were purchased from ScienCell (CA, USA) and cultured in Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin in T25 culture flasks at 37 °C in a humidified incubator with 5% CO₂. When NPCs reached 90% confluence, they were detached using 0.25% trypsin/1 mM EDTA and subcultured into larger flasks for expansion. Third- and fourth-generation NPCs were used for subsequent experiments.

For TNF-α treatment, NPCs were seeded at a density of 5,000 cells per well in 96-well plates with 200 μL of complete culture medium and incubated overnight at 37 °C with 5% CO₂. The next day, the medium was removed, and cells were washed with PBS. Fresh medium containing TNF-α (100 ng/mL) was then added for 24 h (Duan Kuang et al., 2023). Cells were incubated under these conditions for the indicated time points before further analysis.

Briefly, NPCs were seeded in 96-well plates at a density of 5,000 cells per well and incubated at 37 °C with 5% CO₂ overnight. After the designated treatments, 10 μL of CCK-8 reagent was added to each well, followed by incubation at 37 °C for 2 h. The absorbance at 450 nm was measured using a microplate reader.

Total RNA was extracted from NPCs using TRIzol reagent (Invitrogen, CA, USA). The RNA concentration and purity were assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific, CA, USA). For cDNA synthesis, 1 μg of total RNA was reverse-transcribed using the PrimeScript RT reagent Kit following the manufacturer’s protocol. The reaction was carried out at 37 °C for 15 min, followed by 85 °C for 5 s to inactivate the reverse transcriptase. RT-qPCR was performed using TB Green™ Premix Ex Taq™ on a QuantStudio™ 5 Real-Time PCR System. The thermal cycling conditions were as follows: 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s. Gene expression levels were normalized to GAPDH as an internal control, and relative expression was calculated using the 2^−ΔΔCT method. Each experiment was performed in triplicate to ensure reproducibility.

NPCs were harvested and washed twice with cold PBS, then resuspended in 100 μL of binding buffer. Apoptosis was assessed using the Annexin V-FITC/PI Apoptosis Detection Kit following the manufacturer’s protocol. Briefly, cells were stained with 5 μL of Annexin V-FITC and 5 μL of propidium iodide (PI) and incubated in the dark at room temperature for 15 min. After staining, 400 μL of binding buffer was added to each sample, and apoptosis was analyzed using a flow cytometer within 1 h. The proportion of apoptotic cells was determined using FlowJo software.

NPCs were lysed using RIPA buffer on ice for 30 min. The lysates were centrifuged at 12,000 × g for 15 min at 4 °C, and the supernatants were collected. Protein concentration was determined using the BCA protein assay kit. Equal amounts of protein were mixed with 5× loading buffer, boiled at 95 °C for 5 min, and separated on a 10% or 12% SDS-PAGE gel. A pre-stained protein molecular weight marker (#G2086-250UL, Servicebio, Hubei, China) was loaded in parallel. The proteins were then transferred onto a PVDF membrane using a wet transfer system at 300 mA for 90 min. After blocking with 5% non-fat milk or 5% BSA in TBST for 1 h at room temperature, the membranes were incubated overnight at 4 °C with primary antibodies against BPI (1:1,000, #ab175231, Abcam, UK), CTSG (1:1,000, #ab282105, Abcam, UK), MMP-13 (1:1,000, #ab39012, Abcam, UK), Collagen II (1:1,000, #ab188570, Abcam, UK), and GAPDH (1:10,000, #ab128915, Abcam, UK) as a loading control. The next day, membranes were washed with TBST and incubated with HRP-conjugated secondary antibodies (1:2,000, #ab288151, Abcam, UK) for 1 h at room temperature. Protein bands were visualized using an ECL detection reagent. The relative expression levels of BPI and CTSG were quantified using ImageJ software.

The concentrations of inflammatory cytokines IL-6 and COX-2 in the culture supernatants of NPCs were determined using ELISA kits according to the manufacturer’s instructions. Briefly, after siRNA transfection and TNF-α treatment, cell culture media were collected and centrifuged at 1,000 × g for 10 min at 4 °C to remove debris. The supernatants were then subjected to ELISA analysis using human IL-6 and COX-2 ELISA kits. Optical density values were measured at 450 nm using a microplate reader, and cytokine concentrations were calculated based on standard curves generated from serial dilutions of known standards.

All statistical analyses were performed using R software (version 4.3.2) and GraphPad Prism (version 9.0). Data are presented as the mean ± standard deviation. Differences between two groups were assessed using the unpaired Student’s t-test. Pearson correlation analysis was used to assess relationships between variables. Statistical significance was defined as p < 0.05.

First, three GEO datasets (GSE34095, GSE124272, and GSE147383) were merged as the training set. After removing batch effects, a uniform sample distribution was observed. A total of 408 DEGs were identified from the merged dataset, including 174 downregulated and 234 upregulated genes (Figure 1A). In the validation set GSE150408, 277 DEGs were obtained (102 downregulated and 175 upregulated, Figure 1B). The intersection of these DEGs with 4,463 PDGs yielded 14 overlapping genes (Figure 1C), identified as differentially expressed PDGs in IDD. The detailed gene names are listed in Table 1.

To investigate the functional roles of these 14 PDGs, GO and KEGG enrichment analyses were performed (Figure 2). GO analysis revealed that biological processes (BP) were primarily enriched in defense responses to bacteria (especially Gram-negative bacteria) and fungi, cell killing, and leukocyte-mediated cytotoxicity (Figure 2A). The top 5 cellular components (CC) included primary lysosome, azurophil granule, secretory granule lumen, and cytoplasmic vesicle lumen (Figure 2A). Molecular functions (MF) were enriched in heparin binding, glycosaminoglycan binding, sulfur compound binding, heme binding, and tetrapyrrole binding (Figure 2A). Key KEGG pathways included neutrophil extracellular trap formation, acute myeloid leukemia, and Staphylococcus aureus infection (Figure 2B).

Given the association of these 14 genes with antimicrobial immune responses, immune infiltration analysis was conducted. The CIBERSORT algorithm revealed that only neutrophil infiltration levels showed significant differences between IDD and normal samples, with an upregulation observed in IDD (p < 0.01, Figure 3A). To verify the stability of these findings, immune infiltration analysis was further conducted using two additional algorithms. The results consistently demonstrated significantly higher neutrophil levels in IDD samples compared to normal samples (p < 0.05, Figure 3B). A correlation heatmap was used to explore the relationships between the 14 PDGs and immune cells. Notably, neutrophils exhibited significant negative correlations with KIF11, FCGBP, and CD160, while showing positive correlations with FCGR1A, CYP4F2, CTSG, and BPI (Figure 3C).

To explore the clinical diagnostic value of these 14 PDGs in IDD, we constructed a 5-gene diagnostic model using Lasso regression (Figure 4A), comprising BPI, CD160, CTSG, CYP27A1, and KIF11. This model demonstrated strong diagnostic performance in the training, testing, and validation sets, with AUC values of 0.969, 0.889, and 0.869, respectively (Figure 4B). Next, we examined the expression levels of the five model genes in normal and IDD samples. The results showed that BPI, CTSG, CYP27A1, and KIF11 were significantly upregulated in IDD, while CD160 was significantly downregulated (p < 0.05, Figure 4C).

Subsequently, we further explored the genetic association between the five model genes and IDD through MR analysis. Among the five MR methods applied, only BPI and CTSG exhibited significant causal associations with IDD (Table 2). SNP effect plots indicated a negative directional association between both BPI and CTSG and IDD (Figures 5A,C). Cochran’s Q test and MR-Egger intercept suggested that the heterogeneity and pleiotropy in the IDD-BPI and IDD-CTSG associations were negligible (all p > 0.05, Tables 3, 4). Additionally, leave-one-out analysis confirmed the robustness of the findings, as all effect estimates remained on the same side (Figures 5B,D).

To explore the potential mechanisms of BPI and CTSG in IDD, we performed GSEA. The results indicated that BPI was significantly negatively correlated with oxidative phosphorylation (Figure 6A), while CTSG was negatively correlated with the G2M checkpoint (Figure 6B).

Finally, we validated the expression of the model genes through cellular experiments. TNF-α-induced NPCs were used to simulate IDD-related damage. Compared to the control group, cell viability was significantly reduced in the TNF-α group (p < 0.001, Figure 7A), while the apoptosis rate was markedly increased (p < 0.01, Figure 7B), indicating significant cellular damage in the TNF-α-treated NPCs. Consistent with the bioinformatics analysis, the damaged NPCs exhibited upregulated expression of BPI, CTSG, CYP27A1, and KIF11, while CD160 expression was downregulated (p < 0.05, Figure 7C). Western blotting was performed to further validate the protein expression of the two core genes, BPI and CTSG. As shown in Figure 7D, the protein levels of BPI and CTSG were significantly upregulated in the damaged NPCs compared to the control group (p < 0.001, Figure 7D).

Furthermore, to explore the biological functions of BPI and CTSG in NPCs, gene silencing experiments were performed using siRNA. The knockdown efficiency was confirmed by qRT-PCR and Western blotting, which demonstrated that both mRNA and protein levels of BPI and CTSG were markedly reduced after siRNA transfection (p < 0.01, Figures 8A,B). Flow cytometric analysis revealed that the apoptotic rate of NPCs was significantly decreased following the knockdown of BPI or CTSG compared with the si-NC group (p < 0.001, Figure 8C), suggesting that silencing these genes could alleviate TNF-α-induced apoptosis. Moreover, the secretion of pro-inflammatory cytokines IL-6 and COX-2 was notably decreased in the si-BPI and si-CTSG groups (p < 0.01, Figure 8D), indicating a reduction in inflammatory response. To further investigate the effects of BPI and CTSG silencing on extracellular matrix metabolism, we assessed the expression of MMP-13 and Collagen II. qRT-PCR and Western blotting analyses revealed that inhibition of BPI or CTSG significantly suppressed MMP-13 expression while enhancing Collagen II levels (p < 0.001, Figures 8E,F). These findings indicate that silencing BPI and CTSG mitigates inflammation, inhibits matrix degradation, and promotes ECM synthesis in NPCs, thereby exerting a protective effect against IDD.