Multiomics Integration Analysis Reveals the Regulatory Mechanisms of Efferocytosis in Diabetic Kidney Disease.

The GSE96804 and GSE30122 datasets were merged to remove batch effects (Fig.  3 A-B). A total of 378 DEGs were identified in the DKD group compared with the normal group. The DKD group presented 141 upregulated genes and 237 downregulated genes. A histogram of the adjusted P value distribution for the DEGs is shown in Figure S1 . A heatmap (Fig.  3 C) and volcano plot (Fig.  3 D) were generated to visualize the expression differences and clustering features of the DEGs between the normal and DKD groups. Fig. 3 Identification of DEGs ( A ) UMAP distribution before batch removal; ( B ) UMAP distribution after batch removal; ( C ) Heatmap; ( D ) Volcano plot Identification of DEGs ( A ) UMAP distribution before batch removal; ( B ) UMAP distribution after batch removal; ( C ) Heatmap; ( D ) Volcano plot The intersection of DEGs and 256 ERGs yielded 15 overlapping genes (DKD-ERGs), which were then subjected to enrichment analysis (Fig.  4 A). GO enrichment analysis revealed that biological processes (BP) were predominantly involved in the regulation of inflammatory responses and platelet activation. The cellular components (CC) terms were enriched primarily in the collagen fibril extracellular matrix, cytoplasmic vesicle lumen, and related components. The enriched molecular functions (MFs) were glycosaminoglycan binding, calcium-dependent protein binding, and similar functions (Fig.  4 B). KEGG pathway analysis revealed significant enrichment of efferocytosis, complement and coagulation cascades, apoptosis, IL-17 signaling, and ECM-receptor interactions (Fig.  4 C). Figure  4 D shows the chromosomal regions of these genes, facilitating the identification and analysis of similarities and differences in comparative genomic studies. Figure  4 E shows the PPI networks, demonstrating how the proteins interact. Genes were ranked according to their degree, with a higher degree value corresponding to a darker node color, indicating greater importance of the gene represented by the node. The top 10 genes ranked by degree were selected as key targets, and these 10 targets were ANXA1, S100A9, C3, CASP3, CLU, FN1, APOE, IL33, PLG, and CX3CR1 (Fig.  4 F). Fig. 4 Enrichment analysis and PPI network. ( A ) Venn diagram; ( B ) GO enrichment analysis; ( C ) KEGG enrichment analysis; ( D ) chromosomal locations; ( E ) PPI network; ( F ) top 10 hub genes based on degree. (A higher degree corresponds to a darker color, indicating that the node is more important) Enrichment analysis and PPI network. ( A ) Venn diagram; ( B ) GO enrichment analysis; ( C ) KEGG enrichment analysis; ( D ) chromosomal locations; ( E ) PPI network; ( F ) top 10 hub genes based on degree. (A higher degree corresponds to a darker color, indicating that the node is more important) Next, feature genes were further selected via machine learning using the expression matrix of the 15 genes. Figure  5 A shows the shrinkage trajectories of the regression coefficients (coef) of each gene as the penalty coefficient λ changes; the vertical dashed line indicates the λ value at which the 10-fold cross-validation error is minimized. Figure  5 B presents the six genes (ANXA1, C3, CASP3, IL33, GPR18, and S100A9) corresponding to nonzero coefficients at this λ value. Figure  5 C shows the boxplot of the “shadow variables” from BORUTA: the green boxes represent the confirmed important genes, and the red boxes represent the discarded variables. Figure  5 D summarizes the final 14 genes retained. Figure  5 E presents the mean decrease accuracy ranking and cumulative importance curve calculated by the RF. When the curve plateaus and the cumulative importance reaches 0.85, a total of 12 genes are included. By integrating the results from the LASSO, RF, BORUTA, and degree analyses, four hub genes—ANXA1, CASP3, IL33, and C3—were ultimately identified (Fig.  5 F). These genes are suggested to play pivotal roles in DKD. Fig. 5 Machine learning used to screen for hub genes. ( A ) LASSO regression coefficient path diagram; ( B ) LASSO regression cross-validation curve; ( C - D ) BORUTA algorithm; ( E ) RF algorithm; ( F ) Venn diagram of the degree of overlap among three machine learning methods Machine learning used to screen for hub genes. ( A ) LASSO regression coefficient path diagram; ( B ) LASSO regression cross-validation curve; ( C - D ) BORUTA algorithm; ( E ) RF algorithm; ( F ) Venn diagram of the degree of overlap among three machine learning methods In the training dataset, hub genes were upregulated in the DKD group ( P  < 0.05) (Fig.  6 A). The AUC of the hub genes ranged from 0.74 to 0.85, indicating a certain diagnostic value (Fig.  6 B). As DKD typically culminates in renal fibrosis as a late-stage pathological outcome, we examined the relationships among hub genes and fibrosis markers (FN1, COL1A2, COL1A1, and ACTA2). In addition to ANXA1, CASP3, IL33, and C3 were significantly positively correlated with fibrosis-related genes (Fig. S2 ). In the GSE104948 dataset, as well as in the Ju CKD Glom and Ju CKD TubInt validation datasets, hub genes were significantly upregulated ( P  < 0.05), which is consistent with previous results (Fig.  6 C-E, G-H). In terms of clinical indicators, in the Ju CKD Glom dataset, CASP3, IL33, and C3 were negatively associated with the GFR and positively correlated with Scr (Fig.  6 F). In the Ju CKD TubInt dataset, the hub genes were negatively correlated with the GFR and positively correlated with Scr (Fig.  6 I). Fig. 6 Validation and clinical significance. ( A ) Hub gene expression in the training dataset; ( B ) ROC curves in the training dataset; ( C ) Hub gene expression in the GSE104948 dataset; ( D ) Hub gene expression in the Ju CKD Glom dataset; ( E ) ROC curves in the Ju CKD Glom dataset; ( F ) Correlations between hub genes and clinical parameters in the Ju CKD Glom dataset; ( G ) Hub gene expression in the Ju CKD TubInt dataset; ( H ) ROC curves in the Ju CKD TubInt dataset; ( I ) Correlations between hub genes and clinical parameters in the Ju CKD TubInt dataset. (The error bars represent the means ± standard deviations (means ± SDs, * P  < 0.05, ** P  < 0.01, *** P  < 0.001, **** P  < 0.0001) Validation and clinical significance. ( A ) Hub gene expression in the training dataset; ( B ) ROC curves in the training dataset; ( C ) Hub gene expression in the GSE104948 dataset; ( D ) Hub gene expression in the Ju CKD Glom dataset; ( E ) ROC curves in the Ju CKD Glom dataset; ( F ) Correlations between hub genes and clinical parameters in the Ju CKD Glom dataset; ( G ) Hub gene expression in the Ju CKD TubInt dataset; ( H ) ROC curves in the Ju CKD TubInt dataset; ( I ) Correlations between hub genes and clinical parameters in the Ju CKD TubInt dataset. (The error bars represent the means ± standard deviations (means ± SDs, * P  < 0.05, ** P  < 0.01, *** P  < 0.001, **** P  < 0.0001) DKD animal models were constructed to evaluate the hub genes. One week after STZ injection, the FPG level in the DKD group exceeded 16.7 mmol/L, indicating effective diabetes induction (Fig. 7 A). Compared with those in the control group, the DKD group fed a HFD for 20 weeks presented significantly elevated FPG, kidney weight-to-body weight ratio (kW/BW), BUN, Scr, and milligrams of albumin per gram of urinary creatinine (mAlb/uCr) levels (Fig. 7 B-E) [ 32 ]. Figure 7 F compares the histopathological characteristics of renal tissues in the control and DKD groups using HE staining, PAS staining, and Masson staining. The control group presented clear glomerular and tubular structures without pathological changes, whereas the DKD group presented significantly enlarged glomeruli, thickened tubular basement membranes, and marked interstitial fibrosis according to Masson’s trichrome staining, confirming successful model establishment. Fig. 7 Evaluation of the DKD mouse model. ( A ) Blood glucose changes from 0 to 20 weeks; ( B ) Kidney weight/body weight ratio at 20 weeks; ( C ) BUN at 20 weeks; ( D ) Scr at 20 weeks; ( E ) mA/uCr at 20 weeks; ( F ) HE, PAS, and Masson staining (50 μm). The error bars represent the means ± standard deviations (means ± SDs, n  = 6, * P  < 0.05, ** P  < 0.01). Evaluation of the DKD mouse model. ( A ) Blood glucose changes from 0 to 20 weeks; ( B ) Kidney weight/body weight ratio at 20 weeks; ( C ) BUN at 20 weeks; ( D ) Scr at 20 weeks; ( E ) mA/uCr at 20 weeks; ( F ) HE, PAS, and Masson staining (50 μm). The error bars represent the means ± standard deviations (means ± SDs, n  = 6, * P  < 0.05, ** P  < 0.01). Next, we performed validation by multiomics data from mice. Figure 8 A shows the differences in the transcription of the 4 hub genes. Anxa1, Il33, and C3 were notably upregulated in the DKD group ( P 0.05). Additionally, Anxa1 was correlated with BUN, whereas Il33 and Scr were correlated with each other (Fig. 8 B). Figure 8 C shows the relative protein abundance of Annexin 1, Caspase-3, complement C3, and interleukin-33 in DKD renal tissues. Compared with those in the control group, the levels of all these proteins, except Caspase-3, were significantly elevated in the DKD group. Furthermore, Annexin 1 was correlated with BUN, whereas Caspase-3, Interleukin-33, and Complement C3 were correlated with Scr and mA/uCr (Fig. 8 D), suggesting that these genes could be crucial in the pathophysiology of DKD. Efferocytosis is an important pathway for the generation of specialized proresolving mediators (SPMs) [ 18 , 39 ], and SPMs, in turn, can enhance efferocytosis and promote inflammation resolution [ 40 , 41 ]. On the basis of metabolomics data, we further analyzed the levels of classic SPMs—Resolvin D1, Resolvin D2, Maresin 1, and Resolvin E2—in DKD mice, with the aim of performing a reverse assessment of efferocytosis intensity. The results revealed that the levels of SPMs were significantly lower in the DKD group than in the control group, suggesting an imbalance of SPMs in patients with DKD (Fig. 8 E). These SPMs were correlated with renal function (Fig. 8 F). These findings further support the possibility of a defect in efferocytosis function in DKD. Fig. 8 Multiomics validation of hub genes ( A ) Boxplots of transcriptomics levels; ( B ) Correlations between transcriptomics levels of hub genes and renal function; ( C ) Boxplots of proteomics levels for hub genes; ( D ) Correlations between proteomics levels of hub genes and renal function; ( E ) Boxplots showing SPMs levels in metabolomics; ( F ) Correlations between SPMs and renal function. (The error bars represent the means ± standard deviations (means ± SDs, n  = 3, * P  < 0.05, ** P  < 0.01). Multiomics validation of hub genes ( A ) Boxplots of transcriptomics levels; ( B ) Correlations between transcriptomics levels of hub genes and renal function; ( C ) Boxplots of proteomics levels for hub genes; ( D ) Correlations between proteomics levels of hub genes and renal function; ( E ) Boxplots showing SPMs levels in metabolomics; ( F ) Correlations between SPMs and renal function. (The error bars represent the means ± standard deviations (means ± SDs, n  = 3, * P  < 0.05, ** P  < 0.01). On the basis of the relevant literature, directly detecting circulating IL-33 is not feasible [ 42 ]. Therefore, we focused on the serum levels of ANXA1, Caspase-3 and C3. For the clinical study, I collected serum samples from 22 T2DM patients and 66 DKD patients, and the expression of ANXA1 and Caspase-3 in the serum was measured using ELISA. There were no significant differences in age, sex, or body mass index (BMI) among the four groups ( P > 0.05). Compared with those in the T2DM and DKD III groups, the Scr, BUN, urinary albumin‒creatinine ratio (UACR), and 24-hour urinary protein (24-UTP) levels were significantly elevated, whereas the eGFR was significantly reduced in the DKD IV and DKD V groups ( P < 0.05) (Table S2 ). Compared with those in the T2DM and DKD III groups, serum ANXA1 levels were significantly lower in the DKD IV and DKD V groups (Fig. 9 A). However, no significant difference was observed between the DKD IV and DKD V groups. The serum ANXA1 level in DKD patients was negatively correlated with the Scr level and positively correlated with the eGFR (Fig. 9 B-D). Compared with those in the T2DM group, the serum Caspase-3 levels were significantly increased in the DKD III, DKD IV and DKD V groups. Compared with those in the DKD III group, the serum Caspase-3 levels were significantly increased in the DKD IV and DKD V groups (Fig. 9 E). Serum Caspase-3 in DKD patients was negatively correlated with Scr and 24-UTP and positively correlated with the eGFR (Fig. 9 F-H). Fig. 9 Serum ANXA1, Caspase-3 and C3 levels and correlation with renal function in DKD patients. ( A ) Serum ANXA1 levels; ( B ) correlation between serum ANXA1 and Scr; ( C ) correlation between serum ANXA1 and eGFR; ( D ) correlation between serum ANXA1 and 24-UTP; ( E ) serum Caspase-3 levels; ( F ) correlation between serum Caspase-3 and Scr; ( G ) correlation between serum Caspase-3 and eGFR; ( H ) correlation between serum Caspase-3 and 24-UTP; ( I ) serum C3 levels; ( J ) correlation between serum C3 and Scr; ( K ) correlation between serum C3 and GFR; ( L ) correlation between serum C3 and 24-UTP. (Compared with the T2DM group, * P  < 0.05, **  P< 0.01; compared with the DKD III group, # P  < 0.05, ## P  < 0.01; compared with the DKD IV group, ▲ P  < 0.05, ▲▲ P  < 0.01.) Serum ANXA1, Caspase-3 and C3 levels and correlation with renal function in DKD patients. ( A ) Serum ANXA1 levels; ( B ) correlation between serum ANXA1 and Scr; ( C ) correlation between serum ANXA1 and eGFR; ( D ) correlation between serum ANXA1 and 24-UTP; ( E ) serum Caspase-3 levels; ( F ) correlation between serum Caspase-3 and Scr; ( G ) correlation between serum Caspase-3 and eGFR; ( H ) correlation between serum Caspase-3 and 24-UTP; ( I ) serum C3 levels; ( J ) correlation between serum C3 and Scr; ( K ) correlation between serum C3 and GFR; ( L ) correlation between serum C3 and 24-UTP. (Compared with the T2DM group, * P  < 0.05, **  P< 0.01; compared with the DKD III group, # P  < 0.05, ## P  < 0.01; compared with the DKD IV group, ▲ P  < 0.05, ▲▲ P  < 0.01.) A retrospective study (clinical Study II) included 128 DKD patients for C3 detection, including 29 patients in the DKD III group, 68 patients in the DKD IV group, and 31 patients in the DKD V group. There were no significant differences in age, sex, or BMI among the three groups ( P  > 0.05). Compared with the DKD III group, the DKD IV and DKD V groups had significantly higher Scr and 24-UTP levels, whereas the eGFR was significantly lower. Compared with the DKD IV group, the DKD V group had significantly higher Scr and 24-UTP levels and significantly lower eGFRs (Table S3). Conversely, serum C3 levels were significantly lower in the DKD Ⅴ group than in the DKD Ⅲ and DKD IV groups (Fig.  9 I), indicating that as the disease progresses, renal deposition of complement C3 increases, leading to a decrease in circulating C3 levels. Serum C3 was negatively correlated with Scr ( r = −0.371, P  < 0.001) (Fig.  9 J) and positively correlated with the eGFR ( r  = 0.379, P  < 0.001) ( Fig.  9 K) but was not correlated with the 24UTP (Fig.  9 L ) . Further single-cell analysis was conducted by the KIT database. Figure  10 A displays the t-SNE plot of single-cell sequencing analysis in DKD, showing cell clusters representing 12 distinct cell types, including PODO (podocytes), PCT (proximal convoluted tubules), ENDO (endothelial cells), and PEC (parietal epithelial cells). ANXA1 is mainly localized in PODO, PEC, and ENDO (Fig.  10 B). CASP3 is mainly distributed in CD-ICA (collecting duct-intercalated cell type A) and CD-ICB (collecting duct-intercalated cell type B) (Fig.  10 C). IL33 is predominantly found in ENDO (Fig. 10D). C3 is distributed mainly in the PEC and PCT (Fig.  10 E). Fig. 10 Single-cell sequencing expression profile analysis. ( A ) t-SNE plot of single-cell sequencing distribution; ( B ) bubble plot of ANXA1 expression levels; ( C ) bubble plot of CASP3 expression levels; ( D ) bubble plot of IL33 expression levels; ( E ) bubble plot of C3 expression levels Fig. 11 Construction and functional analysis of the DKD prediction model based on hub genes. ( A ) Nomogram; ( B ) calibration curve; ( C ) ROC curve of risk scores; ( D ) comparison of risk scores between the normal and DKD groups; ( E ) Sankey diagram; ( F ) correlation between the risk score and GFR in Ju CKD Glom; ( G ) correlation between the risk score and the GFR in Ju CKD TubInt; ( H - I ): volcano plot and heatmap of differential genes based on risk score grouping; ( J ) G: KEGG pathway enrichment analysis Single-cell sequencing expression profile analysis. ( A ) t-SNE plot of single-cell sequencing distribution; ( B ) bubble plot of ANXA1 expression levels; ( C ) bubble plot of CASP3 expression levels; ( D ) bubble plot of IL33 expression levels; ( E ) bubble plot of C3 expression levels Construction and functional analysis of the DKD prediction model based on hub genes. ( A ) Nomogram; ( B ) calibration curve; ( C ) ROC curve of risk scores; ( D ) comparison of risk scores between the normal and DKD groups; ( E ) Sankey diagram; ( F ) correlation between the risk score and GFR in Ju CKD Glom; ( G ) correlation between the risk score and the GFR in Ju CKD TubInt; ( H - I ): volcano plot and heatmap of differential genes based on risk score grouping; ( J ) G: KEGG pathway enrichment analysis To better predict the likelihood of patients developing DKD, we developed a nomogram based on four diagnostic genes (Fig.  11 A). The calibration curve indicates that the model has relatively high accuracy (Fig.  11 B). The AUC of the efferocytosis risk score was 0.84 (95% CI: 0.77–0.91), suggesting that the model has relatively good discriminatory power (Fig.  11 C). The risk score was significantly greater in the DKD group (Fig.  11 D). Figure  11 E shows the correspondence between the risk score median grouping and DKD status, with the high-risk group predominantly consisting of DKD patients and the low-risk group mainly comprising normal individuals, further demonstrating the guiding significance of the risk score for clinical stratification. The GSE104948 , Ju CKD TubInt, and Ju CKD Glom datasets were subsequently validated, with DKD patients showing significantly higher risk scores than those in the normal group did, and patients with higher scores presented a more pronounced decline in renal function, indicating a close correlation between the risk score and the degree of kidney injury (Fig.  11 F-G; Fig. S3 A-C). There were 458 DEGs detected between the low-risk and high-risk groups, with 225 downregulated and 233 upregulated (Fig.  11 H). ANXA1, CASP3, IL33, and C3 were elevated in the high-risk group (Fig.  11 H-I ) . Next, GO and KEGG enrichment analyses were conducted on 458 DEGs (Fig.  11 J; Fig. S3 D). KEGG pathway enrichment was observed for the PI3K-Akt signaling pathway, cytokine‒cytokine receptor interaction, phagosome, apoptosis, glycolysis/gluconeogenesis, and NF-kappa B signaling (Fig.  11 J). The results of the consensus clustering analysis indicate that the optimal number of clusters is k = 2 (Fig.  12 A, Fig. S4 A). PCA revealed that Cluster 1 and Cluster 2 corresponded to distinct DKD patients (Fig.  12 B). The ANXA1, CASP3, IL33, and C3 expression levels were higher in Cluster 2 (Fig.  12 C). This finding suggests that Cluster 2 may have more active characteristics in the efferocytosis process. To further elucidate the immune microenvironments of these two clusters, we performed immune infiltration analysis using the CIBERSORT (Fig.  12 D). Correlation analysis revealed that ANXA1, CASP3, IL33, and C3 are closely associated with immune cells, with CASP3 showing a weaker correlation than ANXA1, IL33 and C3 (Fig.  12 E). Further consensus clustering analysis was performed on the Ju CKD Glom and Ju CKD TubInt datasets (Fig. S4 B-C). Compared with Cluster 1 patients, Cluster 2 patients had higher Scr levels and lower GFR (Fig.  12 F-G). Fig. 12 Consensus clustering analysis. ( A ) Consensus clustering matrix (k = 2); ( B ) PCA; ( C ) expression of hub genes in different clusters; ( D ) immune infiltration between two clusters; (E) correlation heatmap between hub genes and immune cells; ( F ) Scr and GFR levels in different clusters of the Ju CKD Glom dataset; ( G ) Scr and GFR levels in different clusters of the Ju CKD TubInt dataset ( * P  < 0.05, ** P  < 0.01, *** P  < 0.001). Consensus clustering analysis. ( A ) Consensus clustering matrix (k = 2); ( B ) PCA; ( C ) expression of hub genes in different clusters; ( D ) immune infiltration between two clusters; (E) correlation heatmap between hub genes and immune cells; ( F ) Scr and GFR levels in different clusters of the Ju CKD Glom dataset; ( G ) Scr and GFR levels in different clusters of the Ju CKD TubInt dataset ( * P  < 0.05, ** P  < 0.01, *** P  < 0.001). The drug‒gene interaction results from the DGIdb database revealed 102 potential target drugs/compounds for DKD treatment. The top 30 drugs/compounds are shown in Fig.  13 , with 4 targeting C3, 4 targeting IL33, 8 targeting ANXA1, and 14 targeting CASP3. Notably, ANXA1 and CASP3 are associated with a greater number of drugs; in contrast, C3 and IL33 have fewer drug associations. Additionally, PEGCETACOPLAN, as a complement C3 inhibitor, has received regulatory approval (shown in green), whereas the majority remain “unapproved,” suggesting that further research is needed for their clinical translation. Fig. 13 Potential DKD therapeutic drug-gene interaction network identified by DGIdb database Potential DKD therapeutic drug-gene interaction network identified by DGIdb database

The expression microarray datasets were retrieved and downloaded from the Gene Expression Omnibus (GEO) database ( https://www.ncbi.nlm.nih.gov/geo/ ) and “diabetic kidney disease” or “diabetic nephropathy” was used as the search term. The selection criteria included “ Homo sapiens ,” tissue samples, and transcriptomic data. The training datasets were GSE96804 and GSE30122 , whereas the validation dataset was GSE104948 (Table 1 ). The Nephroseq database ( http://v5.nephroseq.org/ ) integrates publicly available kidney gene expression profiles and clinical information. The group was selected as “diabetic nephropathy” and the tissue type was selected as “kidney.” Ultimately, the glomerulus dataset (Ju CKD Glom) and tubular dataset (Ju CKD TubInt) were chosen for validation (Table 1 ). The GeneCards database ( https://www.genecards.org/ ) was searched using “efferocytosis” as the keyword, resulting in 144 genes with a relevance score ≥ 0.2. In addition, 137 ERGs were extracted from previous reviews [ 17 – 19 ]. Ultimately, 256 ERGs were identified (Table S1 ). Table 1 Overview of the datasets Dataset Database Platform Type GSE96804 GEO GPL17586 41 DKD cases and 20 normal cases GSE30122 GEO GPL571 19 DKD cases and 50 normal cases GSE104948 GEO GPL22945 7 DKD cases and 18 normal cases Ju CKD Glom Nephroseq / 12 DKD cases and 21 normal cases Ju CKD TubInt Nephroseq / 17 DKD cases and 31 normal cases Overview of the datasets The Sangerbox online analysis platform ( http://www.sangerbox.com/ ) was utilized for this analysis [ 20 ]. Batch effects in the GSE96804 and GSE30122 datasets were removed using the combat algorithm, and gene expression matrices were obtained. DEG analysis was performed by the limma package [ 21 ]. The criterion for identifying DEGs was adj. P 1.5. An FC threshold > 1.5 was set for DEGs, which is consistent with the screening criteria used in recent similar studies [ 22 ]. This threshold effectively identifies biologically significant expression changes while maintaining statistical rigor. Additionally, an adj. P < 0.05 was used, with the false discovery rate (FDR) corrected P value calculated by the Benjamini‒Hochberg method. Volcano plots and heatmaps were generated to visualize DEG expression. The DEGs obtained were intersected with 256 ERGs to form DKD efferocytosis-related genes (DKD-ERGs). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for DKD-ERGs were carried out by the clusterProfiler package [ 23 ].The selected DKD-ERGs were uploaded to the STRING database ( https://string-db.org/ ) to construct a protein‒protein interaction (PPI) network. We subsequently used the cytoHubba plugin in Cytoscape 3.9.1 to identify the top 10 genes [ 24 ]. In the network, nodes represent genes, and the degree value of a node is determined by the number of edges connected to it. A higher degree value indicates greater importance of the node. The top 10 targets ranked by degree value were identified as key targets. Given the high redundancy and limited diagnostic specificity of gene sets obtained from traditional differential analysis, we further employed a machine learning framework for secondary screening to identify feature genes with both robustness and potential clinical translational value [ 25 , 26 ]. The DKD-ERGs were filtered using three machine learning methods: least absolute shrinkage and selection operator (LASSO) regression, random forest (RF), and the BORUTA algorithm. LASSO regression improves the interpretability of the model through regularization [ 27 ]. The LASSO algorithm was implemented using the glmnet package [ 28 ], with the response variable set as binomial and the alpha parameter set to 1. The model was validated using 10-fold cross-validation. RF enhances the predictive accuracy by constructing multiple decision trees. RF was implemented using the randomForest package [ 29 ], with 5-fold cross-validation used to assess model performance. The Boruta algorithm was implemented using the Boruta package [ 30 ], which selects significant features by evaluating the relative importance of features and random variables. The genes selected on the basis of degree value were intersected with those identified by machine learning methods to obtain hub genes. The expression and predictive value of the hub genes were observed in the training set. The predictive performance was assessed using the pROC package to calculate the area under the curve (AUC) [ 31 ]. The hub genes were subsequently evaluated on the basis of the GSE104948 dataset and the Nephroseq database, and the correlations between genes and clinical parameters (serum creatinine (Scr) and the glomerular filtration rate (GFR)) were explored. All animal experimental procedures were performed in accordance with ARRIVE guidelines. A total of 24 male SPF-grade C57BL/6J mice, with an average weight of 20 ± 2 g and aged 6–8 weeks, were purchased from Beijing Vital River Laboratory Animal Technology Co. Ltd. (SCXK (Beijing) 20210006). This study received approval from the Ethics Committee of Beijing University of Chinese Medicine (BUCM-2023120104-4282). Diabetes was induced using a 60% high-fat diet (HFD) followed by streptozotocin (STZ) administration. The mice in the control group ( n = 12) were fed a conventional diet, but those in the model group were fed a HFD ( n = 12). After 4 weeks, the mice in the model group were given intraperitoneal injections of STZ (50 mg/kg/day) for 5 consecutive days, whereas the control mice received an equivalent volume of citrate buffer via intraperitoneal injection. Mice with a fasting plasma glucose (FPG) concentration > 16.7 mmol/L in three consecutive tail-tip blood draws were considered to have successfully established a diabetic model. The mice were then fed a HFD for 20 weeks. An increase in urinary microalbumin levels is an important indicator of DKD [ 32 ], and morphological kidney lesions can be directly observed using HE, PAS and Masson staining (Fig. 2 ). Under 1% pentobarbital sodium anesthesia administered intraperitoneally (ip), blood was collected, and the serum was separated by centrifugation at 3000 rpm for 10 min and then stored at −80 °C. The kidneys were carefully excised and weighed, and the capsule was removed to isolate the cortex. One part of the tissue was preserved in 4% paraformaldehyde for pathological staining, while the other part was rapidly snap-frozen in liquid nitrogen and sent to Shanghai Majorbio Biopharm Technology Co., Ltd. for transcriptomic, proteomic, and metabolomic sequencing. The detailed sample preparation, instrument parameter settings, and data collection and analysis can be found in the Supplementary material . Fig. 2 Flowchart of the construction of the STZ combined with HFD-induced DKD mouse model Flowchart of the construction of the STZ combined with HFD-induced DKD mouse model Further analysis of the changes in ANXA1, CASP3, and C3 in the serum of DKD patients was conducted. The detection of serum ANXA1 and CASP3 is based on clinical Study I, whereas the detection of serum C3 is based on clinical Study II. Clinical Study I included patients diagnosed with type 2 diabetes mellitus (T2DM) and DKD who were admitted to Dongzhimen Hospital, Dongcheng Campus, Beijing University of Chinese Medicine, from May 2021 to March 2023. The study was approved by the Ethics Committee of Dongzhimen Hospital, Beijing University of Chinese Medicine (2022DZMEC-062-03), and all participants provided written informed consent. The inclusion criteria, exclusion criteria, and DKD staging criteria are presented in the Supplementary material . On the second day of hospitalization, fasting venous blood samples (2 mL) were collected from each participant in the morning. After standing at room temperature for 30 min, the samples were centrifuged at 3000 rpm for 10 min at 4 °C. The resulting serum was carefully separated, aliquoted into labeled EP tubes, and stored at −80 °C for future analysis. The serum ANXA1 levels were measured using enzyme-linked immunosorbent assay (ELISA) kits by the Annexin A1 ELISA Kit (ab222868), whereas the serum CASP3 levels were measured using the Human CASP3 (Caspase-3) QuickTest ELISA Kit (FineTest, QT-EH0546). In this study, we retrospectively collected data from 128 DKD patients who underwent serum C3 testing and were hospitalized in the Nephrology Department of Dongzhimen Hospital, Beijing University of Chinese Medicine, from January 2018 to December 2022. This study was approved by the Ethics Committee of Dongzhimen Hospital, Beijing University of Chinese Medicine (2023DZMEC-481-01), and adhered to the principles outlined in the Declaration of Helsinki. Owing to the retrospective nature of the study, the Ethics Committee of Dongzhimen Hospital, Beijing University of Chinese Medicine waived the need to obtain informed consent. The inclusion criteria, exclusion criteria, and DKD staging criteria are presented in the Supplementary material . Single-cell sequencing data from DKD patients and healthy controls, which are available in the Kidney Integrative Transcriptomics (K.I.T.) database ( http://humphreyslab.com/SingleCell/ ), were subsequently analyzed. Wilson et al. conducted unbiased snRNA-seq on cryopreserved human DKD tissues. On the basis of this database, our study extracted single-cell sequencing data for three hub genes from diabetic nephropathy lesion samples and visualized the results [ 33 ]. Next, we constructed a nomogram for the hub genes by the rms package [ 34 ]. The nomogram analysis predicts outcomes on the basis of key factors in the study and visually presents the complex relationships between multiple variables and the outcomes. After obtaining the risk score, we classified the patients into high-risk or low-risk groups on the basis of the median risk score. Differential expression analysis and enrichment analyses were performed on the basis of new groups. The criterion for identifying DEGs was adj. P 1.5. On the basis of the above model, validation was conducted on both the validation set and the Nephroseq database. Additionally, correlation analysis between the risk score and both the serum creatinine level and the eGFR was performed to further enhance the clinical applicability of the findings. Consensus clustering is a method for identifying molecular subtypes on the basis of approximate clustering numbers and is used to discover DKD subgroups related to the expression of apoptosis-associated regulatory factors through the k-means method. On the basis of the expression data of ANXA1, CASP3, IL33, and C3 from 60 DKD patients in the training set, an unsupervised clustering analysis of DKD samples was performed using the ConsensusClusterPlus package [ 35 ]. The similarity among the samples was determined via Euclidean distance, which was subsequently subjected to K-means clustering. Next, 500 iterations were conducted with a resampling proportion of 0.8. PCA was then employed to validate different efferocytosis molecular subtypes. To evaluate immune infiltration across different subtypes, we employed the CIBERSORT algorithm (with the “PERM” parameter set to 100 and a P value threshold of 0.05). The Drug–Gene Interaction Database (DGIdb) ( https://www.dgidb.org/ ), an online repository, integrates drug‒gene interaction data from multiple sources [ 36 ]. On the basis of previous studies and interaction scores, we selected the top 30 drugs as potential therapeutic agents for DKD [ 37 , 38 ]. Statistical analysis and plotting were performed using R Studio 4.4.3 and GraphPad Prism 9.5.1. For normally distributed data, the mean ± standard deviation was used, whereas nonnormally distributed data are presented as the median and interquartile range. For normally distributed data, independent sample t tests (for two groups) or one-way ANOVA (for ≥ 3 groups) were used. For nonparametric data, the Mann‒Whitney U test (for two groups) or the Kruskal‒Wallis test (for ≥ 3 groups) was used. For correlation analysis, Pearson’s correlation coefficient was used for normally distributed variables; for nonnormally distributed variables, Spearman’s rank correlation was used. To address the issue of multiple testing, we applied the Benjamini‒Hochberg method to adjust the two-sided p values and control the FDR. All the statistical tests were two-tailed, with a significance threshold set at P  < 0.05.

This study demonstrated the potential relevance of the efferocytosis process in DKD and identified 4 hub genes (ANXA1, CASP3, IL33, and C3), providing new insights into the early diagnosis and personalized treatment of DKD. These findings provide important theoretical foundations and practical guidance for further understanding the immunopathological mechanisms of DKD and exploring new therapeutic targets. Future research will further validate the clinical application of these genes and explore how regulating the efferocytosis process could improve the prognosis of DKD.

DKD exhibits significant clinical variability among individuals. Currently, there are insufficient sensitive and specific biomarkers for screening high-risk patient populations, highlighting the urgent need to explore more efficient diagnostic tools, predictive markers, and therapeutic strategies [ 43 ]. Cell death and the effective clearance of dying cells are key mechanisms for maintaining homeostasis in multicellular organisms [ 17 ]. Failed efferocytosis often leads to autoimmune or chronic inflammatory diseases [ 18 , 44 ]. Immune cell infiltration, which involves primarily macrophages (specialized phagocytes), can be observed in the renal tissue of DKD patients at various stages [ 45 ]. Failed efferocytosis may contribute to the progression of DKD [ 16 ] (Fig. 14 ). Fig. 14 Role of efferocytosis in DKD Role of efferocytosis in DKD This work methodically reveals the potential function of efferocytosis in DKD. Owing to the complexity of efferocytosis regulation in DKD, in the chronic inflammatory microenvironment, efferocytosis may not only be dependent on changes in gene expression levels but also be regulated by the tissue microenvironment [ 46 ]. Additionally, posttranslational modifications of receptor signaling pathways, such as the activation of TAM receptors and Rho GTPases, may also play a role in its regulation [ 17 , 18 , 47 ]. Ultimately, 15 genes were identified from the intersection of the DEGs and ERGs. Although the overlap is limited, this study successfully identified key targets, including ANXA1, CASP3, C3, and IL33, through multiomics integration. Multiomics sequencing was conducted on a DKD mouse model, validating findings at the transcriptomic and proteomic levels, whereas metabolomic analysis was employed to investigate the SPMs involved in the regulation of efferocytosis. Moreover, the hub genes also exhibited high diagnostic efficacy in external validation and were closely associated with kidney function markers. The application value of these hub genes in clinical diagnosis was further evaluated by measuring the circulating levels of ANXA1, CASP3, and C3 in DKD patients. Consensus clustering revealed that efferocytosis is closely related to immune status, suggesting that efferocytosis, as a process for clearing apoptotic cells, plays a critical role in immune homeostasis, tissue repair, and inflammation regulation. ANXA1 is one of the key SPMs involved in efferocytosis and plays a pivotal role in resolving inflammatory responses and restoring tissue integrity [ 18 , 48 ]. ANXA1 is upregulated in the renal tissues of patients with DKD, and high levels of renal ANXA1 are closely associated with elevated SCr levels and reduced eGFRs [ 32 , 49 ]. These findings suggest that ANXA1 may serve as a biomarker for predicting the severity of DKD [ 32 ]. Animal experiments further confirmed that AnxA1 gene knockout significantly exacerbates kidney injury. Conversely, the overexpression of AnxA1 in the kidney significantly alleviates renal damage. This phenomenon may be associated with the suppression of NF-κB p65 [ 32 ], the activation of the AMPK/PPARα/CPT1b signaling pathways [ 50 ] and Akt signaling [ 51 ]. Studies have indicated that serum ANXA1 levels are significantly lower in patients with DKD than in compared to individuals without kidney disease and healthy controls [ 49 , 52 ]. This aligns with our findings, where circulating ANXA1 levels gradually decrease as DKD progresses, and serum ANXA1 is associated with decreased eGFR and increased 24UTP. This discrepancy may be attributed to differences in sample selection, detection methods, and disease progression stages across studies. Additionally, we further analyzed the levels of classic SPMs in DKD kidney tissues using metabolomics. Our findings revealed the dysregulation of SPMs in DKD, highlighting the potential therapeutic role of SPMs in treating DKD. Furthermore, the ANXA1 mimetic, Ac2-26, is also considered a promising candidate drug for the treatment of DKD [ 53 ]. CASP3 primarily initiates programmed cell death by degrading key cellular substrates [ 54 ]. It is also involved in recruiting macrophages and regulating efferocytosis. When apoptotic cells are not efficiently engulfed by macrophages, they undergo secondary necrosis and release intracellular components that act as damage-associated molecular patterns, thereby triggering secondary inflammatory responses [ 55 ]. Studies have shown that CASP3 expression is upregulated in the renal tissues of DKD patients, potentially exacerbating renal injury by inducing apoptosis [ 56 ]. Drugs commonly used to treat DKD, such as saxagliptin [ 57 ], LCZ696 (valsartan/sacubitril) [ 58 ], and calcium dobesilate [ 59 ], can alleviate the pathological outcomes of DKD by inhibiting CASP3-mediated apoptosis. These findings suggest that therapeutic strategies targeting CASP3 hold promise for ameliorating renal injury in DKD. Clinical evidence has confirmed that blood CASP3 levels are significantly higher in patients with diabetes [ 60 , 61 ] and chronic kidney disease [ 62 ] than in healthy individuals. Therefore, we measured serum CASP3 levels in DKD patients. The results indicated that serum CASP3 levels were elevated in patients with DKD compared with those with T2DM and were significantly greater in patients with stage IV and stage V DKD than in those with stage III DKD. These findings suggest that serum CASP3 may have potential predictive value for DKD. IL-33 is a cytokine belonging to the IL-1 family [ 63 ]. When cells undergo damage or stress, IL-33 is swiftly secreted by the affected cells, activating both innate and adaptive immune mechanisms, which then drive proinflammatory reactions as well as tissue repair processes [ 64 ]. Compared with that in the kidneys of healthy organ donors, the level of IL-33 is notably increased in DKD patients [ 65 ]. The upregulation of IL-33 is negatively correlated with the eGFR. Blocking IL-33 signaling diminishes eosinophil infiltration in a DKD mouse model and alleviates glomerular injury and albuminuria, ultimately enhancing the repair of glomerular endothelial injury. A clinical study is also underway to assess the safety and efficacy of the tozorakimab (an IL-33 monoclonal antibody) in DKD patients [ 66 ]. However, there remains controversy regarding the circulating levels of IL-33 in DKD patients. Zhang et al. reported that the serum IL-33 concentration was greater in DKD patients than in healthy controls [ 67 ]. However, some small-sample studies suggest that serum IL-33 may not be a sensitive marker for early-stage DKD [ 68 ]. The expression of ST2L, the membrane-bound receptor for IL-33, is increased in the kidney tissues of individuals with diabetes [ 69 ]. Serum sST2 levels are elevated in CKD patients and are correlated with disease severity, but IL-33 levels are not significantly different [ 70 ]. New research provides an explanation for this contradiction, as direct detection of circulating IL-33 is not feasible [ 42 ]. Circulating IL-33 often binds to its receptor, soluble ST2 [ 71 ]. Further exploration of its mechanisms of action is necessary, along with improved experimental approaches to assess the feasibility of using circulating IL-33 in DKD evaluation.Further exploration of its mechanisms of action is necessary, along with improved experimental approaches to assess the feasibility of using circulating IL-33 in DKD evaluation. The complement system is a crucial part of the immune system and plays a pivotal role in enhancing phagocyte activity to eliminate microbes and damaged cells [ 72 ]. Excessive activation of the complement system at the local level contributes to the onset of DKD. C3 is a central component in complement activation, and transcriptomic analysis revealed that the relative expression of C3 in the glomerular tissue of DKD patients is six times greater than that in healthy individuals, closely correlating with the severity of glomerulosclerosis [ 73 ]. Deposition and activation of C3 in the glomeruli and renal tubulointerstitium accelerate the progression of DKD [ 74 – 76 ]. Renal complement C3 levels are negatively correlated with the GFR in DKD patients and may serve as an immune-related biomarker for DKD [ 77 ]. Data from DKD patients in this study revealed that as DKD progresses, circulating C3 levels gradually decrease. Consistent with our study, Zhang et al. performed a long-term follow-up study on 171 DKD patients diagnosed through renal biopsy and reported that those with reduced circulating C3 levels presented poorer renal function and more severe glomerular damage than patients with normal C3 levels did [ 78 ]. High C3 levels are an independent marker for nondiabetic renal disease in patients with T2DM [ 79 ]. The reduced circulating C3 levels in DKD patients may be a result of complement activation, leading to excessive consumption of C3 and its deposition in the kidneys or excretion in the urine [ 74 , 80 ], along with the massive production of complement cleavage products [ 81 ]. Given that efferocytosis may be a key regulatory factor in DKD, the nomogram risk score model shows good discriminative ability and predictive performance across different datasets. The risk score is significantly correlated with clinical indicators such as Scr and the eGFR, further indicating the potential application value of this scoring model in assessing DKD risk and disease severity. In a hyperglycemic environment, local inflammatory activation and immune infiltration in renal tissues, along with persistent, unresolved microinflammation, further exacerbate damage to the glomeruli and tubulointerstitium. This vicious cycle is a key factor driving the progression of DKD and related cardiovascular diseases [ 82 , 83 ]. Immune cells play an indispensable role in maintaining renal homeostasis [ 84 ]. Correlation analysis revealed significant associations between hub genes and immune cells, revealing a complex and tightly interconnected immune cell network in DKD pathology. These genes may participate in regulating the efferocytosis process, thereby influencing immune cell infiltration in the kidneys or periphery and leading to pathological changes in renal tissue. The efferocytosis process is strongly tied to immune infiltration in DKD, with immune infiltration playing a complex role in DKD pathogenesis. Many pathological mechanisms are involved in DKD development. However, over the past few decades, apart from new anti-diabetic drugs with significant cardiovascular and renoprotective effects, such as dapagliflozin, no new direct therapies for DKD have successfully been brought to market. A total of 102 potential drugs or compounds for DKD treatment were identified. The clinical translation of these drugs requires further investigation.

This study has several limitations. (1) This study is primarily based on bioinformatics analysis and animal model validation and lacks support from large-scale clinical trials. (2) This study utilized an STZ combined with HFD-induced DKD mouse model, whose disease progression, immune microenvironment, and metabolic characteristics differ from those of DKD. Transcriptomic analysis indicated that CASP3 was significantly upregulated in the kidney tissues of DKD patients, but no statistically significant difference was observed in the mouse model. This discrepancy may be due to differences in species-specific gene regulatory networks, disease duration, or methods of drug induction, and further investigation is warranted. (3) The overlap between DEGs and ERGs is relatively low, which may be attributed to the inclusion of posttranslational modification data, such as phosphoproteome data, thereby potentially overlooking crucial regulatory nodes. Future studies could integrate single-cell sequencing and spatial transcriptomics to better elucidate the functional heterogeneity of macrophage efferocytosis within specific renal microenvironments.

Diabetic kidney disease (DKD) is a critical microvascular complication associated with diabetes and one of the primary causes of end-stage renal disease (ESRD) [ 1 ]. Traditionally, hyperglycemia has long been considered the primary cause of DKD; however, recent clinical trial results suggest that hyperglycemia may not be the main pathogenic factor [ 2 , 3 ]. Strict glucose control may even increase mortality in DKD patients [ 4 ]. The development of DKD is driven by a variety of factors, such as inflammation, oxidative stress, autophagy, cell death, and the renin‒angiotensin‒aldosterone system [ 5 ]. Current treatments mainly focus on glucose control, renal protection, and blood pressure management. Nevertheless, the clinical prognosis of DKD remains unfavorable, primarily because of the absence of personalized and effective targeted therapies [ 6 ]. Early identification and effective risk assessment of DKD are crucial. Microalbuminuria is currently the most widely recognized marker for early screening of DKD, but it is crucial to recognize that renal damage can progress even in the absence of detectable microalbuminuria [ 7 ]. Consequently, there is an urgent demand for the development of innovative diagnostic methods and therapeutic targets to improve the clinical outcomes of patients with DKD. The interplay between immune cells and resident kidney cells is crucial in DKD [ 8 – 10 ]. Efferocytosis, the mechanism through which innate immune cells engulf apoptotic cells, is critical for maintaining immune balance, regulating chronic inflammation, and facilitating tissue repair [ 11 ]. In diabetes and its related complications (such as nonhealing wounds, osteoporosis, cardiomyopathy, periodontitis, and retinopathy), dysregulated efferocytosis is regarded as a central factor driving chronic inflammation and immune system dysfunction [ 12 – 15 ]. Whole-genome analysis in patients with DKD has highlighted the crucial role of impaired apoptotic cell clearance and the resulting inflammatory damage in the progression of renal disease [ 16 ]. How dysfunction of efferocytosis-related genes (ERGs) influences the occurrence and development of DKD remains an underresearched area. Therefore, systematically identifying genes linked to the efferocytosis process and uncovering their mechanisms in DKD will enhance our understanding of its molecular pathogenesis and offer new biomarkers for early diagnosis and personalized therapeutic strategies. In this study, we employed bioinformatics methods to compare the expression profiles of ERGs. Machine learning was subsequently applied to identify potential genes. External datasets and transcriptomic and proteomic data from a DKD mouse model were utilized as validation datasets. A risk model was constructed on the basis the identified hub genes, and two distinct subtypes of DKD were identified on the basis of ERGs. The immune infiltration levels of these subtypes were further analyzed. This research has potential for identifying novel targets for the diagnosis and therapeutic intervention of DKD (Fig.  1 ). Fig. 1 Schematic of the research workflow Schematic of the research workflow

Below is the link to the electronic supplementary material. Supplementary File 1 (DOCX 1.43 MB) Supplementary File 1 (DOCX 1.43 MB) Supplementary File 2 (DOCX 12.5 KB) Supplementary File 2 (DOCX 12.5 KB)