Multi-Omics and Experimental Validation Reveal the Targets and Mechanisms of Urolithin A in diabetic kidney disease

1. Introduction

Diabetic kidney disease (DKD) is characterized by both its high prevalence (30%-40% of people with diabetes) and its role as the leading cause of end-stage renal disease (ESRD). It substantially diminishes patients’ quality of life and places a heavy healthcare burden on medical systems [1,2]. Although current management strategies include strict glycemic and blood pressure control along with novel agents such as SGLT2 inhibitors and RAAS blockers, many patients still develop persistent proteinuria and progressive renal function decline, often advancing to ESRD [3]. DKD pathogenesis is multifactorial in nature, driven by metabolic disturbances, hemodynamic abnormalities, oxidative stress, inflammatory responses, and genetic factors [4]. The lack of highly sensitive and specific biomarkers often leads to delayed diagnosis until irreversible stages. Therefore, elucidating the molecular mechanisms of DKD, identifying new therapeutic targets, and developing targeted interventions are crucial for early precision-based prevention and delaying the progression of the disease.

Urolithins are bioactive metabolites derived from the gut microbiota-mediated metabolism of dietary polyphenols such as ellagitannins and ellagic acid [5]. Since their discovery in 1980, the physiological functions of urolithins have garnered increasing interest [6]. Evidence indicates that urolithins, including Urolithin A (UA) and Urolithin B, confer benefits in regulating oxidative stress, suppressing inflammation, enhancing mitophagy, and modulating energy metabolism [7,8]. In diabetic complications, UA has attracted attention due to its potent antioxidant and anti-inflammatory properties, suggesting a protective role in DKD. Jiang et al. demonstrated that ellagic acid and UA significantly reduced blood glucose, serum creatinine, and TNF-α levels in diabetic mice, inhibited the renal HMGB1/TLR4/NF-κB pathway, attenuated oxidative damage, and improved renal histopathology, indicating a potential mechanism for delaying DKD progression [9]. In podocytes exposed to high glucose, UA played a protective role by lowering reactive oxygen species (ROS), boosting nephrin expression and endocytosis, and counteracting foot process effacement, as reported by Kotewicz et al. [10]. Furthermore, Kalinowski et al. showed that UA modulates the TGF-β pathway, ameliorates podocyte epithelial-mesenchymal transition and fibrosis, and retards DKD progression [11]. Although these studies suggest that UA exerts renal protection through anti-inflammatory, antioxidant, and podocyte-protective mechanisms, its specific molecular targets and detailed mode of action in DKD remain unclear. A deeper understanding of UA’s targets and signaling networks is a prerequisite for developing novel DKD therapies.

Notably, accumulating single-cell transcriptomic evidence has revealed that proximal convoluted tubule (PCT) cells, which constitute a major component of the renal parenchyma, undergo profound subtype-specific transcriptional remodeling under diabetic conditions [12]. Proximal tubular epithelial cell dysfunction is increasingly implicated in the pathogenesis and progression of DKD, and some studies suggest that proximal tubular injury may even precede microalbuminuria and early glomerulopathy [13]. In disease states, the distinction between anatomic and physiological PCT subtypes is blurred or erased, making it critical to characterize PCT state identities at single-cell resolution [14]. However, despite growing recognition of PCT subtype heterogeneity in DKD, the cell-type-specific expression of UA’s molecular targets within distinct PCT subpopulations and their roles across PCT differentiation stages have not been characterized, leaving the cellular basis of UA’s renoprotective actions poorly understood.

New avenues for targeted therapy and diagnosis are emerging, enabled by high-throughput sequencing and bioinformatics, which have deepened our understanding of complex disease mechanisms and identified numerous risk genes [15]. Network pharmacology systematically elucidates interactions among multi-component drugs, targets, and pathways, revealing potential mechanisms of action [16]. As a key computational drug design tool, molecular docking predicts ligand–target binding modes and affinities, providing structural support for target validation [17]. Recent integration of single-cell transcriptomics, high-throughput omics, and multimodal bioinformatics has identified key differentially expressed genes in DKD, some of which are proposed as biomarkers for early diagnosis or prognosis [18]. However, due to the complexity and cellular heterogeneity of DKD, targeted mechanistic studies investigating natural compounds such as UA remain scarce.

Based on this background, this study combines transcriptomic and single-cell data with network pharmacology to identify UA targets in DKD and employs molecular docking to assess binding affinities. Furthermore, we examined the cell-type-specific expression and intercellular communication of these targets at single-cell resolution, providing a theoretical basis for the renal protective mechanisms of UA. The overall workflow is summarized in Figure 1.

2. Results

2.1. 18 Candidate genes affecting DKD development

There were 4,256 DEGs between the DKD and control samples in the training set, including 2,860 up-regulated genes and 1,396 down-regulated genes (Figure 2a). A total of 57 target genes for UA were obtained (Figure 2b). To further identify DEGs associated with UA, 4,256 DEGs and 57 target genes were intersected, and 18 intersection genes were obtained as candidate genes (Figure 2c). Moreover, GO enrichment analysis was performed on candidate genes. Among the 15 significantly enriched results (P < 0.05) were the one-carbon metabolic process, protein autophosphorylation, protein tyrosine kinase activity, and membrane raft (Figure 2d, Table S1). Furthermore, the pathways for KEGG enrichment were nitrogen metabolism, cell cycle, MAPK signaling pathway, FoxO signaling pathway, and so on (P < 0.05) (Figure 2e, Table S2). In the PPI network, there were 13 interacting proteins. Among them, KDR interacts with other proteins, such as CCNA2, CA9, and ESR2 (Figure 2f). Overall, the comprehensive analysis has offered a wealth of information that could guide further research efforts aimed at deciphering the underlying mechanisms of DKD and devising innovative therapeutic approaches.

