Identification and analysis of cellular senescence-associated signatures in diabetic kidney disease by integrated bioinformatics analysis and machine learning

Background Diabetic kidney disease (DKD) is a common complication of diabetes that is clinically characterized by progressive albuminuria due to glomerular destruction. The etiology of DKD is multifactorial, and numerous studies have demonstrated that cellular senescence plays a significant role in its pathogenesis, but the specific mechanism requires further investigation. Methods This study utilized 5 datasets comprising 144 renal samples from the Gene Expression Omnibus (GEO) database. We obtained cellular senescence-related pathways from the Molecular Signatures Database and evaluated the activity of senescence pathways in DKD patients using the Gene Set Enrichment Analysis (GSEA) algorithm. Furthermore, we identified module genes related to cellular senescence pathways through Weighted Gene Co-Expression Network Analysis (WGCNA) algorithm and used machine learning algorithms to screen for hub genes related to senescence. Subsequently, we constructed a cellular senescence-related signature (SRS) risk score based on hub genes using the Least Absolute Shrinkage and Selection Operator (LASSO), and verified mRNA levels of hub genes by RT-PCR in vivo. Finally, we validated the relationship between the SRS risk score and kidney function, as well as their association with mitochondrial function and immune infiltration. Results The activity of cellular senescence-related pathways was found to be elevated among DKD patients. Based on 5 hub genes (LIMA1, ZFP36, FOS, IGFBP6, CKB), a cellular senescence-related signature (SRS) was constructed and validated as a risk factor for renal function decline in DKD patients. Notably, patients with high SRS risk scores exhibited extensive inhibition of mitochondrial pathways and upregulation of immune cell infiltration. Conclusion Collectively, our findings demonstrated that cellular senescence is involved in the process of DKD, providing a novel strategy for treating DKD.


Introduction
Diabetic kidney disease (DKD) is one of the most common complications of diabetes and remains the leading cause of endstage renal disease (ESRD) in the world population (1). The current treatments for DKD are largely confined to optimal glucose and blood pressure control, as well as renin-angiotensin-aldosterone system blockades, which have limited therapeutic efficacy (1). It is therefore imperative to gain a deeper understanding of the pathogenesis of DKD, identify reliable biomarkers for high-risk patients, and develop novel therapeutic approaches to prevent or reverse the progression of DKD.
Cellular senescence is a cellular state that involves a halt in proliferation, resulting in inflammation, impaired tissue repair, irreversible tissue damage, and organ dysfunction (2). Recent research has shown that renal parenchymal cells are induced to undergo cellular senescence in the context of diabetes, by various pathogenic stimuli, including hyperglycemia (3) and the accumulation of advanced glycation end products (AGEs) (4). This leads to the deterioration of kidney function mainly by mediating a complex pro-inflammatory response called senescence-associated secretory phenotype (SASP) (3). The SASP includes the secretion of various molecules including cytokines, chemokines, and growth factors, and has been suggested as a possible source of inflammatory factors in DKD (5). However, the linkage of cellular senescence and DKD is complex and not yet fully understood at present. Therefore, it is vital to clarify the role and mechanism of cellular senescence in the pathogenesis of DKD.
To systematically assess the correlations between cellular senescence and the pathogenesis of DKD, we developed a cellular senescence-related signature (SRS) risk score using biomarkers that are closely associated with cellular senescence, and evaluated its potential significance in predicting renal function. Subsequently, we grouped patients according to the SRS risk score and compared mitochondrial function and immune cell infiltration between the high and low the SRS score groups. Our study sheds new light on the regulatory mechanisms of cellular senescence in the process of DKD.

Data acquisition
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 senescenceassociated 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.

Data preprocessing
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).

Pathway and functional enrichment analysis
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
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).

Construction of the co-expression network and key module identification by weighted gene co-expression network analysis
We used the Weighted Gene Co-Expression Network Analysis (WGCNA) algorithm (17) to screen for cellular senescenceassociated 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.

Identification of intersection 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 senescenceassociated genes from the CellAge database, and module genes using a Venn diagram.

