The microarray data of the mRNA expression profile related to DKD are retrieved from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). More details of the collected datasets are presented in Supplementary Table 1 . Cellular senescence-related pathways were derived from the Molecular Signatures Database (6). A total of 279 cellular senescence-associated genes were downloaded from the CellAge database https://genomics.senescence.info/cells/) (7). Clinical data for DKD patients were downloaded from the Nephroseq v5 online database (http://v5.nephroseq.org) (8). The mitochondrial genes and pathways were obtained from the MitoCarta3.0 database (http://www.broadinstitute.org/mitocarta) (9) and the Reactome database (https://reactome.org) (10). Our workflow is illustrated in Supplementary Figure 1.

To combine the data from the above 5 datasets, we first removed batch effects using the surrogate variable analysis (SVA) algorithm (11). The distribution patterns of samples were visualized by box plots and principal component analysis (PCA).

Kyoto Encyclopedia of Genes and Genomes (KEGG) (12) and Gene Ontology (GO) (13) enrichment analyses were applied using the R package clusterProfiler (14). Gene set enrichment analysis (GSEA) (15) was also performed to identify the underlying pathways; the threshold for significant terms was adjusted p-value < 0.05.

Gene set variation analysis (GSVA) is a nonparametric unsupervised analysis method mainly used to evaluate the gene set enrichment results of sequencing; GSVA allows the assessment of potential changes in pathway activity in each sample. The GSVA package in R software was used for the analysis, and the enrichment scores of pathways in all samples were calculated (16).

We used the Weighted Gene Co-Expression Network Analysis (WGCNA) algorithm (17) to screen for cellular senescence-associated module genes based on the enrichment scores of pathways obtained from GSVA. The “goodSamplesGenes” function was utilized to identify and remove outliers. To ensure that the co-expression network followed a scale-free distribution, we calculated a soft-thresholding power with the “pickSoftThreshold” function. We used the dynamic tree-cutting method to identify different modules, setting a minimum number of 100 genes per module. The module with the highest correlation with cellular senescence-associated pathways was identified as the senescence-associated module genes.

The differentially expressed genes (DEGs) between the normal and DKD groups were implemented via the “limma” R package (18) and visualized with heatmaps and volcano plots by the “ggplot2” package in R software (19). The intersection genes were obtained by the overlapping DEGs between control and DKD, cellular senescence-associated genes from the CellAge database, and module genes using a Venn diagram.

Two machine learning algorithms, Support Vector Machine-Recursive Feature Elimination (SVM-RFE) (20) and Random Forest (RF) algorithms (21), were employed to further screen the cellular senescence-associated signature genes in intersection genes. SVM-RFE is a sequence backward selection algorithm, which has superior classification performance for high-dimensional datasets (22). The SVM-RFE algorithm was implemented using the “e1071”, “kernlab”, and “caret” packages in R software for feature dimensionality reduction. The RF algorithm is an ensemble method that combines many decision trees and makes a single decision on behalf of the ensemble by combining the results of multiple classifiers together. The RF algorithm was implemented using the “randomForest” package in R software. Ultimately, genes overlapping among the machine learning algorithms were regarded as hub genes. We also developed a receiver operating characteristic (ROC) curve to assess the predictive capacity of the hub genes. The area under the ROC curve (AUC) value was calculated using the pROC package (23) to estimate the predictive utility of the hub genes.

The Least Absolute Shrinkage and Selection Operator (LASSO) (24) was used to construct optimal SRS in DKD. The risk score formula was as follows: ∑ i = 1 n β i ∗   exp ( i ) , where exp represented the gene expression value, and β represented the LASSO coefficient.

GSVA based on the single-sample gene set enrichment analysis (ssGSEA) algorithm was used to quantify the infiltration of 28 immune cell types (15).

All mouse studies were conducted according to protocols approved by the Ethics Committee of The Sixth Affiliated Hospital, Sun Yat-sen University Male C57BL/6 mice aged 4 weeks old were used for experiments. High-fat diet and streptozotocin (HFD/STZ) were used to produce the DKD model following the previous study (25). The mice were randomly divided into two groups: the control group and the HFD/STZ group. The control group was fed a normal chow diet, and the HFD/STZ group was given HFD. After 4 weeks, the mice in the HFD/STZ group were intraperitoneally injected with a dose of STZ (100 mg/kg) to induce diabetes. One week later, the mice were considered diabetic if their random blood glucose levels exceeded 16.7 mmol/L. The control group received citrate buffer and was processed in parallel with the diabetic mice. Mice were sacrificed 16 weeks later.