2.2. A total of 4 biomarkers that treat DKD were discovered

First, based on the candidate genes, 5 candidate feature genes were acquired by LASSO regression analysis with the log (lambda) value being -2.784865 (Figure 3a). In addition, through the SVM-RFE algorithm, based on the importance ranking, a total of 8 characteristic genes were obtained (Figure 3b). Finally, 5 intersection genes of the 2 algorithms (HSPA1A, ESR2, NUAK1, CA9, and PLK1) were obtained (Figure 3c). Although all five genes were significantly dysregulated in the training set (P < 0.001; Figure 3d), validation set analysis confirmed significant expression differences for HSPA1A, ESR2, NUAK1, and PLK1 only (P < 0.01; Figure 3e), all of which displayed consistent expression trends. Among them, ESR2, NUAK1, and PLK1 were upregulated in DKD, while HSPA1A was downregulated in DKD compared with control samples. Therefore, these 4 genes were utilized as candidate biomarkers.

The four genes HSPA1A, ESR2, NUAK1, and PLK1 all demonstrated significant diagnostic value for DKD, with AUC values greater than 0.7 following ROC curve analysis (Figure 3f-g). These outcomes showed that the diagnostic efficacy of these 4 genes was satisfactory, and thus they could serve as biomarkers for subsequent analysis.

2.3. Functional enrichment and immunoregulation of biomarkers in DKD

Among the four biomarkers, ESR2 was associated with 100 enriched pathways based on GSEA. A visualization of the top 10 pathways featured ascorbate and aldarate metabolism, butanoate metabolism, and the TCA cycle (Figure S1a, Table S3). HSPA1A was enriched in 85 pathways, including ascorbate and aldarate metabolism, citrate cycle (TCA cycle), and propanoate metabolism (Figure S1b, Table S4). NUAK1 was associated with 101 enriched pathways, including ascorbate and aldarate metabolism, butanoate metabolism, and the citrate cycle (TCA cycle) (Figure S1c, Table S5). Similarly, PLK1 enrichment was observed in 99 pathways, which also featured the same key metabolic processes (Figure S1d, Table S6). The above results indicated that there was significant overlap in the enriched pathways among the biomarkers, all of which were enriched in pathways such as ascorbate and aldarate metabolism, butanoate metabolism, citrate cycle (TCA cycle), and degradation of valine, leucine, and isoleucine, and fatty acid metabolism. This suggested that the drug targets had certain similarities in terms of function. Collectively, these findings laid a solid foundation for elucidating the molecular mechanisms of DKD.

The composition of 22 immune cell types in DKD and control groups is presented in Figure S1e. After excluding cells with a proportion of 0 in all samples, 18 immune cells were finally analyzed for differences. There were 8 immune cells with marked differences between DKD and control. The outcomes revealed that compared with the control, the DKD samples had a higher proportion of CD8 T cells and a lower proportion of naive B cells (P < 0.01) (Figure S1f). Among them, naive B cells and T follicular helper cells showed the highest positive correlation (cor = 0.57, P < 0.001); naive B cells and resting mast cells showed the most pronounced negative correlation (cor = -0.60, P < 0.001) (Figure S1g). ESR2, NUAK1, and PLK1 showed strong positive correlations with CD8 T cells (cor > 0.7, P < 0.001) but strong negative correlations with T follicular helper cells (cor < -0.5, P < 0.001). In contrast, HSPA1A exhibited an inverse relationship, correlating negatively with CD8 T cells (cor = -0.65, P < 0.001) and positively with T follicular helper cells (cor = 0.52, P < 0.01) (Figure S1h, Table S7). The above results suggest that biomarkers might influence the immune cell infiltration of DKD, and this might provide a reference for the clinical treatment of DKD.

2.4. Biomarker-UA binding affinity and clinical renal correlation

Molecular docking revealed a strong binding affinity between UA and the biomarkers HSPA1A, ESR2, NUAK1, and PLK1, as evidenced by binding energies of -8.8, -8.8, -7.6, and -7.0 kcal/mol, respectively (Table S8). All values were below the -5 kcal/mol threshold, confirming stable interactions. The binding structure analysis revealed that UA infiltrated into the binding site of each biomarker, and there were abundant hydrogen bond donors and acceptors around the vicinity of the binding sites (Figure S2a). A significant positive correlation was observed between PLK1 expression and serum creatinine levels in patients (cor = 0.7005, P = 0.01117). No other biomarkers, however, showed significant correlations with glomerular filtration rate (GFR) (Figure S2b-c). This could potentially offer a valuable reference for the clinical management of DKD.

2.5. Identification of 11 renal cell types from scRNA-seq data

Quality control was performed in the single-cell dataset GSE131882, and after the removal of some cells, the nFeature_RNA and nCount_RNA distribution maps after quality control were displayed (Figure S3a). Thereafter, 2,000 HVGs were identified and the 10 most important genes were labeled, such as PLA2R1, SLC26A4, and PTPRQ (Figure S3b). The top 10 principal components with statistical significance were determined according to the inflection point of the Elbow plot for subsequent single-cell analyses (P < 0.05) (Figure S3c-d). As shown, cells were divided into 17 clusters (Figure S3e); all clusters were annotated as 11 cell types, such as proximal convoluted tubules (PCT), endothelial cells (ENDO), podocytes (PODO), and mesangial cells (MES) (Figure S3f). A dot plot depicting marker gene expression across the various cell types is shown in Figure S3g.

2.6. Functional enrichment and cell communication in DKD

Functional enrichment analysis of all cells was conducted to identify associated biological pathways. For example, COX reactions were upregulated in collecting duct-intercalated A cells (CD-ICA), but downregulated in collecting duct-intercalated B cells (CD-ICB), PCT, and ENDO. Potassium transport channels were upregulated in loops of Henle (LOH), collecting duct-principal cells (CDPC), and distal convoluted tubules (DCT) (Figure S4a).

