DN datasets were obtained from the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/), encompassing bulk-tissue messenger RNA (mRNA) sequence datasets (GSE30122, GSE30528, GSE30529, GSE96804, and GSE142153). Among these datasets, GSE30529 and GSE142153 served as independent validation sets. This comprehensive compilation included 92 DN samples juxtaposed with 93 NC samples. Furthermore, 2,794 human genes encoding proteins, with robust evidence supporting mitochondrial localization, and 2,001 genes associated with senescence were extracted utilizing GeneCards with a relevance score threshold of 1.0 (https://www.genecards.org/).

The merged datasets, namely GSE30122, GSE30528, and GSE96804, comprised the training set. The R package “sva” was used to remove batch effects from a training set. The “limma” R package (P < 0.05 and |log2 fold change (FC) | > 0.585) was utilized to identify DEGs in the training set.

WGCNA assumes paramount significance in bioinformatic analysis and is widely adopted for trait and gene association studies. To uncover important DNRGs, we first performed WGCNA based on the expression matrix of the bulk RNA-seq data. Genes located within the module exhibiting the strongest correlation with the DN phenotype were designated as DNRGs, reflecting their close association with DN based on co-expression patterns. These DNRGs were subsequently intersected with DEGs to identify key candidate genes for further analysis. The R package “WGCNA” was employed to construct a co-expression network.

Clustering analysis of the DN sample was conducted by “ConsensusClusterPlus” package, utilizing expression levels of DEGs. The optimal number of clusters was determined by evaluating the consensus matrix, consensus cumulative distribution function (CDF), and trace plots.

The MRGs were identified by intersecting genes encoding mitochondrial-related proteins, DNRGs, and DEGs. Similarly, SRGs were identified by intersecting genes encoding senescence related-proteins, DNRGs, and DEGs. To mitigate the high false discovery rates observed in previous studies, the identified genes were subjected to further filtering using the Wilcox test between DN and NC groups (P < 0.05). Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed.

MRGs and SRGs underwent further intersection and screening utilizing four machine learning algorithms to identify key genes. Support Vector Machine (SVM) was implemented using the “kernlab” R package. The random forest (RF) algorithm was implemented using the “randomForest” R package. Additionally, the XGBoost (XGB) algorithm, known for its efficient and scalable implementation of gradient boosting, was implemented using the “xgboost” package.

Mitochondrial dysfunction and senescence phenotype scores were calculated using the single-sample gene set enrichment analysis (ssGSEA) algorithm with the “GSVA” R package. Differences in phenotype scores between groups were assessed, and correlations among phenotype scores were examined. Subsequently, DN and control samples were segregated based on the median phenotype scores for subsequent GSEA involving all genes.

The performance of the diagnostic model for DN was assessed using the “pROC” package. Further analyses were conducted to identify the five most important features contributing to the model’s prediction accuracy. Moreover, the model was validated using an independent dataset, GSE142153, to evaluate its robustness and generalizability.

The scRNA-seq dataset GSE183276, comprising 9,235 cells from healthy kidneys and 27,929 cells from the DN cohort, was analyzed using the standard “Seurat” pipeline to assess single-cell characteristics. Cells were filtered based on stringent criteria: those with fewer than 400 genes, more than 5000 total genes, and over 30% mitochondrial genes were excluded. After quality control filtration, the R package “harmony” was employed to mitigate batch effects among samples. Subsequently, cell cluster annotation was performed based on previous research findings. Uniform Manifold Approximation and Projection (UMAP) was employed for visualization. To further clarify the intercellular communication patterns within the kidney microenvironment, we utilized the “CellChat” R package, which infers ligand–receptor interactions and predicts signaling networks across cell types.

The human renal tubular epithelial cell line (HK-2) (WHELAB C1116) was procured from Shanghai Whelab Bioscience Limited (Shanghai, China). HK-2 cells were cultured in Dulbecco’s Modified Eagle Medium/Nutrient Mixture F-12 (DMEM/F12) (Gibco, China) supplemented with 10% fetal bovine serum (Sigma-Aldrich, China) and 100 U/ml Penicillin-Streptomycin Liquid (Solarbio, Beijing, China) at 37 °C in a humidified incubator with 5% CO2. When the cell density reached 60-70%, glucose solution was added to the complete medium to achieve a final glucose concentration of 40 mmol/L. The cells were then incubated under the same conditions for 72 hours to establish a DN cell model (26). For the control group, an equal volume of phosphate-buffered saline solution was added to the complete medium, and the incubation conditions and duration were similar to those of the DN cell model group. After 72 hours, cells from the DN cell model and the control groups were collected for subsequent experiments.

siRNAs and the transfection reagent RNATransMate (E607402) were obtained from Sangon Biotech (Shanghai, China). HK-2 cells were seeded in six-well plates and cultured for 12 hours to reach a cell density of 60-70%. Next, appropriate amounts of siRNAs and RNATransMate were added to the six-well plates according to the manufacturer’s instructions to achieve a final siRNA concentration of 40 nM. After 4–6 hours of transfection, the medium was replaced with either a complete medium or a medium for constructing the DN cell model, and cells were cultured for an additional 72 hours. Subsequently, these cells were used for reverse transcription-quantitative polymerase chain reaction (RT-qPCR). The primer sequences used are provided in Supplementary Table 1 .

Total RNA was extracted from HK-2 cells using an RNA extraction kit (Magen, Shanghai, China) following the manufacturer’s protocol. The isolated RNA was reverse-transcribed into complementary DNA (cDNA) using the PrimeScript RT Master Mix (Takara, Japan). RT-qPCR was performed using the ABI Applied Biosystems PowerUp SYBR Green Master Mix (Thermo Fisher Scientific, USA) on a LightCycler 96 real-time PCR system (Roche Diagnostics Gmbh, Switzerland). β-actin served as the internal reference gene. Relative gene expression levels were calculated using the comparative cycle threshold (2^-ΔΔCt) method. The primer sequences used are provided in Supplementary Table 1 .

HK-2 cells were lysed in radioimmunoprecipitation assay buffer (Beyotime, Shanghai, China) supplemented with a protein phosphatase inhibitor. The lysates were centrifuged at 14,000 × g for 10 minutes to extract total protein. Proteins were separated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis and subsequently transferred onto polyvinylidene fluoride membranes. Membranes were then blocked with 5% Skim milk for the specified time. Subsequently, membranes were incubated at 4°C overnight with primary antibodies, followed by incubation with secondary antibodies at room temperature for 1 hour. The details of the antibodies used in this study are provided in Supplementary Table 2 . Protein bands were visualized using enhanced chemiluminescent reagents HY CEZMBIO West ECL-A and ECL-B (Hubei, China). The resulting bands were scanned, and their intensities were analyzed and quantified using ImageJ software and normalized to the corresponding internal controls.

After reaching a culture density of 70%, HK-2 cells were digested with trypsin and centrifuged at 3,000 × g for 5 minutes. Cells were resuspended in an electron microscopy fixation solution (Servicebio, Wuhan, China) and gently dispersed. Samples were fixed at room temperature for 30 minutes in the dark and then stored at 4 °C. Ultrathin sections were prepared and then stained with 2% uranyl acetate saturated in ethanol in the dark. Finally, the ultrathin sections were visualized under a transmission electron microscope (HITACHI HT7800/HT7700).

Differences in continuous variables were assessed using the Wilcoxon test. Pearson correlation analysis is used to assess the correlation between two continuous variables. Data were visualized using the R package “ggplot2”. All statistical tests were two-sided, and P < 0.05 was considered statistically significant. Significance levels were denoted as follows: ‘ns’ for not significant, *P < 0.05, **P < 0.01, and ***P < 0.001.

The datasets information and workflow of the presented study are shown in Figure 1 . In this study, we integrated three bulk RNA-seq datasets—GSE30122, GSE30528, and GSE96804—as a combined training cohort for differential gene expression analysis. Principal component analysis (PCA) revealed a clear distinction between DN and NC samples based on the expression profile of the training set ( Figure 2A ). Differential expression analysis revealed 2,176 DEGs between the DN and NC samples ( Figure 2B ). Subsequently, WGCNA analysis was conducted, wherein the sample clustering diagram displayed discernible patterns ( Figure 2C ), and a soft-thresholding power (β) of 14 yielded a scale-free topology with a correlation coefficient exceeding 0.85, supporting robust network construction ( Figure 2D ). Employing the dynamic tree cut algorithm, multiple gene modules were initially identified ( Figure 2E ), subsequently, we determined the correlations between each gene significance (GS) and module membership (MM), and retained four modules with relatively strong associations (GS > 0.2 and MM > 0.6) for further analysis. Among them, the brown module showed a significantly high correlation (R = 0.49) and significant P value (P = 3 × 10⁻⁹) with DN ( Figure 2F ). Consequently, genes within the brown module were designated as DNRGs. Notably, GO ( Figure 2G ) and KEGG ( Figure 2H ) functional enrichment analyses revealed that the differentially expressed DNRGs were primarily associated with pathways related to mitochondrial dysfunction, cellular aging, insulin resistance, and inflammation. Collectively, these results demonstrate that genes clustered within the brown module were strongly associated with DN and were predominantly involved in mitochondrial dysfunction, senescence, insulin resistance, and inflammatory processes, providing a robust gene set for subsequent biomarker selection and mechanistic investigations.

To investigate the diverse expression patterns of DEGs associated with mitochondrial dysfunction in DN, consensus clustering was applied. This analysis identified two distinct clusters, as shown in the consensus matrix (k = 2, Figure 3A ). Moreover, each cluster showed a score exceeding 0.8 at k = 2 ( Figure 3B ). The robustness of clustering was further confirmed based on minimal fluctuations of the consensus CDF curves across various consensus indexes ( Figure 3C ). Consequently, 60 DN samples were classified into two clusters: Cluster 1 (C1; 29 samples) and Cluster 2 (C2; 31 samples). PCA effectively separated the clusters ( Figure 3D ). WGCNA was performed to identify gene co-expression modules related to DN subtypes. Outlier samples identified through hierarchical clustering were excluded ( Figure 3E ). Following the selection of genes with the top 25% variance and exclusion of abnormal samples, a scale-free network was constructed with a soft threshold of 16, yielding a scale-free R²-value of 0.975 ( Figure 3F ). The WGCNA algorithm was employed to explore key genes associated with DN. The sample clustering diagram revealed four distinct co-expression modules ( Figure 3G ).

Among the identified gene modules, the brown module exhibited a significant correlation with mitochondrial dysfunction, with a coefficient of 0.81 and a P value of 3 × 10⁻¹⁵ ( Figure 3H ). In addition, the GS and MM were highly correlated in the brown module ( Figure 3I ). From this module, 461 genes were identified and selected for further analysis. Due to their significant positive association with mitochondrial dysfunction, these genes were designated as MRGs ( Supplementary Table 3 ).

To explore the diverse expression patterns of the DEGs associated with senescence in DN, consensus clustering was performed. This analysis identified two optimal clusters (k = 2), as shown in the consensus matrix ( Figure 4A ), with consensus scores exceeding 0.8 ( Figure 4B ) and stable consensus CDF curves ( Figure 4C ). The 60 DN samples were divided into C1 (33 samples) and C2 (27 samples), which were separated by PCA ( Figure 4D ). Outliers were excluded using a sample clustering tree ( Figure 4E ). Genes with the top 25% variance were selected and used to construct a scale-free network through the soft threshold of 16 ( Figure 4F ). The WGCNA test identified four co-expression modules ( Figure 4G ), with the brown module showing the strongest correlation with senescence (coefficient: 0.87, P = 3 × 10⁻¹⁹) ( Figure 4H ). A strong correlation was found between GS and MM in the brown module ( Figure 4I ), which identified 726 SRGs ( Supplementary Table 4 ).

To explore the potential interplay between mitochondrial dysfunction and senescence in DN, we analyzed the overlap between their respective gene modules. After intersecting the above MRGs, SRGs, and DEGs, 259 hub MRGs ( Figure 5A ) and 273 hub SRGs ( Figure 5B ) were identified. Subsequently, these genes were screened between DN and NC groups by the Wilcoxon test ( Figures 5C, D ). Further exploration of gene functions was conducted through GO and KEGG pathway enrichment analyses. The GO terms enriched among the hub MRGs included T cell activation, regulation of mitochondrial organization, and leukocyte homeostasis ( Figure 5E ). Notably, KEGG pathways encompassed the AGE-RAGE signaling pathway, phagosome, and leukocyte transendothelial migration, indicating a potential association between mitochondrial dysfunction and inflammaging in DN ( Figure 5F ). Similarly, GO terms enriched for hub SRGs involved regulation of T cell activation, leukocyte chemotaxis, aging, and mitochondrial organization ( Figure 5G ). KEGG analysis further highlighted enrichment in pathways such as the AGE-RAGE signaling pathway, chemokine signaling pathway, nuclear factor-kappa B signaling pathway, and leukocyte transendothelial migration, suggesting a potential link between senescence, mitochondrial dysfunction, and inflammation in DN ( Figure 5H ). Collectively, these enrichment results demonstrated significant functional overlap between hub MRGs and hub SRGs.

Furthermore, the ssGSEA analysis was performed to generate phenotype scores for mitochondrial dysfunction and senescence. As illustrated in Figure 6A , DN samples exhibited significantly higher phenotype scores for mitochondrial dysfunction and senescence (P < 0.001). The correlation between these two phenotype scores was remarkably strong (R > 0.9995) ( Figure 6B ). Based on these scores, hub MRGs and SRGs were divided into high- and low-expression groups. Differential expression and GSEA analyses were performed for both gene sets. Notably, mitochondrial dysfunction and increased senescence groups were significantly enriched in inflammatory pathways, such as macrophage activation, cytokine signaling, T cell activation, and chemokine responses ( Figures 6C, D ). These convergent findings further verify the strong relationship between mitochondrial dysfunction, senescence, and inflammation in the pathogenesis of DN.

To further investigate the potential mechanisms underlying mitochondrial dysfunction and senescence in DN, we intersected hub MRGs and hub SRGs identified previously. A total of 259 shared genes were identified from this intersection ( Figure 7A ), from which key genes with potential diagnostic value for DN were screened. Of note, 34 genes were selected for Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis based on the minimum λ value ( Figure 7B ; Supplementary Table 5 ). Subsequently, four machine-learning approaches RF, SVM, XGB, and generalized linear model (GLM) were employed to construct diagnostic models. Cumulative residual distribution plots ( Supplementary Figure 1A ) and residual boxplots ( Figure 7C ) revealed that SVM had smaller residual values, indicating a high model accuracy. Moreover, results shown in the Supplementary Figure 1B illustrate the top ten variables ranked by root mean square error (RMSE) for each model. Next, the receiver operating characteristic (ROC) curves were utilized to explore the diagnostic performance of the models, with three models demonstrating excellent discriminative power, with an area under the curve (AUC) exceeding 0.95 ( Figure 7D ). Among them, SVM revealed the best predictive accuracy and reliability, making it the most suitable diagnostic model for DN. The five most significant variables in the SVM model- CLDN1, TYROBP, HDAC9, CASP3 and RCN1, were selected as the hub genes for the diagnosis of DN. The expression of the identified diagnostic genes was further verified in both GSE142153 ( Figures 7E, G ) and GSE30529 ( Figures 7F, H ) datasets, consistent with the results from the training cohort. Notably, the diagnostic model exhibited robust predictive performance in both validation sets, with AUC values exceeding 0.8. Correlation analysis revealed potential relationship between the expression of HDAC9 and SASP markers ( Supplementary Figure 2 ). Furthermore, RT-qPCR was also conducted to validate the five hub genes in high glucose-stimulated HK-2 cell and NC samples ( Supplementary Figure 3 ), consistent with our previous bioinformatics findings.

Next, we characterized the single-cell landscape of DN. After rigorous quality control and batch effect correction, a total of 37,164 cells were yielded ( Supplementary Figure 4 ). Unsupervised clustering revealed 13 distinct cell types ( Figure 8A ), including B cells, distal convoluted tubule cells (DCT), endothelial cells (Endo), fibroblasts (Fib), intercalated cells (IC), loop of Henle cells (LOH), macrophages (Mac), monocytes (Mono), natural killer T cells (NKT), principal cells (PC), podocytes (Podo), proximal tubule (PT) cells, and T cells. The cell types were annotated based on their canonical marker genes, and the relative abundance of each cell type in shown in Figure 8B . Notably, PT cells were the most abundant in both DN and control tissues, with a predominance of PT cells originating from DN samples. The distribution of these cell types, along with their respective tissue origins, was visualized in the UMAP plot ( Figure 8C ).

Next, we analyzed DEGs between DN and NC samples across the 13 cell types based on their transcriptional profiles. Most DEGs were linked to “human diseases” pathways. Notably, the PT cells exhibited the highest number of downregulated DEGs compared to all other cell types ( Figure 8D , lower panel). Moreover, we investigated the activity of mitochondrial dysfunction- and senescence-related gene sets using AUCell scoring ( Figures 8E, F ). Both AUC scores were significantly elevated in PT cells, indicating that these pathways were mainly involved in DN pathogenesis, particularly within PT cells. In addition, HDAC9 expression was significantly upregulated in PT cells from DN samples ( Supplementary Figure 5 ). To investigate cell-to-cell communication patterns, we applied the “CellChat” R package and identified that PT cells, LOH, DCT, PC, and B cells shared common intercellular signaling patterns in DN. Specifically, macrophage migration inhibitory factor (MIF) and growth differentiation factor (GDF) signaling pathways were enriched in pattern 2 ( Figure 8G ). Previous studies have implicated both pathways in regulating mitochondrial dysfunction and senescence (27, 28). Collectively, these findings suggest that HDAC9 may play a role in the progression of DN by regulating mitochondrial dysfunction, enhancing inflammatory responses, and promoting cellular senescence, especially in PT cells.

To identify the role of HDAC9 and PT cells in DN, HDAC9 inhibitor and high glucose-stimulated HK-2 cell models were used. As expected, the expression of HDAC9 was significantly upregulated in HK-2 cells following exposure to a high glucose environment compared to untreated HK-2 cells, which is in line with findings from our previous study ( Figure 9A ). As shown in Figure 9B , HK-2 cells exposed to high-glucose conditions exhibit significant mitochondrial structural abnormalities. These include a reduction in mitochondrial matrix density, decreased electron density, and the appearance of vacuole-like structures. Additionally, the intermembrane spaces are enlarged, and the cristae become stubby or swollen, leading to smaller, rounded mitochondria with disorganized or even absent cristae. These morphological changes reflect mitochondrial dysfunction. Senescent cells secrete various profibrotic and proinflammatory factors collectively referred to as the senescence-associated secretory phenotype (SASP). Cellular senescence plays a crucial role in the pathogenesis of DN. Established senescence markers, including p53, p21, and p16, are commonly employed to identify and characterize senescent cells (29, 30). Here, we observed a significant increase in the expression of p53, p21, and p16 ( Figure 9C ), as well as significantly upregulated mRNA levels of SASP components (IL-6, IL-1β, IL-8, and TNFα) in high glucose-stimulated HK-2 cells ( Figure 9D ). The increase in these SASP factors highlights the potential involvement of cellular senescence in DN progression.

To further verify the role of HDAC9 in DN, siRNA-mediated knockdown of HDAC9 was performed. Data shown in Figure 9E confirm the successful knockdown of HDAC9 expression in HK-2 cells by with a knockdown efficiency of approximately 80%. Following HDAC9 silencing, the expression levels of SASP components (IL-6, IL-1β, IL-8, and TNFα) ( Figures 9F–I ), as well as senescence markers p53, p21, and p16 ( Figures 9J–L ), were significantly downregulated in high glucose-treated HK-2 cells compared to the corresponding controls.