Paraffin-embedded kidney sections were stained using commercial kits (Servicebio, Wuhan, China) for hematoxylin and eosin (H&E) staining, periodic acid–Schiff (PAS) staining, and Masson’s trichrome staining.

At the termination of the animal experiment, we extracted total RNA from the kidney tissues of the mice using RNAiso Plus (Takara, Otsu, Japan). Then, reverse transcription was performed to synthesize cDNA (Novoprotein, Shanghai, China). The polymerase chain reaction was conducted using NovoStart SYBR qPCR SuperMix Plus (Novoprotein, Shanghai, China). The primer sequences are shown in Supplementary Table 2 .

Correlations between the 5 hub genes, the SRS risk score, and mitochondrial-related gene set were calculated using the Mantel test. The “cowplot” R package was used to create a graphical display of any correlations and their combinations.

All statistical tests were implemented utilizing R version 4.2.1 (https://www.r-project.org/) and GraphPad Prism 8.0. Wilcoxon or Student’s t-test was utilized for analyzing the difference between the two groups. The correlation between the variables was determined using Pearson’s or Spearman’s correlation test. All statistical p values were two-sided, and p < 0.05 was regarded as statistical significance.

After eliminating the batch effects by the SVA algorithm (11), we merged and normalized 5 DKD datasets. Box plots ( Supplementary Figures 2A, C ) and principal component analysis (PCA) ( Supplementary Figures 2B, D ) were used to visualize the distribution patterns between normal and DKD samples. The results confirmed that the batch effects had been successfully removed.

We observed a significant distinction between normal and DKD samples, with 687 DEGs identified, including 334 upregulated genes and 352 downregulated genes, as shown in the volcano plot and heatmap ( Figures 1A, B ). GSEA results revealed that the “KYNG_DNA_DAMAGE_UP”, “AGING_KIDNEY_UP”, and “FRIDMAN_SENESCENCE_UP” pathways were activated ( Figures 1C–E ), all of which are associated with cellular senescence. Furthermore, since the senescence process is typically accompanied by an inflammatory response (5), we also observed activation of the “HALLMARK_INFLAMMATORY_RESPONSE” pathway within the hallmark gene sets of the DKD group ( Figure 1F ). These findings indicate that the senescence process plays a role in the development of DKD.

Furthermore, we utilized GSVA to calculate the enrichment score of the aforementioned senescence-associated pathways to evaluate potential changes in pathway activity across each sample, revealing that senescence-associated pathway activity was elevated in the DKD group ( Figures 2A, B ). We then performed correlation analysis to validate the association between senescence-associated pathway activity and renal fibrotic markers, demonstrating a strong positive correlation between pathway activity and the expression of fibrotic genes, such as FN1, COL1A2, COL1A1, and ACTA2 ( Figure 2C ). Subsequently, we constructed a coexpression network using the enrichment score of senescence-associated pathways obtained via GSVA. We first calculated a soft threshold and established a scale-free topology model, setting the soft threshold to 5 ( Figure 2D ). The genes in the expression profile were then clustered ( Figure 2E ), with the yellow and turquoise modules displaying the strongest correlation with senescence-associated pathways ( Figure 2F ). Ultimately, we identified 2,938 module genes through WGCNA, which were identified as senescence-associated module genes in DKD.

To comprehensively characterize the expression pattern of cellular senescence-related genes, we utilized the CellAge database to obtain 279 human genes related to cellular senescence. We then intersected these genes with senescence-associated module genes and DEGs between control and DKD using a Venn diagram, which revealed 11 overlapping genes ( Figure 3A ). To further narrow down the range of key senescence-related genes, we subjected these 11 genes to machine learning analysis. SVM-RFE and RF models were independently established based on the control and DKD groups. The SVM-RFE algorithm identified 8 feature genes ( Figures 3B, C ), while the RF algorithm generated a sequence of the 11 genes ( Figures 3D, E ). Finally, we identified the 5 most important explanatory variables (LIMA1, ZFP36, FOS, IGFBP6, and CKB) from the machine learning models as hub genes ( Figure 3F ).

We further verified the expression of the 5 hub genes in the kidney of the normal control and DKD groups ( Figure 4A ). Additionally, we performed the ROC curve to evaluate the diagnostic efficacy of the 5 hub genes in DKD; the ROC curve indicated that all of the hub genes were highly accurate diagnostic markers for DKD (AUC between 0.8726 and 0.935) ( Figure 4B ). Subsequently, we utilized HFD/STZ to construct a DKD mouse model to verify the expression of the hub genes. Firstly, we validated the successful construction of the DKD model through histological analysis, which revealed that compared to the control group, the HFD/STZ group exhibited greater mesangial expansion and irregular thickening of the glomerular membrane. Additionally, Masson staining showed the formations of blue-stained extracellular collagen in the glomerulus ( Figure 4C ). Next, we performed RT-PCR to validate the expression of the hub genes in mouse kidney tissues. Our results indicated that the mRNA abundance of LIMA1 and IGFBP6 was higher in the HFD/STZ group, while that of ZFP36, FOS, and CKB was lower in the HFD/STZ group ( Figure 4D ), consistent with our bioinformatics analysis. Using the 5 hub genes, we employed LASSO to determine the number of factors by introducing shrinkage penalties and limiting the coefficients. Through continuous selection and simulation of the number of features, we established an SRS risk score for DKD patients ( Figures 4E, F ). The explicit formula of the SRS risk score was as follows: LIMA1 expression * (3.792075646) + ZFP36 expression * (−1.673662939) + FOS expression * (−0.150739822) + IGFBP6 expression * (2.350866436) + CKB expression * (−2.398663582). We have further validated the diagnostic efficiency of the SRS risk score in DKD patients and found that it is highly effective for diagnosis (AUC = 0.995) ( Figures 4G, H ).

We divided DKD patients into low-risk and high-risk groups based on the median of the SRS risk score. Subsequently, we examined the association between the risk score and renal function using clinical data from the Nephroseq v5 online tool. As depicted in Figure 5A , the high-risk group demonstrated a significantly lower glomerular filtration rate (GFR) than the low-risk group, and there was a negative correlation between the risk score and GFR levels ( Figure 5B ). We also investigated the link between the risk score and the expression of fibrotic genes. As illustrated in Figures 5C–F , the risk score exhibited a positive correlation with the expression of FN1, COL1A2, COL1A1, and ACTA2. These findings collectively validate the effectiveness of the SRS risk score as a means of assessing renal function and the degree of renal fibrosis.

In the low-risk and high-risk groups, a total of 541 DEGs were identified, with 215 showing downregulation and 326 showing upregulation patterns as demonstrated by the volcano plot and heatmap in Supplementary Figures 3A, B . These DEGs were primarily associated with biological processes related to senescence-associated secretory phenotype (SASP), which has been proposed as a possible source of inflammatory factors in DKD caused by cellular senescence. Such SASP-associated biological processes include “POSITIVE REGULATION OF CYTOKINE PRODUCTION”, “TUMOR NECROSIS FACTOR PRODUCTION”, “CHEMOKINE PRODUCTION”, and “TRANSFORMING GROWTH FACTOR BETA PRODUCTION” ( Figure 6A ). Furthermore, pathway enrichment analysis revealed that DEGs were enriched in key pathways implicated in DKD pathogenesis, such as the “AGE-RAGE SIGNALING PATHWAY IN DIABETIC COMPLICATIONS”, “RENIN-ANGIOTENSIN SYSTEM”, “MAPK SIGNALING PATHWAY”, “PI3K-AKT SIGNALING PATHWAY”, and “TGF-beta SIGNALING PATHWAY”, among others ( Figures 6B–D ). We also conducted GSEA to identify underlying pathways or processes in hallmark gene sets obtained from the Molecular Signatures Database for patients with low or high scores. The results indicated that “EPITHELIAL MESENCHYMAL TRANSITION” and “TNF-α SIGNALING VIA NF-κb” were activated in the high-score group, while the activities of mitochondrial pathways, such as “FATTY ACID METABOLISM” and “OXIDATIVE PHOSPHORYLATION”, were suppressed ( Figure 6E ).

The previous analysis suggests that pathways related to mitochondrial function such as oxidative phosphorylation (OXPHOS) and fatty acid metabolism were significantly suppressed in the high SRS score group. Furthermore, mitochondrial dysfunction is a key hallmark of cellular senescence (26). Therefore, we further investigated the relationship between mitochondrial function and cellular senescence in our subsequent research. We obtained mitochondrial genes and pathways from the MitoCarta3.0 and Reactome database and assessed the differential expression of mitochondrial genes between the low- and high-score groups. Our analysis revealed that almost all differentially expressed mitochondrial genes were downregulated in the high SRS score group, primarily localized to the mitochondrial matrix and mitochondrial inner membrane (MIM) ( Figure 7A ). Based on the differential gene analysis between low- and high-score groups, we conducted GSEA enrichment analysis on mitochondrial pathways obtained from the MitoCarta3.0 database, indicating a significant suppression of pathways involved in OXPHOS and metabolism in the high SRS score group ( Figure 7B ). To further explore differences in mitochondrial pathways between the high- and low-score groups, we conducted GSVA analysis on mitochondrial pathways. We selected a log fold change threshold of 0.2 and a p-value below 0.05 to identify 41 pathways with significant differences, among which 39 pathways were downregulated in the high-score group ( Figure 7C ). Most of the downregulated pathways are related to OXPHOS and metabolism. In addition, the Mantel t-test confirmed an extremely strong correlation between the SRS score and the 5 hub genes with the OXPHOS and carbohydrate/lipid/amino acid metabolism gene sets ( Supplementary Figures 4A–D ). These findings suggest that mitochondrial function was widely suppressed in the high-score group.

To further investigate the mechanism behind the decreased activity of OXPHOS-related pathways in the high-scoring group, we analyzed the relationships between the expression of the SRS risk score, 5 hub genes, and genes related to the mitochondrial respiratory chain complex. As genes related to mitochondrial respiratory chain complex V were undetected in the expression profiles, we focused solely on genes related to complex I–IV. The correlation values revealed a clear negative correlation between the most respiratory chain complex I–IV genes and the SRS risk score, LIMA1, and IGFBP6 expression levels, while ZFP36, FOS, and CKB showed a positive correlation ( Figure 7D ).

Continuing our investigation, we delved deeper into the correlation between the SRS score, 5 hub genes, and mitochondrial metabolism genes. Using the MitoCarta3.0 database, we obtained genes involved in mitochondrial metabolic pathways associated with various nutrient metabolisms. Our correlation analysis revealed a consistent regulatory effect of the SRS score and the 5 hub genes on genes regulating mitochondrial metabolism-related genes including macronutrients (carbohydrates/lipids/amino acids) and micronutrients (vitamins and minerals). Our findings indicated that the SRS risk score, LIMA1, and IGFBP6 expression levels showed a negative correlation with most mitochondrial metabolism-related genes, while ZFP36, FOS, and CKB showed a positive correlation ( Figure 7E ).

Previous studies have suggested that both mitophagy and mitochondrial biogenesis play important roles in the senescence process (27, 28). Our GSVA indicated that the activity of the mitophagy pathway was decreased in the high-risk score group ( Figure 7C ). We downloaded gene sets related to PINK1/PRKN-mediated mitophagy, receptor-mediated mitophagy, and mitochondrial biogenesis from the Reactome database. Unfortunately, we did not find any clear patterns or correlations between the SRS risk score and the 5 hub genes with the genes related to the mitophagy gene sets ( Supplementary Figures 4E, F ). Regarding the gene sets related to mitochondrial biogenesis, we found that only SIRT3 was negatively correlated with both the SRS risk score and the expression levels of LIMA1, ZFP36, FOS, and IGFBP6. Meanwhile, CKB showed a positive correlation with SIRT3 expression. However, there was no significant correlation between the SRS risk score, the 5 hub genes, and other genes related to mitochondrial biogenesis as shown in Supplementary Figure 4G .

One of the hallmarks of DKD is immune remodeling (29). To determine whether the SRS risk score accurately reflected the immune status of DKD, we evaluated immune cell infiltration in DKD using ssGSEA. We observed distinct immune infiltrate patterns among patient risk groups stratified by the SRS risk score. Compared to the low-risk group, most innate and adaptive immune cells showed higher levels of infiltration in the high-risk group ( Figure 8A ). Furthermore, there were notable interactions between immune cell populations across kidney tissues affected by DKD ( Figure 8B ). Correlation analysis revealed that infiltration levels of multiple types of immune cells, including T-cell subsets, different developmental or functional stages of B cells, and mast cells, were positively correlated with the SRS risk score, while CD56dim natural killer cells, monocytes, and CD56bright natural killer cells were negatively correlated with the risk score ( Figure 8C ). These findings suggest that immune system disorders in the process of DKD may be closely related to cellular senescence.