The cell communication network was shown in the figure (Figure S4b-c). The left side showed the number of cell communications. In the control samples, ENDO and other cell types had a large number of interactions, indicating their important role in cell-cell crosstalk. In the DKD samples, PODO had numerous interactions with other cell types. On the right side was the cell communication probability; in the control samples, only ENDO and DCT had a higher communication probability with the rest of the cells, demonstrating their higher active degree in cell communication. In DKD samples, there was a high probability of communication between DCT and CD-ICA. Receptor-ligand interaction plots showed that in control samples, PODO-ENDO had the highest probability of receptor-ligand pairing with VEGFA-VEGFR1, VEGFA-VEGFR1R2, and VEGFA-VEGFR2 (Figure S4d). In DKD samples, PODO-ENDO still showed an extremely high interaction probability with VEGFA-VEGFR1, suggesting that this molecular pair might play an important role in cell-cell interactions (Figure S4e). By understanding the mechanism of signal communication between cells, it was hoped to open up a new direction for the treatment of DKD.

2.7. Biomarkers associated with proximal tubule differentiation in DKD

Comparative analysis of cell abundance between the DKD and control groups identified a total of 9 cell types with significant abundance differences. Among them, the proportions of CDPC, LOH, etc., in the disease group were significantly higher than those in the control group, while the proportions of DCT and PCT in the DKD group were significantly lower than those in the control group (P < 0.05) (Figure 4a). In addition, based on the expression of biomarkers in differential cells, after removing cells with no biomarker expression, 5 differential cells remained, including LOH, CDPC, DCT, PCT, and complement factor H (CFH). Among them, both biomarkers NUAK1 and PLK1 were differentially expressed in PCT (P < 0.05), and literature reports indicated that PCT played an important role in DKD (Figure 4b). Analysis of the GSE209781 dataset confirmed differential expression of HSPA1A and NUAK1 in PCT cells (P < 0.01), with expression trends for all biomarkers being consistent with the GSE131882 dataset (Figure 4c). So the PCT was used as the key cell for the subsequent analysis.

Re-annotation resulted in the subdivision of PCT cells into three subtypes (PTH1R PCT, HAO2 PCT, CLCNKA PCT). A dot plot was generated to visualize the marker gene expression across these different cell subtypes (Figure 4d-e). Key cell pseudotime analysis showed that the cell developmental trajectory included 9 states: state 1 was in the early developmental stage, and parts of state 9 were in the late developmental stage (Figure 4f). The cell pseudotime density map showed that HAO2 PCT and PTH1R PCT had higher density at the cell developmental starting point during a developmental cycle, while CLCNKA PCT were more at the developmental starting point and then in large numbers at the developmental endpoint (Figure 4g). The expression of biomarkers in different developmental states of key cells showed that HSPA1A, ESR2, and NUAK1 were mainly highly expressed in the late stage of cell development, while PLK1 was highly expressed in the early stage of cell development (Figure 4h). The above results suggested that the biomarkers might affect DKD by interacting with key cells, which provided a reference for the treatment of DKD.

2.8. General condition and biochemical parameters in diabetic mice

At 12 and 16 weeks, db/db mice exhibited significantly increased body weight and kidney-to-body weight ratio compared to db/m controls (Figure 5a-b). UA intervention reduced body weight and attenuated renal hypertrophy in db/db mice. Biochemical analysis showed markedly elevated serum creatinine, blood urea nitrogen, and total cholesterol levels in db/db mice at 12 weeks (P < 0.05). These metabolic abnormalities were further exacerbated by week 16, with additional significant increases in 24-hour urinary protein, serum albumin, and triglycerides. UA treatment effectively ameliorated all these diabetic metabolic alterations, significantly reducing each parameter measured (Figure 5c-h).

2.9. UA treatment ameliorates renal morphological alterations in diabetic mice

Examination of renal tissues from the db/m group by morphological staining showed no obvious abnormalities. In contrast, 12 and 16-week db/db mice, H&E staining showed glomerular hypertrophy with compensatory tubular dilation and interstitial inflammation; Masson’s staining demonstrated diffuse-to-focal glomerular basement membrane thickening and increased extracellular matrix; PAS staining confirmed glomerular expansion, tubular degeneration/atrophy, and mesangial matrix accumulation. Transmission electron microscopy revealed ultrastructural alterations in diabetic mice, including extensive foot process effacement, reduced podocyte number, and significant glomerular basement membrane thickening. UA intervention markedly ameliorated these histopathological and ultrastructural impairments at both time points versus age-matched db/db controls (Figure 6).

2.10. Expression level of hub genes in diabetic mice renal tissue

RT-qPCR analysis of hub gene expression in renal cortical tissue revealed that, compared to the db/m group, PLK1 and NUAK1 were significantly up-regulated, while HSPA1A was significantly down-regulated in the db/db group (all P < 0.05). These findings are consistent with the transcriptomic data. After UA intervention, expression of PLK1 and NUAK1 in the db/db+UA group was significantly reduced, whereas HSPA1A expression markedly recovered, compared to the db/db group. No significant difference in ESR2 mRNA expression was observed among groups (P > 0.05) (Figure 7a-d). The expression level of PLK1 was significantly correlated with serum creatinine, a key biomarker of renal impairment. NUAK1 showed a consistent expression pattern with PLK1 and exhibited binding affinity with urolithin A in molecular docking assays. Based on emerging evidence supporting roles for both genes in diabetic kidney disease [33,34], we validated their protein expression in kidneys from db/db mice.Western blot analysis revealed significantly elevated protein levels of NUAK1 and PLK1 in db/db mice compared with db/m controls, which were attenuated by UA treatment (Figure 7e-h). A similar trend was observed by Immunohistochemical analysis (Figure 7i-j).