Screening hub genes of DKD by machine learning
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.

Development and validation of the SRS based on the LASSO
The Least Absolute Shrinkage and Selection Operator (LASSO) (24) was used to construct optimal SRS in DKD. The risk score formula was as follows: o n i=1 bi * exp(i), where exp represented the gene expression value, and b represented the LASSO coefficient.

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

Animal experimental design
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.

RT-PCR
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.

Correlation analysis between two matrices
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.

Statistical analysis
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.

Removal of batch effects
After eliminating the batch effects by the SVA algorithm (11)

Identification of DEGs and enrichment analysis between control and DKD
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.

Analysis of module closely related to cellular senescence in 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 senescenceassociated pathways ( Figure 2F). Ultimately, we identified 2,938 module genes through WGCNA, which were identified as senescence-associated module genes in DKD.

Hub genes identified by machine learning models
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).

Construction of SRS risk score based on the 5 hub genes
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 bluestained 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).

Validation of the role of the SRS risk score in DKD
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 lowrisk 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. 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 senescenceassociated 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-

R A G E S I G N A L I N G P A T H W A Y I N D I A B E T I C 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-a SIGNALING VIA NF-kb" were activated in the high-score group, while the activities of mitochondrial pathways, such as "FATTY ACID METABOLISM" a n

d " O X I D A T I V E P H O S P H O R Y L A T I O N "
, were suppressed ( Figure 6E).

Evaluation of mitochondrial pathway
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 The relationship between the SRS risk score and GFR, fibrotic markers. (A) The box plots of GFR between the low-risk and high-risk groups. (B) The correlation analysis between the SRS risk score and GFR. (C-F) The correlation analysis between the SRS risk score and the expression of fibrotic markers (FN1, COL1A2, COL1A1, and ACTA2). *** p < 0.001.
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 highscore 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). Assessment of mitochondrial pathways in the low-risk and high-risk groups. (A) Heatmap of differentially expressed mitochondrial genes in the lowrisk and high-risk groups. (B) GSEA of mitochondria-related pathways in the low-risk and high-risk groups. (C) GSVA of mitochondrial pathways demonstrated by box plots in the low-risk and high-risk groups. (D) Correlations between the SRS risk score, the 5 hub genes and mitochondrial respiratory chain complex-related genes. (E) Correlations between the SRS risk score, the 5 hub genes and mitochondrial metabolism-related genes. * p < 0.05, ** p < 0.01, *** p < 0.001.
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

Evaluation of immune cell infiltration
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.

