The microarray datasets of glomeruli-associated DN samples and normal control glomeruli samples were systematically extracted from GEO (https://www.ncbi.nlm.nih.gov/geo/). The selection criteria included: 1) Homo sapiens; 2) expression profiling by array; 3) the experiment included patients with DN and controls; 4) availability of data for glomerulus. GSE30122 dataset consisting of 9 DN samples and 26 control samples was ultimately used as the training cohort, while GSE1009 and GSE104948 datasets were used as the test cohort to verify candidate genes. Meanwhile, a total of 866 cellular senescence-related genes (SRGs) were retrieved and downloaded from CellAge database ( Supplementary Table S1 ) (14). The details of these datasets are shown in Table 1 . The analytical flowchart of this study is illustrated in Figure 1 .

R statistical software was employed to normalize and analyze the raw microarray data. The analysis of differentially expressed genes (DEGs) between DN patients and controls were conducted by the “limma” package in R software with cut-off values of p< 0.05 and |log2FC|> 1. Probes without matching gene symbols were excluded, and only the probes with the highest signal values for the same molecule were retained. DEGs were visualized in volcano maps and heat maps by conducting R packages “ggplot2” and “ComplexHeatmap” respectively. Gene set enrichment analysis (GSEA) was performed on the DEGs to explore biological signaling pathways by R package “clusterProfiler”. The gene set “c5.go.all.v7.5.1.symbols.gmt” derived from the Molecular Signatures Database was chosen as reference. FDR <0.25 and the p. adjust value <0.05 were considered as statistically significant. Venn diagrams facilitated the screening intersection of DEGs and cellular senescence-associated genes, yielding the final DEGs associated with senescence (hereinafter referred to as DESRGs).

We deployed “ggplots”, “stats” and “car” packages in R software alongside the Mann Whitney U test to depicit the expressions of DESRGs between DN patients and controls. GO and KEGG enrichment analysis were processed via R package “clusterProfiler” with a threshold p-value of < 0.05 to determine the biological functions of the genes and relevant pathways. Through STRING online database (http://string-db.org), the PPI network was constructed with a combined score >0.4 to evaluate the protein-protein interactions of DESRGs (15). Software Cytoscape facilitated PPI network analysis and visualization.

The LASSO regression and Random Forest (RF) methods were adopted to determine hub genes from the 62 DESRGs via “glmnet” and “randomForest” packages respectively (16, 17). Marker genes were selected from variables corresponding to the value of the penalty parameter lambda.min by using 10-fold cross-validation in LASSO regression analysis. The RF model is an ensemble machine-learning algorithm for prediction (18). We set the number of decision trees to 500 and used parametric optimization to find the optimal parameters as “mtry=7”. Genes overlapping between the LASSO model and the top 20 genes from the RF model (with MeanDecreaseGini value >0.2) were defined as signature genes.

The receiver operating characteristic (ROC) curve was employed to evaluate the predictive accuracy of signature genes by the “pROC” R package in both training cohort (GSE30122) and test cohorts (GSE1009 and GSE104948). Visualization of results was implemented with the ggplot2 package and the genes with AUC > 0.7 were considered diagnostic.

Mice were randomly allocated into two groups: (i) a control group with normal blood glucose levels (Control group, n = 9-15), and (ii) streptozotocin (STZ)-induced diabetic mice (DN group, n = 9-15). Following a 12-hour fast, STZ (Sigma Chemicals, St. Louis, MO) was administered intraperitoneally (IP) at a dose of 50mg/kg, dissolved in 100mM, for five consecutive days. Diabetes was induced in mice aged 6-8 weeks, weighing 18-20g. The control group received citrate buffer (pH 4.5) for five consecutive days. Control mice received an equivalent volume of citrate buffer. After two weeks, blood was collected via submandibular vein puncture, and fasting blood glucose levels were measured using the HemoCue B-Glucose kit (HemoCue AB, Angelholm, Sweden). Mice with fasting blood glucose levels ≥ 12mmol/L were considered diabetic and included in further studies. After 20 weeks of STZ treatment, final body weight and kidney weight were measured. Blood samples were collected before euthanizing the mice under anesthesia. Some mice in various groups did not survive the induction process. All procedures were conducted following the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health, as well as the Animal Welfare Act guidelines.

Serum and urine creatinine values were quantified employing high-performance liquid chromatography (HPLC, Beckman Instruments, Fullerton, CA, USA). The urinary albumin value was determined using an immunoassay method (Bayer, Elkhart, IN, USA). The urinary albumin-to-creatinine ratio (UACR) was computed as the quotient of urine albumin and urine creatinine concentrations (μg/mg). All assays were executed according to the manufacturer’s protocols.

RNA was extracted from diabetic kidney tissue using the TRIzol reagent. The extracted RNA was converted into cDNA while removing possible genomic DNA contamination using the PrimeScript™ RT reagent Kit with gDNA Eraser. Real-time RT-PCR was performed using the Q7 RT-PCR Detection System and SYBR^® premix Ex Taq. Relative quantification (RQ) was calculated based on the 2^(-ΔΔCt) method, normalizing outcomes to β-actin mRNA as the control. The primer sequences used in the experiments are listed in Supplementary Table S2 . The overarching objective of this experiment was to gauge specific gene expression levels in diabetic kidney relative to various samples.

Diabetic kidney specimens were rinsed with pre-cooled PBS before homogenization. The homogenized tissues were then lysed in RIPA buffer (Beyotime, P0013F) containing protease inhibitor mixture (Roche Diagnostics, 4693116001), phosphatase inhibitor mixture (Roche Diagnostics, 4906845001), and PMSF (Beyotime, ST506). After a 30-minute ice bath incubation, the lysates were centrifuged at 13,000g for 15 minutes and the supernatant were stored at -80°C. Total protein concentration was estimated using the BCA Protein Assay Kit (Thermo Fisher Scientific, 23225). Protein samples (10-20μg) were diluted with 5x loading buffer, heated in a metal bath for 10-15 minutes, separated by 12.5% SDS polyacrylamide gel electrophoresis (SDS-PAGE), and then transferred to PVDF membranes (Millipore) for antibody probing. Immunoreactive proteins were detected using an enhanced chemiluminescence detection system (Millipore). We deployed a gamut of primary antibodies including Anti-LOX (CST, #58135), Anti-FOXD1 (Invitrogen, #PA5-35145), Anti-GJA1 (CST; #83649), Anti-BTG3 (Sangon Biotech; D220325) and Anti-GAPDH (KANGCHEN, KC-5G4). The secondary antibodies used were peroxidase-conjugated goat anti-mouse IgG (H+L) (33201ES60, Yeasen) and peroxidase-conjugated rabbit anti-goat IgG (H+L) (33701ES60, Yeasen). This Western blot procedure was conducted to discern specific protein expression levels in diabetic kidney samples, leveraging targeted antibodies for detection.

Eucommia ulmoides Oliver (EUO), one of the traditional Chinese medicine, is renowned for its effective treatment of kidney deficiency, lipid reduction and anti-obesity, attributable to the high concentrations of chlorogenic acids (CGA), quercetin, iridoid and α-linolenic acid in its leaves and bark (19, 20). CGA, also abundant in coffee beans, was reported to have substantial biological activities possibly mediating coffee’s beneficial impacts on glucose regulation and the development of type 2 diabetes (21). Zheng et al. presented that CGA have a positive effect on senescence delays and lifespan extension of Caenorhabditis elegans (22). However, a comprehensive review of the effects and mechanisms of CGA is currently unavailable. Therefore, our investigation sought to demystify the therapeutic potentials of CGA through the molecular docking.

The molecular structure of CGA was derived from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/), processed by the LigPrep module in Schrödinger software to generate the 3D conformations (23). AlphaFold2 (https://alphafold.ebi.ac.uk/) was utilized to predict the crystal structures of potential proteins (24). The obtained protein crystals were processed using the Protein Preparation Wizard of Schrödinger, including protein preprocess, regenerate states of native ligand, H-bond assignment optimization, protein energy minimization and remove waters. SiteMap and Receptor Grid Generation modules were applied to predict active sites of targeted proteins (25). We then performed the extra-precision docking (XP docking) of ligand compound CGA and active sites of the potential proteins. The lower XP Gscore represents the lower binding free energy between the compound and the protein, which implies higher binding stability (26). Molecular mechanics generalized Born surface area (MM-GBSA Analysis) was used to confirm the ligand-protein interaction and the binding free energy (27).

R software (version 4.2.1) and GraphPad Prism version 8.1 were applied in this study for all statistical analysis and graph visualization. A p-value of less than 0.05 was considered statistically significant.

A total of 756 DEGs were obtained from the comparisons between DN patients and control samples in the training cohort GSE30122, including 172 significantly upregulated and 584 significantly downregulated genes. Volcano plot and heatmap were utilized to depict these DEGs ( Figures 2A, B ). GSEA results ( Figure 2C ) revealed that the primary involvement of DEGs was in the immune responses processes. The top three correlations enriched depending on normalized enrichment score (NES) values were humoral immune response, activation of immune response, and positive regulation of immune response ( Figure 2D ).

We applied Venn diagram to determine 62 overlapped genes after intersecting the 756 DEGs with 866 cellular senescence-associated genes derived from the CellAge database ( Figure 3A ). Among these 62 DESRGs, 15 genes were upregulated and 47 genes were downregulated ( Figures 3B , 4A ). The correlations between DESRGs were demonstrated through Spearman correlation test in Figure 4B , suggesting all DESRGs interacted with each other. The color intensity and shape magnitude in Figure 4C highlighted the strength and direction of the relationships between each DESRG, with darker colors indicating stronger correlations.

To explore the biological functions and enrichment pathways, GO and KEGG analysis were performed. As shown in the results of GO and KEGG analysis ( Figures 5A, B ), the functional enrichments of 62 DESRGs were markedly involved in epithelial cell proliferation, reproductive structure and system development, regulation of epithelial cell proliferation, MAPK signaling pathway and Rap1 signaling pathway. The significantly enriched genes were revealed in Figure 5C . The PPI network was constructed through STRING database to further analyze the interactions between 62 DESRGs. A total of 48 nodes and 97 edges were identified ( Figure 5D ).

LASSO regression analysis was applied among all selected DESRGs to screen out 9 candidate genes, including DHCR24, GJA1, SFN, P2RY14, MARCKS, LOX, FOXD1, TPR and BTG3 ( Figures 6A, B ). Utilizing RF combined with feature selection, 20 genes were ranked based on their relative relevance with an estimated out-of-bag (OOB) error rate of 0% ( Figures 6C, D ). The intersection of the 9 candidate genes from the LASSO analysis with the top 20 feature genes from RF yielded four signature genes, namely BTG3, LOX, FOXD1 and GJA1 ( Figure 6E ).

ROC curves were operated in the training cohort to assess the diagnostic potential of signature genes. FOXD1, LOX, GJA1 and BTG3 showcased outstanding diagnostic values for differentiating DN patients from healthy individuals, with AUC values of 1.000, 1.000, 0.996 and 1.000 respectively ( Figures 7A–D ). We also created ROC curves in the test cohorts GSE1009 and GSE104948. In the GSE1009, these four characteristic genes all had significant diagnostic performance with AUC values of 0.889 in FOXD1, 1.000 in LOX, 1.000 in GJA1 and 1.000 in BTG3 ( Figures 7E–H ). In the GSE104948, LOX, GJA1 and BTG3 showed great predictive ability with respective AUC values of 0.937, 0.913, and 0.865, while FOXD1 had AUC of 0.468 ( Figures 7I–L ). These broadly similar findings in both training cohort and test cohorts suggested these four feature genes as potential diagnostic markers for DN.

To further validate our findings, we conducted experiments using a mouse model of DN. In comparison to the controls, the levels of LOX, FOXD1, and GJA1 were decreased while Kim-1 was significantly elevated in the DN mice. However, the expression of BTG3 in the DN mice did not show a significant change ( Figure 8A ). We then performed correlation analysis and observed that LOX, FOXD1, and GJA1 exhibited a negative correlation with Kim-1, while BTG3 showed no significant correlation with Kim-1 ( Figure 8B ). Subsequently, Western blot on diabetic kidney tissue confirmed these findings, indicating that the protein levels of LOX, FOXD1 and GJA1 also significantly decreased while Kim-1 in DN mice increased ( Figure 8C ).

We then assessed renal function in serum and urine samples from the mouse model by measuring serum creatinine (Scr) and the urine albumin-creatinine ratio (uACR) and performed the correlation analysis with the relative expressions of selected feature gens (Relative to GAPDH). Supplementary Figure S1A indicated that Scr and uACR in the DN group were significantly higher than in the control group, confirming the successful modeling. LOX and FOXD1 exhibited significantly negative correlations with Scr, while BTG3 and GJA1 demonstrated negative correlations without significance ( Supplementary Figure S1B ). LOX, FOXD1 and GJA1 showed significantly negative correlations with uACR while BTG3 had no significance ( Supplementary Figure S1C ). Supplementary Figure S1D depicted the correlation between the quantified protein expression levels of feature genes and Kim-1 ( Figure 8C ) detected in mouse kidney tissue. FOXD1 and GJA1 displayed negative correlations with the relative expressions of Kim-1. These findings were consistent with our initial analysis predictions, suggesting the potential of LOX, FOXD1 and GJA1 as indicators for DN in aging mice. Moreover, these genes might prospectively become therapeutic targets to manage elderly DN patients.

LOX, FOXD1 and GJA1, identified in the mice experiments, were subjected to molecular docking analysis alongside CGA to evaluate their potential therapeutic impacts on aging DN patients. The binding affinity and binding-free energy were gauged using XP Gscore (kcal/mol) and MM-GBSA dG Bind (kcal/mol) metrics respectively, to determine the ligands’ positions. The results of XP Gscore and MM-GBSA dG Bind between CGA and the three key targets were presented in Table 2 . LOX, FOXD1 and GJA1 were all demonstrated low binding free energy with MM-GBSA dG Bind values of lower than -30 kcal/mol, indicating the strong stability of ligand-protein binding. Among all three key proteins, GJA1 showcased the most stable and sufficient binding with CGA by its lowest XP Gscore of -7.045. FOXD1 followed with an XP Gscore of -6.273, demonstrating its fairly substantial binding stability.

Figures 9A, B depicted how CGA interacted with GJA1 by forming a hydrogen bond with SER86 and establishing hydrophobic interactions with surrounding residues including LEU90, PHE30, LEU26, VAL14, LEU11, etc. The interaction between CGA and FOXD1, where CGA created hydrogen bonds with ASN176 and LEU215, was illustrated in Figures 9C, D . CGA was additionally observed to interact with LOX as it formed hydrogen bonds with ASP188, GLY165, and ASP164, along with double hydrogen bonds with both ASP276 and ARG19, as shown in Figures 9E, F .