2.11. UA restored the altered expression of NUAK1 and PLK1 in HK-2 cells exposed to HG

Using CCK-8 assays, we determined that 10 μM UA was the optimal concentration for subsequent experiments (Figure 8a). Similar effects were observed in HK-2 cells. Immunofluorescence staining demonstrated increased fluorescence intensity of NUAKk1 and PLK1 under HG conditions, which was reversed by UA treatment (Figure 8b-e). These results suggest that UA confers protective effects in HK-2 cells and db/db mice, through downregulating NUAK1 and PLK1.

3. Discussion

DKD is a common and severe microvascular complication of diabetes with a complex pathogenesis. Currently, specific diagnostic methods and effective treatments remain limited. Therefore, exploring novel molecular mechanisms, identifying key genes and biomarkers, and developing small molecule inhibitors are crucial for early DKD intervention [19,20]. Advances in high-throughput sequencing and artificial intelligence have enabled biomarker and therapeutic target discovery using multi-omics data. Previous RNA sequencing studies in DKD relied mainly on bulk tissue transcriptomes, lacking resolution for renal cellular heterogeneity. This study integrated transcriptomic and single-cell sequencing with network pharmacology, identifying 18 differentially expressed genes linked to DKD and UA. KEGG enrichment indicated these genes were involved in AMPK and FoxO signaling pathways. Further machine learning and protein-protein interaction analyses identified four hub genes: PLK1, NUAK1, HSPA1A, and ESR2. ROC curve analysis and correlation with clinical markers (e.g., serum creatinine, eGFR) supported their validity as biomarkers for DKD diagnosis and prognosis. GSEA revealed these genes collectively regulate metabolic reprogramming and redox homeostasis, suggesting UA ameliorates DKD through modulation of energy metabolism and antioxidant responses via these targets. Experimental validation confirmed upregulation of PLK1 and NUAK1 and downregulation of HSPA1A in a diabetic model, which were reversed by UA treatment.

NUAK1 encodes a serine/threonine kinase within the AMP-activated protein kinase (AMPK) family and regulates key processes including cell survival, metabolism, and tumorigenesis [21,22,23]. It acts as a TGF-β-induced pro-fibrotic kinase, upregulated in multiple organ fibrosis models, where it promotes fibrotic scarring through YAP and TGF-β/SMAD pathway activation. Activated YAP/TAZ promotes NUAK1 transcription, creating a pathological positive feedback loop (NUAK1-YAP/TAZ-NUAK1) that continuously sustains fibrosis [24]. Additionally, NUAK1 coregulates reactive oxygen species (ROS) production and inflammatory responses with NF-κB signaling, suggesting that dual inhibition of NUAK1 and NF-κB may attenuate diabetic kidney disease (DKD) progression [25,26]. Recent investigations have identified NUAK1 as a key differentially expressed gene in DKD through bioinformatic analysis, with validation across multiple experimental systems—including human renal proximal tubular epithelial cells (HK-2), high-fat diet/streptozotocin-induced DKD mouse models, d-galactose-induced aging mouse models, and human peripheral blood mononuclear cells—confirming its upregulated expression. Further studies revealed that E26 transformation–specific 1 (ETS1) directly binds to the NUAK1 promoter, driving its transcriptional activation in DKD, and mediates renal tubular senescence and oxidative stress via the ROS/p53 axis, thereby expanding the known YAP/TAZ feedback pathway. Additionally, these studies identified Asiatic acid (AA) as a novel NUAK1 inhibitor through molecular docking and dynamics simulations; AA demonstrated stable binding affinity and was validated in both in vitro and in vivo experiments to effectively suppress NUAK1 expression and downstream pathological processes, significantly ameliorating renal injury in DKD models. The NUAK1-targeted strategy based on AA provides a promising direction for DKD drug development [27], and these findings are closely aligned with our own research results. Notably, our pseudotime trajectory analysis further revealed that NUAK1 expression peaks at late stages of PCT differentiation, suggesting that its pathological upregulation in DKD may specifically impair the terminal maturation of proximal tubule cells. This temporal specificity positions NUAK1 not merely as a fibrotic mediator, but as a stage-specific regulator of PCT cell fate, and provides a cellular context for understanding how UA-mediated suppression of NUAK1 may restore normal tubular differentiation programs [28].

PLK1 is defined as a widely expressed serine/threonine kinase in eukaryotic cells, with critical functions in mitosis, cell cycle progression, and DNA damage response [29,30,31]. PLK1 promotes renal tubulointerstitial fibrosis by modulating ATP6V1A phosphorylation and the autophagy/lysosomal pathway [32]. The PLK1/PTEN/PGK1 signaling axis represents a core mechanism linking metabolic reprogramming to fibrotic progression, with PLK1 inhibition attenuating extracellular matrix deposition [33]. Under high-glucose conditions, BRD4 recruits histone acetyltransferase P300 to enhance PLK1 transcription via H3K27 acetylation, subsequently activating NLRP3 inflammasome and triggering pyroptosis and fibrosis [34]. Upregulated PLK1 expression in diabetic mice and a significant positive correlation with serum creatinine in DKD patients were both observed, supporting the pathway relevance. Conversely, Wen et al. demonstrated that E3 ubiquitin ligase FBXO45 binds PLK1 and suppresses ferroptosis through the PLK1/GPX4/SOX2 axis, with PLK1 inhibition alleviating renal injury [35]. In contrast to NUAK1, pseudotime analysis in the present study identified PLK1 as predominantly expressed in early-stage PCT differentiation, implicating it in the regulation of proximal tubule progenitor cell states. The divergent temporal expression of PLK1 and NUAK1 along the PCT differentiation axis suggests that these two kinases govern non-overlapping checkpoints in tubular cell maturation [36], and that their concurrent dysregulation in DKD may disrupt the full spectrum of PCT developmental progression. Further investigation using larger clinical cohorts and mechanistic studies is warranted to elucidate PLK1’s precise contributions and pathway crosstalk.