Discussion
Our study constructed an SRS risk score using 5 biomarkers related to cellular senescence based on the GEO datasets. We found that patients in the high-scoring group had lower GFR and higher expression of fibrotic genes compared to those in the low-scoring group. Additionally, we observed significant suppression of mitochondrial function in the high-scoring group. Notably, the SRS risk score was positively correlated with the degree of immune cell infiltration in DKD. These results provide new insight into the role of cellular senescence and associated genes in the pathogenesis of DKD. Cellular senescence was first reported in the 1960s (30) and is defined as a persistent cell cycle arrest that limits their proliferative life span. Emerging evidence for accelerated senescence and SASP has been reported in DKD (3,31). Recent studies have shown that drugs targeting senescent cells can alleviate proteinuria in obese insulin-resistant mice (32), suggesting that anti-cellular senescence strategies hold promise for future interventions aimed at delaying, preventing, or treating renal dysfunction. In our research, the 5 hub genes (LIMA1, IGFBP6, ZFP36, FOS, and CKB) related to cellular senescence were identified. Importantly, these hub genes have significant implications in multiple senescence-related fields. For example, ZFP36 contributes to inflammation-associated lung damage by interacting with the p53/p21 pathway, which is a core pathway in senescence (33). In addition, FOS plays a crucial role in the senescence process of granulosa cells in the ovary (34), and is responsible for UV-related skin aging (35). CKB is involved in cigarette smoke-induced bronchial epithelial cell senescence (36). However, the roles of these hub genes in the regulation of cellular senescence in DKD are less well understood. It should be noted that further research is necessary to validate the effects of these identified genes in DKD.
Clear evidence indicates that AGEs formation generation plays a central role in DKD pathology mechanisms. The binding of AGEs with receptor for AGEs (RAGE) provokes oxidative stress and chronic inflammation in renal tissues, resulting in progressive renal diseases (37). These processes involve activations of signaling pathways nuclear factor-kappa B (NF-kB), PI3K/Akt, and MAPK/ERK (38). Our study found that the AGE/RAGE, PI3K/Akt, and ERK pathways were enriched in the DEGs grouped by the SRS risk score; these results suggest that cellular senescence is involved in the regulation of the core pathways in DKD.
Along with cellular senescence, mitochondrial dysfunction is an essential "hallmark" and drives and maintains cellular senescence (26,39). However, there is still no consensus on which type of mitochondrial dysfunction is the primary driving factor of cellular senescence. Previous studies have emphasized that mitochondrial OXPHOS dysfunction is a common feature of senescent cells (40), and our study confirms this finding again in the DKD model. Furthermore, we have demonstrated that genes related to mitochondrial respiratory chain complexes are widely downregulated in patients with higher SRS score. Mitochondria integrate nutrients as fuels to generate energy for the cell (41), and the metabolic disorders of these nutrients are involved in the occurrence and development of DKD (42-44); however, few studies have suggested that nutrient metabolism disorder is closely related to cellular senescence in the process of DKD, and it is worth investigating whether targeting mitochondrial metabolism could potentially restore mitochondrial function during cellular senescence. Additionally, numerous studies have suggested that mitophagy plays a role in regulating cellular senescence and is involved in the development of DKD (31). Our GSVA indicated that the activity of the mitophagy pathway was decreased in the high SRS score group. Unfortunately, we did not find any clear patterns or correlations between our SRS risk score and the 5 hub genes with the genes related to the mitophagy gene sets. This may be due to the fact that the regulation of mitophagy often occurs at the post-transcriptional protein modification level (45, 46).
Although immune dysregulation has been recognized as a crucial driving factor in the pathogenesis of DKD, the potential regulation of the immune system remains largely unknown. Our study identified numerous immune-related biological processes and pathways that were significantly enriched in the differentially expressed DEGs grouped by the SRS risk score. Notably, most immune cells were increased in the high-risk score group, suggesting that cellular senescence may trigger an immunological response in DKD. Multiple studies have reported that immune cell infiltration contributes to the development of DKD. For instance, group 2 innate lymphoid cells have been implicated in the pathogenesis of renal fibrosis in DKD (47). The deposition of macrophages in the kidney is closely associated with decreased renal function in DKD patients (48), while mast cells are reported to participate in renal interstitial fibrosis during DKD progression (49). Memory CD8+ T cells were also found to be significantly elevated in kidney tissues in various types of chronic kidney disease (CKD), including DKD, leading to podocyte injury and glomerulosclerosis (50). Based on our findings, targeting senescent cells could potentially represent a novel therapeutic approach for improving immune dysfunction in DKD patients.
This study has several limitations. First, we built and evaluated our SRS from public databases. To confirm its clinical value, more prospective real-world evidence is needed. Second, the SRS construction based solely on a single signature is inescapable because many significant genes in DKD may have been ignored. Third, the association between SRS and mitochondrial pathway, immune response, needs to be studied experimentally.
In conclusion, an innovative SRS can distinguish normal control and DKD patients and can be utilized to predict kidney function. Moreover, the SRS was found to be linked to mitochondrial function and immune cell infiltration.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement
The animal study was reviewed and approved by The Sixth Affiliated Hospital, Sun Yat-sen University.

Author contributions
TZ designed and supervised the entire work. YL conceived the study and performed the literature search and bioinformatics analysis. LZ performed the experiments and wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding
This work was supported by grants from the National Natural Science Foundation of China (NSFC 82270872) and the Natural Science Foundation of Guangdong (2022A1515012249).