HSPA1A encodes heat shock protein 70 (HSP70), a molecular chaperone functioning under stress conditions in protein folding, repair, and cell survival [37]. Its expression pattern in DKD remains controversial: clinical studies have associated elevated serum HSP70 levels with albuminuria and kidney injury extent in diabetic patients [38,39,40,41]. While some studies reported increased HSP70 expression in diabetic renal tissues linked to inflammatory responses [42,43], others suggested protective roles against cell death [43,44]. However, our data demonstrate HSPA1A downregulation in the DKD model, contrasting with earlier upregulation reports. This discrepancy may reflect the context-dependent duality of HSPA1A function: intracellular HSPA1A inhibits NF-κB signaling with anti-inflammatory effects, whereas extracellular HSPA1A promotes inflammation through TLR2/TLR4-mediated pathways [45,46]. Moreover, HSPA1A expression varies among renal cell types and its function shifts with disease stage and microenvironment [47]. Thus, the observed pattern likely reflects cell type-specific distribution and functional heterogeneity within the kidney. Further investigation integrating cell type specificity, disease progression, and subcellular localization is warranted to elucidate HSPA1A mechanisms in DKD.

ESR2 encodes estrogen receptor beta (ERβ), a nuclear receptor transcription factor implicated in various chronic diseases, including DKD [48]. Previous studies suggest ESR2 may confer renal protection by regulating lipid metabolism, upregulating PPARα, and inhibiting TGF-β-mediated fibrosis and extracellular matrix accumulation [49,50,51]. Based on the well-established renoprotective role of ESR2, a downregulation of ESR2 expression was anticipated in DKD. However, no significant changes in levels of ESR2 mRNA were detected in the renal cortex of db/db mice. These results suggest that ESR2 mediates renal protection mainly through post-translational mechanisms rather than transcriptional regulation. It is plausible that ESR2 activity is modulated via post-translational modifications (e.g., phosphorylation and acetylation) and influenced by ligand availability as well as specific cellular microenvironments. Moreover, bulk tissue analysis may miss changes restricted to specific cellular subpopulations. This limitation could explain the discrepancy between the observed data and initial expectations. Thus, the role of ESR2 in DKD likely operates predominantly at the level of protein activity, subcellular localization, and pathway crosstalk, rather than transcriptional regulation. Further investigation using single-cell sequencing and protein functional assays is warranted to elucidate the precise mechanisms of ESR2 across different cell types and disease stages in DKD.

Immune profiling revealed significant alterations in the renal immune microenvironment of DKD patients, characterized by dysregulation of immune cells including activated mast cells, CD8⁺ T cells, naïve B cells, Tregs, and Tfh cells, consistent with prior reports implicating these subsets in DKD pathogenesis [52,53,54]. Substantial infiltration of CD4⁺ and CD8⁺ T cells was observed within the diabetic renal interstitium, with enhanced production of inflammatory cytokines IFN-γ and TNF-α [55]. CD8⁺ T cells mediate renal injury through direct cytotoxicity and indirect inflammatory effects via proinflammatory cytokine secretion [56]. Adoptive transfer of Tregs can suppress CD8⁺ T cell infiltration, thereby ameliorating insulin resistance and attenuating DKD progression. Activated CD8⁺ T cells release effector molecules including granzyme B, TNF-α, and IFN-γ, exacerbating renal dysfunction and fibrosis [53]. Correlation analysis demonstrated significant positive associations between CD8⁺ T cell infiltration and expression of ESR2, NUAK1, and PLK1 (correlation coefficient > 0.7), whereas HSPA1A showed a notable negative correlation (correlation coefficient = -0.65). These genes were also correlated with T cells (Tfh) and naïve B cells, suggesting their involvement in reshaping the DKD immune microenvironment. This study delineates abnormal immune cell recruitment and gene interactions in DKD, providing theoretical foundation for immune-targeted therapeutic strategies.

Single-cell RNA sequencing further resolved 11 cell types, with proximal tubule epithelial cells (PCT) subdivided into three distinct subtypes: PTH1R⁺, HAO2⁺, and CLCNKA⁺ PCT. In an STZ-induced diabetic mouse model, expression levels of PTH1R (parathyroid hormone 1 receptor) and PTHrP (PTH-related protein) were significantly upregulated in both glomeruli and renal tubules, correlating with increased urinary albumin excretion. PTHrP appears to mediate diabetic renal cellular hypertrophy by regulating the cell cycle and TGF-β1 signaling pathway [57]. Although research on the HAO2⁺ PCT subtype in diabetic kidney disease (DKD) remains limited, HAO2 has been employed as a proximal tubule marker in numerous single-cell RNA sequencing studies for identifying cellular subpopulations. Notably, we identified a CLCNKA-expressing subpopulation within PCT, which is of particular interest given that CLCNKA has been consistently used as a canonical marker of the ascending thin limb and thick ascending limb in prior single-cell studies of renal fibrosis and injury, with clear demarcation from proximal tubule populations [58]. Its emergence within PCT in DKD therefore represents a transcriptional boundary crossing that is not observed under homeostatic conditions.This phenomenon is consistent with the broader concept of pathological phenotypic plasticity in renal tubular cells: single-nucleus sequencing studies of acute kidney injury have documented ectopic activation of segment-restricted genes within injured PTCs, interpreted as a stress-induced transcriptional reprogramming response rather than cell misclassification [59]. In the context of DKD, the appearance of CLCNKA⁺ PCT cells may similarly reflect a maladaptive ion transport reprogramming triggered by chronic hyperglycemia-induced metabolic stress, potentially contributing to tubular electrolyte dysregulation.This subpopulation has not been reported in prior DKD single-cell studies, and its functional characterization warrants dedicated experimental investigation [60]. Pseudotemporal trajectory analysis revealed stage-specific expression patterns: HSPA1A, ESR2, and NUAK1 were upregulated in late PCT development, whereas PLK1 peaked early, indicating temporal regulatory roles in tubular cell differentiation. These findings reinforce the central role of PCT cells in DKD, contributing to metabolic reprogramming, injury response, and fibroinflammatory progression [61,62,63]. In summary, this multi-level analysis—spanning immune contexture and single-cell transcriptomics—elucidates key cellular and molecular features of DKD, providing a foundation for targeting immune dysfunction and cell state-specific mechanisms in future therapeutic strategies.

Taken together, the findings of this study delineate a multi-level disease mechanism linking UA’s molecular targets to proximal tubule cell fate in DKD. At the transcriptomic level, NUAK1 and PLK1 were identified as UA-responsive kinases with stage-specific expression across the PCT differentiation trajectory — PLK1 governing early progenitor states and NUAK1 regulating late maturation — suggesting that their concurrent upregulation in DKD disrupts the full developmental continuum of PCT. In vivo and in vitro experiments confirmed that UA suppresses both kinases, restoring renal architecture and metabolic parameters.This mechanistic axis, whereby UA exerts renoprotection through suppression of NUAK1 and PLK1 — kinases with stage-specific roles in PCT differentiation — represents a framework that extends beyond the previously recognized antioxidant and anti-inflammatory actions of UA, offering a cell-state-level explanation for its renoprotective effects.These findings offer novel insights into DKD pathogenesis and suggest promising directions for therapeutic strategy development. Several limitations, however, should be considered. First, given the multifactorial nature of DKD, this work represents only an initial screening of candidate targets; extensive mechanistic investigations are required to elucidate their precise biological roles and molecular pathways. Second, the analyses relied on public datasets with limited sample sizes, which increases the risk of false discoveries and may affect the generalizability of the results. Although partial validation was achieved using independent cohorts and mRNA-level data, further functional assays and validation in large-scale, prospective clinical cohorts are essential. Addressing these aspects will be a primary objective of our ongoing research.

4. Conclusions

In conclusion, four hub genes were identified as diagnostic biomarkers for DKD through integrated multi-omics approaches, demonstrating potential for predicting disease progression. This study provides new insights into the underlying mechanisms and offers promising targets for the diagnosis and targeted therapy of DKD.

5. Materials and methods

5.1. Data collection

A total of three datasets relevant to DKD were retrieved from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Dataset GSE142025 (platform GPL20301), comprising renal tissues from 27 DKD patients and 9 controls, served as the training set. For validation, dataset GSE96804 (platform GPL17586), encompassing glomerular samples from 41 DKD patients and 20 controls was utilized. Additionally, the single-cell dataset GSE131882 (platform GPL24676) included three renal tissue samples from early DKD cases and three samples from healthy controls. The single-cell dataset GSE209781 (platform GPL24676) was also obtained, including three renal tissue samples from early DKD cases and three samples from healthy controls. The data were downloaded on May 28, 2025.

5.2. Differential expression analysis

Identification of DEGs in the training set was conducted with ‘DESeq2’ (v1.42.0) [64] using a P-adjust < 0.05 and |log₂FC| > 0.5. Visualization was achieved by generating a volcano plot with ‘ggplot2’ (v3.5.1) [65] and a heatmap with ‘ComplexHeatmap’ (v2.18.0) [66].

5.3. Identification and functional analysis of urolithin A-related DEGs

The canonical SMILES of UA was obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Following its acquisition, it was submitted to the SwissTargetPrediction database (https://www.swisstargetprediction.ch/) to screen for potential target proteins in Homo sapiens, applying a probability threshold of > 0. Next, the obtained target protein names were converted to target gene names via the UniProt database (https://www.uniprot.org/). Moreover, the UA-target gene network was visualised using the ‘Cytoscape’ software (v3.9.1) [67].

Identification of UA-related DEGs in DKD was performed by intersecting the UA target genes with DEGs using the ‘ggvenn’ package (v0.1.10) [68], with the overlapping genes defined as candidates. Subsequent Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted with the ‘clusterProfiler’ package (v4.10.1) [69] (P < 0.05) to display the most significantly enriched terms. A protein-protein interaction (PPI) network was generated by submitting candidate genes to Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/, confidence score > 0.4) followed by visualization in Cytoscape (v3.9.1).

5.4. Identification and validation of biomarkers

In the training set, feature genes were screened from the candidate pool using both least absolute shrinkage and selection operator (LASSO) regression (‘glmnet’ v4.1.8) [70] and support vector machine recursive feature elimination (SVM-RFE) (‘e1071’ v1.7.14) [71], and their intersection was obtained via a Venn diagram (‘ggvenn’ v0.1.10) [68]. Following this, the expression consistency of these feature genes was assessed across the training and validation sets using the Wilcoxon test (P < 0.05). Finally, genes exhibiting significant and concordant differential expression in both datasets were designated as candidate biomarkers.

The discriminatory power of candidate biomarkers for DKD versus control samples was evaluated across datasets using the ‘pROC’ package (v1.18.5) [72]. Receiver operating characteristic (ROC) curves were plotted and area under the curve (AUC) values computed in both training and validation sets, with genes achieving an AUC > 0.7 in both designated as biomarkers for subsequent analysis.

5.5. Gene set enrichment analysis (GSEA)

To explore the functional pathways in relation to biomarkers in DKD, the GSEA of biomarkers was performed. The reference gene set ‘c2.cp.v2023.2.Hs.symbols.gmt’ was downloaded from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb). For each biomarker, Spearman correlation coefficients with all other genes were computed across all training set samples using the ‘psych’ package (v2.4.3) [73]. Following this, genes were ranked in descending order of their correlation coefficients. This ranked list was subjected to GSEA using the ‘clusterProfiler’ package (v4.10.1) [69], applying thresholds of |NES| > 1, P.adjust < 0.25, and P < 0.05. Significantly enriched pathways were finally visualized employing the ‘enrichplot’ package (v1.22.0).

5.6. Immune infiltration analysis

We employed the CIBERSORT algorithm (v0.1.0) [74] to evaluate the infiltration patterns of 22 immune cell types in the training set [75]. Wilcoxon test identified cells with significant abundance differences between conditions (P < 0.05). Following this, Spearman correlations between differential immune cells and biomarkers were further analyzed using the ‘psych’ package (v2.4.3) [73], applying thresholds of |cor| > 0.3 and P < 0.05.

5.7. Molecular docking

In order to ascertain the binding affinity between the UA molecule and biomarkers, molecular docking simulations were carried out on UA and proteins encoded by the biomarkers. Initially, the chemical structure of UA in the SDF format was retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Subsequently, the 3D protein structures of proteins corresponding to the biomarkers were downloaded from the Protein Data Bank (PDB) database (https://www.rcsb.org/). Thereafter, CB-Dock2 (https://cadd.labshare.cn/cb-dock2/php/index.php) was utilized to compute the potential docking sites within the target proteins, followed by performing docking simulations and visualizing the results.

5.8. NephroSeq database analysis

Leveraging the NephroSeq V5 database (https://nephroseq.org/resource/login.html), we interrogated clinical patient data to assess correlations between biomarker expression levels and key renal function indices—serum creatinine and glomerular filtration rate (GFR), applying thresholds of |cor| > 0.3 and P < 0.05.

5.9. Single-cell data quality control and high-variable gene screening

To ensure the accuracy and reliability of subsequent analyses, data from the single-cell dataset were initially integrated using the ‘Seurat’ package (v5.1.0) [76] for ‘Harmony’ batch effect correction. Subsequently, a comprehensive quality control assessment was performed on the sample data. Given that mitochondria had already been filtered out in the raw data, the quality control criteria were established as follows: cells with cellular unique molecular identifier (UMI) counts (nCount_RNA) surpassing 16,000 were excluded, and cells with gene expression counts (nFeature_RNA) falling within the range of less than 400 or exceeding 4,000 were also removed.

A two-step procedure was employed to identify genes with a high intercellular coefficient of variation. First, data were normalized using the NormalizeData function in the ‘Seurat’ package (v5.1.0). Subsequently, the FindVariableFeatures function with the ‘vst’ method was applied to select these genes. The top 2,000 highly variable genes (HVGs) were displayed, and the top 10 most variable genes were highlighted using the LabelPoints function for visualization.

5.10. Cell dimension-reduction clustering, and annotation

For dimensionality reduction and visualization, we first scaled the data using the ScaleData function from ‘Seurat’ (v5.1.0). Principal component analysis (PCA) was then performed on the top 2,000 HVGs via RunPCA for clustering. Significant components (P < 0.05) were identified with JackStrawPlot, and the outcomes were visualized using ElbowPlot.

To determine the optimal number of cell clusters, we performed unsupervised clustering with the FindNeighbors and FindClusters functions, using a resolution of 0.6 and 10 dimensions. To visualize the clusters, we generated a uniform manifold approximation and projection (UMAP) plot via the RunUMAP function. To annotate the cell types, we identified marker genes from relevant published literature [12] and then generated a dot plot to display their expression levels across the different cell types.

5.11. Functional enrichment analysis and cell communication analysis

To explore the interactions among different cells, an analysis was conducted on ligand-receptor pairs as well as molecular interactions within the single-cell dataset, thereby presenting the intercellular communication networks and the receptor-ligand interaction probability point plots (P < 0.05).

To understand the biological functions in which differential cells were involved, functional enrichment analysis was performed for each cell using the analyze_sc_clusters function of the ‘ReactomeGSA’ package (v1.16.1) [77] (P < 0.05). Subsequently, the enrichment results corresponding to different cells were extracted via the pathways function.

5.12. Identification of key cells and pseudotiming analysis

Based on the proportion of each annotated cell type in the DKD and control groups, differences in cell abundance between samples were determined using the chi-square test. In this test, when the the P-value was < 0.05, cells that exhibited significant differences between different samples were defined as differential cells.

Based on biomarker expression levels in differential cells, cells with significant biomarker expression were selected as key cells (P < 0.05). The expression of the biomarkers was additionally confirmed using the GSE209781 dataset, which also underwent dimensionality reduction, clustering, and annotation of key cells. For trajectory inference, pseudotime analysis on these key cells was carried out with the ‘Monocle’ package (v2.30.1) [78]. We reduced the dimensionality of all cells in the key cell population using the DDRTree algorithm, then performed clustering and pseudotime trajectory analysis on dimensionally reduced key cells using the clusterCells function. Besides, density curves were drawn using the ‘ggpubr’ package (v0.6.0) [79] to reveal the quantitative changes of cells during development for each subset of cells. Finally, the plot_gene_in_pseudotime function was utilized to chart the expression trends of the biomarkers in cell subsets, focusing on their patterns across various developmental stages.

5.13. Animal experiments

Male C57BLKS-db/db mice (7 weeks old, 37.43 ± 2.00 g, n = 36) and db/m mice (23.59 ± 1.04 g, n = 12) were obtained from Changzhou Cavens Experimental Animal Co., Ltd. (SCXK(Su)2021-0013) and housed under specific pathogen-free conditions at Hebei Medical University. Animals were maintained at 22 ± 2 °C with a 12 h light/dark cycle. All mice were raised and fed under SPF conditions and all procedures were approved by the Institutional Animal Care and Use Committee of Hebei Medical University (IACUC-Hebmu-2025029). UA and Sodium carboxymethyl cellulose (CMC-Na) were purchased from MedChem Express (USA). After one week of acclimatization, hyperglycemic db/db mice (random blood glucose ≥ 16.7 mmol/L for two consecutive days) were randomly assigned to three groups (n = 12/group): db/db, db/db + 0.5% CMC-Na (vehicle, db/db+CMC), and db/db + UA (50 mg/kg/d, by gavage). Age-matched db/m mice served as normal controls (n = 12). All animals had ad libitum access to food and water. After 12 and 16 weeks, mice were euthanized under anesthesia. Kidney tissues were collected, snap-frozen in liquid nitrogen, and stored at -80 °C.

5.14. Biochemical assays

At designated time points (weeks 12 and 16), mice (n=6 per group) were fasted for 8 hours. Blood samples were then obtained from the retro-orbital plexus. After centrifugation (3000 rpm, 15 min, 4℃) to isolate serum, fasting blood glucose was measured using an Omron glucometer. Subsequently, serum creatinine (Scr), blood urea nitrogen (BUN), total cholesterol (T-CHO), triglycerides (TG), and albumin (ALB) levels were quantified using commercial kits (Nanjing Jiancheng Bioengineering Institute) as per the manufacturer’s instructions.

5.15. Electron microscopy

Immediately after collecting the renal specimens, the renal tissues were diced into 1-mm3 pieces and fixed in 4% glutaraldehyde (pH7.2–7.4) at 4℃. The renal tissues were dehydrated, infiltrated, and embedded to prepare ultrathin sections for observing the ultrastructural changes in the glomeruli and renal tubules using transmission electron microscopy.

5.16. Pathological and immunohistochemical staining

The kidney tissue from each mouse was fixed with 4% paraformaldehyde and embedded with paraffin wax to make 3-μm sections, which were stained with hematoxylin and eosin, HE, PAS, and Masson. Immunohistochemistry was performed with primary antibodies NUAK1 (1:100, ProteinTech), PLK1 (1:100, ProteinTech), which were incubated overnight at 4℃. The renal sections were subsequently incubated with biotin-conjugated secondary antibody and streptavidin coupled with horseradish peroxidase and counterstained with hematoxylin. The images were obtained using an inverted microscope (Olympus BX63). The quantitative analysis of positive staining was assessed by ImageJ software (NIH).

5.17. Western blot analysis

Total protein was extracted from the renal tissue using RIPA lysis buffer (Solarbio, Beijing, China). The protein concentration was determined using a BCA assay (Seven, China). Equal amounts of protein were separated by sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to polyvinylidene fluoride (PVDF) membranes (Millipore, USA). The membranes were then blocked with 5% skim milk at 37℃ for 1 h and incubated overnight with the primary antibodies NUAK1 (1:500), PLK1 (1:1000), GAPDH (1:5000) at 4℃. The following day, they were incubated with the secondary antibody (1:10000) at 37℃ for 1 h. The results were compared with those of GAPDH.

5.18. Cell culture

A human renal proximal tubular epithelial cell line, HK-2, was purchased from the American Type Culture Collection and cultured in Dulbecco’s Modified Eagle Medium/Nutrient Mixture F-12 (DMEM/F12) medium, supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin in 5% CO₂ conditions at 37℃. The HK-2 cells were divided into normal glucose concentration group (NG group, 5.5 mmol/L glucose), high glucose concentration group (HG group, 30 mmol/L glucose), high glucose concentration group plus UA group (UA group, 30 mmol/L glucose plus UA 10 μM). HK-2 cells were cultured in penicillin/streptomycin and serum-free medium for 6 h, then HG (30 mM), and HG plus UA (10 μM) media. The cells were cultured for 48 h and collected for analysis.

5.19. Cell viability assay

HK-2 cells were inoculated into 96-well plates at a density of 4 × 10³cells/well, each well containing 100 μL of the medium, and cultured at 37℃ in 5% CO₂ conditions. After the cells were adhered to the plate, they were treated with different concentrations of UA and 30 mM HG for 48 h. The cells were then incubated with 10 μL of CCK-8 reagent (HYK0301-500 T, MedChem Express) at 37℃ for 2 h. The optical density was read at 450 nm using an enzyme labeler (Bio-Tek, USA).

5.20. Immunofluorescence

HK-2 cells were incubated overnight with primary antibodies NUAK1 (1:100), PLK1 (1:100) at 4℃ and then incubated with goat anti-rabbit or goat anti-mouse IgG H&L secondary antibody at 37℃ for 2h. Finally, cells were stained with DAPI (SouthernBiotech, USA) for 5 min. Images were obtained using a laser scanning confocal microscope (Leica, Germany).

5.21. RNA extraction and quantitative PCR

Total RNA was extracted from kidney tissues using the RNAsimple Total RNA Kit (TIANGEN), according to the reagent vendor’s instructions. cDNA was synthesized using the FastKing RT Kit (With gDNase, TIANGEN). Quantitative PCR was performed using FastReal PreMix (SYBR Green, TIANGEN) on the Applied Biosystems QuantStudio 3&5 system (Thermo Fisher, USA). Primer sequences are listed in Table S9. GAPDH (Sangon Biotech, China) was used as the internal control, and relative gene expression was calculated using the 2^−ΔΔCT method.

5.22. Statistical analysis

Statistical analysis was performed using GraphPad Prism (v10.0) and R (v4.3.1). Measurement data are expressed as mean ± standard deviation (SD). For multiple group comparisons, one-way ANOVA or t-tests were applied based on data distribution. Differences between two groups were assessed using the Wilcoxon test. A p-value < 0.05 was considered statistically significant. All experiments were independently repeated.