Integrated Analysis of Ferroptosis and Immunity-Related Genes Associated with Diabetic Kidney Disease

Introduction

Chronic diseases such as diabetes are becoming increasingly prevalent each year owing to the steadily growing older population, which poses a significant burden to the public healthcare system.¹ Diabetic kidney disease (DKD) is a serious microcirculatory complication in the late stages of diabetes and is the leading cause of chronic kidney disease.² According to reports, 35–40% of diabetic patients will experience irreparable kidney impairment and advance to DKD.³ The precise pathophysiology of DKD is currently unknown, and there is a lack of effective therapies. Therefore, exploring the mechanisms underlying DKD and identifying drugs to treat these conditions have become the therapeutic priorities.

Ferroptosis is a type of cell death caused by iron-dependent lipid peroxidation.⁴ In DKD, ferroptosis produces large amounts of reactive oxygen species, causing kidney damage through oxidative stress.⁵ Inhibiting ferroptosis could be a strategy to treat DKD.⁶ Ferroptosis offers a new perspective on the pathogenesis of DKD; it is a complicated physiological process regulated by multiple genes and its specific regulatory mechanisms require further investigation.

There is a relationship between ferroptosis and immune cell infiltration. Ferroptotic cells can initiate macrophage recruitment⁷ and affect macrophage polarization via iron metabolism.⁸ Immune cells undergo phagocytosis and iron accumulation. Activated immune cells secrete certain cytokines⁹ and ROS¹⁰ to regulate ferroptosis. Li et al¹¹ and Ma et al¹² suggested a possible connection between immunity and ferroptosis in DKD; however, the precise mechanisms remain unknown. The role of ferroptosis in immunotherapy has received increasing attention.¹³ Altering macrophage polarization could help treat cancer by altering the tumor microenvironment.^14–16 Therefore, further studies on reliable biomarkers connected to the ferroptosis-immune crosstalk in DKD may provide new perspectives on the pathophysiology of the disease and potential treatment options.

This study aimed to perform an integrated analysis of ferroptosis and immune-related differentially expressed genes (FIRDEGs) in DKD. In this study, we searched for differentially expressed genes (DEGs) in the GSE30122, GSE30528, and GSE30529 datasets to investigate the primary pathogenesis of DKD. An integrated analysis of ferroptosis and immune-related DEGs was performed to observe immune cell infiltration and its relationship with ferroptosis. In addition, the expression, function, and regulatory elements of FIRDEGs were observed and investigated. The RNA expression of FIRDEGs was validated in animal kidney samples using reverse transcription-polymerase chain reaction (RT-PCR). The finding offers a possible therapeutic target for treating patients with DKD and provides new information on the mechanism of DKD cell death.

Materials and Methods

Data Acquisition

The gene expression profiles of GSE30122,¹⁷ GSE30528,¹⁷ and GSE30529¹⁷ were downloaded from the Gene Expression Omnibus (GEO)¹⁸ based on the GPL571[HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array platform and annotated using the R package GEOquery.¹⁹ Nineteen kidney samples from patients with DKD and 50 from normal controls in GSE30122, 9 DKD samples and 13 normal control samples in GSE30528, and 10 DKD samples and 12 normal control samples in GSE30529 were included in the study.

Ferroptosis-related genes (FRGs) were obtained from three databases: GeneCards (http://www.genecards.org/), MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/),²⁰ and FerrDb (http://www.zhounan.org/ferrdb/).²¹ Immune-related genes (IRGs) were obtained from GeneCards and MSigDB databases.

Analysis of DEGs

To identify the pathways and biological mechanisms related to DKD, we used the limma package²² to find DEGs in GSE30122, GSE30528, and GSE30529. The threshold was defined as a P-value < 0.05, and |log2FoldChange| > 0.²³ In this study, DEGs with logFoldChange (FC) > 0 and P < 0.05 were considered upregulated genes, whereas those with logFoldChange < 0 and P < 0.05 were considered downregulated genes.

To analyze the FIRDEGs correlated with DKD, co-DEGs and ferroptosis immune-related genes (FIRGs) were analyzed together and the results were plotted as a Venn diagram. The R package pheatmap was used to create the heat map, and ggplot2 was used to draw the volcano map.

Gene Ontology and Pathway Enrichment Analysis

The packages “cluster Profiler”²⁴ in R software were used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for FIRDEGs. P-values were adjusted using the Benjamini–Hochberg (BH) method. P < 0.05 was used as the threshold for statistical significance, and a false discovery rate (FDR) < 0.05 was considered significant.

Distribution of a predefined gene set was evaluated using GSEA to ascertain the contribution of each gene to the phenotype. In this study, the logFC value was used for molecular ordering to evaluate significant enrichment in the predefined gene set, and the clusterProfiler package was used to conduct enrichment analysis on all genes related to the phenotype. The parameters used in GSEA were as follows: seed, 2020; number of calculations, 1000; number of genes contained in each gene set, at least 10; and number of genes contained, at most 500. P-values were adjusted using the BH method. GSEA of the GSE30122 dataset was performed using gene set collections (h.all.v7.4. symbols.gmt) from MSigDB. If P < 0.05 and FDR <0.25, enrichment was regarded as significant.

Protein-Protein Interaction (PPI) Network

PPI networks were investigated using the STRING database²⁵ and visualized using Cytoscape software. CytoHubba²⁶ of Cytoscape was used to score each node gene using five algorithms: Degree Correlation (DEGREE), Maximum Neighborhood Component (MNC), Maximal Clique Centrality (MCC), Edge Percolated Component (EPC), and Density of Maximum Neighborhood Component (DMNC). The top 15 node genes as scored by the algorithms were selected as hub genes.

Multi-Factor Regulatory Network

Interactions between target genes and miRNAs were predicted using two databases: ENCORI²⁷ (https://starbase.sysu.edu.cn/) and miRDB²⁸ (http://mirdb.org/). We used the CHIPBase database²⁹ (https://rna.sysu.edu.cn/chipbase/) and hTFtarget database³⁰ (http://bioinfo.life.hust.edu.cn/hTFtarget) for predicting the associations between transcription factors (TF) and hub genes. The DGIdb database³¹ (http://www.dgidb.org) was used to identify potential metabolites or compounds associated with hub genes. Cytoscape software was used to visualize mRNA-miRNA, mRNA-TF, and mRNA-drug interaction networks.

Evaluation of Immune Cell Infiltration

Gene expression profile data were uploaded to the CIBERSORT website³² (https://cibersort.stanford.edu) to obtain an immune cell infiltration matrix. Differences in immune cell infiltration between normal and DKD groups were assessed using the Wilcoxon test. The correlation between different immune cells was calculated using Spearman’s analysis and visualized using the R package ggplot2. The R package ggplot2 was used to visualize the heat map of the correlation between FIRDEGs and infiltrating immune cells.

Correlation Analysis with Renal Function

The association between FIRDEGs expression and clinical characteristics was validated using the NephroSeq v5 online database (http://v5.nephroseq.org) using Pearson’s correlation analysis.

Receiver Operating Characteristic (ROC) Curve

To identify the underlying mechanisms and related biological characteristics of hub genes in DKD, the correlations between hub genes were displayed using a group comparison map and a correlation heatmap for the GSE30122, GSE30528, and GSE30529 datasets. The R package pROC was used to visualize the ROC curves of hub genes between the normal and DKD groups. The area under the curve (AUC) was calculated to assess the diagnostic effect of hub genes on DKD.

Validation of FIRDEGs in a Mouse Model

Animal Model

Animal models of type 2 diabetes (db/db male mice (n = 6)) and db/m male mice (n = 6) were used as controls. All the mice were fed until they were 20-week-old. Metabolic cages were used to collect 24-hour urine from mice. Glucose meters (Roche Diagnostics, Shanghai, China) were used to monitor fasting blood glucose levels. After the mice were anesthetized, blood was collected through the orbit and centrifuged to preserve the upper serum. After sacrifice, kidney and body weight were recorded. Kidney tissues were collected following cardiac perfusion and frozen at −80°C. The animal experiments were carried out in line with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and were granted permission by the Anhui Medical University Animal Ethics Committee (LLSC20221089).

Chemicals and Reagents

Hematoxylin and eosin (H&E) and periodic acid–Schiff (PAS) staining kits were acquired from Solarbio (Beijing, China). Immunohistochemistry (IHC) kits were obtained from Beijing Zhongshan Company (Beijing, China). The PCR amplification kit was acquired from Shanghai Yeasen Company (Shanghai, China). Mouse monoclonal anti-Gpx4 antibody (1:100, cat number:67763-1-Ig) was provided by Proteintech (Wuhan, China). Mouse monoclonal anti-Acsl4 antibody (1:100, cat number: sc-365230) was provided by Santa Cruz Biotechnology (CA, USA). Mouse urine albumin, serum creatinine, and serum urea nitrogen detection kits were purchased from the Nanjing Jiancheng Bioengineering Institute (Nanjing, China).

Pathological Staining of Kidney Tissues

H&E and PAS staining were performed followed the manufacturer’s methods. IHC assay was performed as previously described.³³ After dewaxing the kidney tissues, the antigens were extracted by boiling in citrate buffer. The sections were blocked with peroxidase and 10% normal goat serum, and then incubated with primary antibodies at 4 °C overnight. They were then incubated with secondary antibodies for 1 h at 37 °C; this was followed by chromogenic detection with 3, 3’-diaminobenzidine.

RT-PCR Assay

Total RNA was extracted from the kidney tissue of mice following the manufacturer’s instructions. Reverse transcription and RT-PCR were performed using the Hifair II First Strand cDNA Synthesis Kit and SYBR Green Kit (Yeason Biotech, Shanghai, China), respectively. The RT-PCR primer sequences are listed in Table 1.

Table 1 Primer Sequences of RT-PCR

Statistical Analysis

All statistical analyses were performed using the R package (V4.1.2). Comparisons involving continuous variables between two groups were performed using Student’s t-test for normally distributed variables and the Mann–Whitney U-test for non-normally distributed variables. Categorical variables were compared using the chi-square or Fisher’s exact tests. Spearman correlation coefficients were calculated to analyze the correlations. All P-values were two-sided, and statistical significance was set at P < 0.05.

Results

Identification of Differential Expression Genes

The study design is presented in Figure 1.

Figure 1 Flow chart of the study design. Abbreviations: FIRDEGs, ferroptosis immune-related differentially expressed genes; GSEA, Gene Set Enrichment Analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI network, Protein-protein interaction network; TF, Transcription factors; Co-DEGs, Common differentially expressed genes; ROC, Receiver operating characteristic curve.

Three datasets, GSE30122, GSE30528, and GSE30529, were used to identify DEGs. Normalized boxplots for each dataset are shown in Supplementary Figure 1. The DEGs were visualized using a volcano map (Figure 2A–C) and heatmap (Figure 2D–F).

Figure 2 Differential gene analysis of diabetic kidney diseases databases. (A) The volcano map of GSE30122. (B) The volcano map of GSE30528. (C) The volcano map of GSE30529. (D) The heatmap of GSE30122. (E) The heatmap of GSE30528. (F) The heatmap of GSE30529. (G) Venn diagram of the Co-DEGs in GSE30122, GSE30528 and GSE30529. (H) Venn diagram of Co-DEGs and FIRGs. Abbreviations: Co-DEGs, Common differentially expressed genes; FIRGs, ferroptosis immune-related genes; DKD, Diabetic kidney disease.

A total of 3521 DEGs were found in GSE30122, including 1320 upregulated and 2201 downregulated genes (Figure 2A). The top ten upregulated DEGs were Ltf, Igj, C3, Serpina3, Evi2a, Col6a3, Cxcl6, Ly96, Cpa3, and Cd53, whereas the top ten downregulated DEGs were B3gnt1, Lpl, Avpi1, Nell1, Hrg, Fgf9, Kat2b, Apold1, Magi2, and Il12ra2 (Figure 2D). GSE30528 contained 3440 DEGs, with 1714 upregulated and 1726 downregulated genes (Figure 2B). The top 10 upregulated DEGs (Olfml3, Serpina3, Igj, C1qb, C1qa, Ltf, Vtcn1, Vsig4, Pycard, and Ap1m2) and the top 10 downregulated DEGs (Arhgap19, Lox, Plce1, C1orf21, Magi2, Lrrc2, Fbxl7, Dach1, Fam65a, and Dpysl3) are shown in Figure 2E. GSE30529 included 4270 DEGs, comprising 1490 upregulated and 2780 downregulated genes in DKD (Figure 2C). Positive correlations among the top 10 DEGs were found for Ltf, Hopx, Psmb9, Igj, C3, C1s, Evi2a, Col6a3, Serpina3, and Tac1, whereas the top 10 downregulated DEGs were Avpi1, Egf, Nell1, Nr0b2, Arsf, Lpl, Thy1, Hrg, B3gnt1, and G6pc (Figure 2F).

Identification of FIRDEGs

When the word “ferroptosis” was used as the search tag, 698 FRGs were retrieved from GeneCards; 64, from MSigDB; and 339, from FerrDb. Together, 853 FRGs were obtained after the union (Supplementary Table 1). When the word “immune” was used as the search tag, 17,731 and 772 IRGs were retrieved from GeneCards and MSigDB, respectively. After intersection of sets, the final number was: 711 IRGs (Supplementary Table 2). After combining FRGs and IRGs, 1500 FIRGs were obtained (Supplementary Table 3). A Venn diagram was created according to the different gene expression patterns of the three datasets (GSE30122, GSE30528, and GSE30529) and 558 co-DEGs were identified (Figure 2G). Finally, the overlap of the co-DEGs and FIRGs yielded 80 FIRDEGs (Figure 2H). The FIRDEGs are shown in Supplementary Table 4.

GO and KEGG Pathway Enrichment Analysis

GO and KEGG enrichment analyses were performed to evaluate the functions of FIRDEGs in patients with DKD. GO analysis consisted of three parts: cellular components (CC), biological processes (BP), and molecular functions (MF). Among them, the JAK−STAT (GO:0007259), MAPK (GO:0051403), cell growth (GO:0016049), and gliogenesis (GO:0042063) pathways were significantly associated with CC. Specific granules (GO:0042581), membrane rafts (GO:0045121), tertiary granules (GO:0070820), and granular membranes (GO:0035579) were significantly associated with BP. For MF, glycosaminoglycan binding (GO:0005539), integrin binding (GO:0005178), cytokine receptor activity (GO:0004896), and amyloid beta binding (GO:0001540) were significant (Figure 3A, C, and E, Table 2). In addition, KEGG pathway analysis revealed some Hippo signaling pathways, comprising ECM-receptor interaction, focal adhesion, lipid and atherosclerosis, chemokine signaling pathway, and PI3K-Akt signaling pathway (Figure 3B, D, and F, Table 3).

Table 2 GO Enrichment Analysis Results of Ferroptosis Immune Related DEGs

Table 3 KEGG Enrichment Analysis Results of Ferroptosis Immune Related DEGs

Figure 3 Go and KEGG Functional enrichment analysis. The bubble chart of (A) GO and (B) KEGG enrichment analyses of FIRDEGs. The circular trellis of (C) GO and (D) KEGG enrichment analyses of FIRDEGs. The bar graph of (E) GO and (F) KEGG enrichment analyses of FIRDEGs. The screening criteria for GO and KEGG enrichment items were P < 0.05 and FDR < 0.05. Abbreviations: FIRDEGs, ferroptosis immune-related differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto encyclopedia of genes and genomes.

Gene Set Enrichment Analysis

We conducted GSEA to identify the most-enriched gene sets in the DKD and normal control groups. GSEA (Supplementary Table 5 and Figure 4A) indicated that the top four most-significantly enriched gene sets related to the DKD group were PI3K pathway (Figure 4B, NES = −1.848, P < 0.01, FDR = 0.041), IL-10 pathway (Figure 4C, normalized enrichment score [NES] = 2.032, P < 0.01, FDR = 0.031), Notch pathway (Figure 4D, NES = −1.456, P < 0.05, FDR = 0.212), and IL-17 pathway (Figure 4E, NES = 1.901, P < 0.01, FDR = 0.031).

Figure 4 Gene set enrichment analysis. (A) GSEA plot showing most enriched gene sets in the DKD group and normal controls. The differential expression genes in GSE30122 were significantly enriched in (B) PI3K pathway, (C) IL10 pathway, (D) Notch pathway and (E) IL17 pathway. The screening criteria for GSEA items were P < 0.05 and FDR < 0.25. Abbreviations: DKD, Diabetic kidney disease; GSEA, Gene set enrichment analysis; NES, normalized enrichment score.

Construction of PPI Network of the FIRDEGs

A PPI network was created for FIRDEGs using the STRING database based on a minimum correlation coefficient greater than 0.400 (Figure 5A). We used five algorithms (MCC, MNC, EPC, degree, and DMNC) to calculate the scores for the node genes. Figure 5B shows the Venn maps of the most common genes of the top 15 FIRDEGs chosen using the five algorithms. Finally, six common hub genes (C-C chemokine ligand 5 [Ccl5], interleukin 18 [Il18], cytochrome b-245 beta chain [Cybb], fc fragment of IgG receptor IIb [Fcgr2b], myeloid differentiation factor 88 [Myd88], and C-C chemokine receptor 2 [Ccr2]) were extracted from the Venn diagrams (Figure 5C).

Figure 5 Protein–Protein Interaction Network. (A) PPI network of FIRDEGs. (B) Venn maps of common genes of the top 15 FIRDEGs were selected using five algorithms: MCC, MNC, EPC, Degree and DMNC. (C) PPI network of the six hub genes. Abbreviations: FIRDEGs, ferroptosis immune-related differentially expressed genes; PPI, protein-protein interaction; Degree, Degree correlation; MNC, Maximum neighborhood component; MCC, Maximal clique centrality; EPC, Edge percolated component; DMNC, Density of maximum neighborhood component.

Multi-Factor Regulatory Analysis

The synthesis and function of protein-coding genes are heavily regulated by TF and miRNAs. We obtained mRNA-miRNA data from the ENCORI and miRDB databases to predict miRNAs interacting with the six hub genes. The intersection of the two databases was visualized using Cytoscape. The mRNA-miRNA interaction network was composed of two hub genes (Il18 and Myd88) (Figure 6A) and 20 miRNA molecules, which constituted 20 pairs of mRNA-miRNA interactions (Supplementary Table 6). The intersection of the CHIPBase and hFtarget databases is shown in Figure 6B. A regulatory network of mRNA-TFs was established involving four hub genes (Myd88, Ccl5, Fcgr2b, and Il18) and 12 TFs, which constituted 17 pairs of mRNA-TF interactions (Supplementary Table 7). The DGIdb database was searched to find potential medications or chemical compounds targeting the hub genes. The mRNA-drug interaction network (Figure 6C) revealed 29 potential medications or chemical compounds linked to five hub genes (Il18, Cybb, Fcgr2b, Myd88, and Ccr2) (Supplementary Table 8).

Figure 6 Multi-factor regulatory analysis. (A) The mRNA-miRNA interaction network of hub genes. (B) The mRNA-TF interaction network of the hub genes. (C) The mRNA–drug interaction network of hub genes. Abbreviation: TF, Transcription factor.

Expression and Correlation Analysis of Hub Genes

To further investigate the expression of six FIRDEGs (Ccl5, Il18, Cybb, Fcgr2b, Myd88, and Ccr2) in DKD, a correlation analysis was performed on their expression levels in GSE30122 (Figure 7A), GSE30528 (Figure 7C), and GSE30529 (Figure 7E) between the normal and DKD groups. Six hub genes (Ccl5, Cybb, Fcgr2b, Myd88, Ccr2, and Il18) in GSE30122, five (Cybb, Fcgr2b, Myd88, Ccl5, and Ccr2) in GSE30528, and four (Ccl5, Myd88, Ccr2, and Fcgr2b) in GSE30529 were statistically significant (P < 0.01). Correlation heat maps among the hub genes in GSE30122 (Figure 7B), GSE30528 (Figure 7D), and GSE30529 (Figure 7F) were plotted to determine the correlations between the six FIRDEGs. The results showed a clear association between FIRDEGs in the three databases.

Figure 7 Expression differences and correlation analysis of hub genes. Group comparison of the expression of hub genes in datasets (A) GSE30122, (C) GSE30528, and (E) GSE30529. Correlation heatmap of hub gene expression in datasets (B) GSE30122, (D) GSE30528, and (F) GSE30529. Vs Normal group, *P < 0.05, **P < 0.01 and ***P < 0.001.

The correlation between FIRDEGs and renal function in clinical patients was analyzed using the NephroSeq database (Figure 8A–L). Ccl5 (R = −0.693, P = 0.026, reporter Woroniecka et al¹⁷) and Il18 (R = −0.711, P = 0.014, reporter Schmid et al³⁴) were negatively correlated with glomerular filtration rate (GFR). Ccr2 (R = 0.675, P = 0.016, reporter Ju et al³⁵), Cybb (R = 0.489, P = 0.047, reporter Ju et al³⁶), Fcgr2b (R = 0.586, P = 0.013, reporter Ju et al³⁶), and Il18 (R = 0.697, P = 0.017, reporter Schmid et al³⁴) were positively correlated with serum creatinine. We constructed ROC curves and calculated the AUC to explore the significance of hub genes in identifying the normal and DKD groups (Figure 9). All six hub genes showed good diagnostic values in the GSE30122, GSE30528, and GSE30529 datasets. Fcgr2b (AUC = 0.957) and Myd88 (AUC = 0.923) in GSE30528 and Ccl5 (AUC = 0.925), Myd88 (AUC = 0.958), and Ccr2 (AUC = 0.983) in GSE30529 had high diagnostic value, with an AUC > 0.9.

Figure 8 Clinical relevance of the expression of ferroptosis immune-related DEGs. Correlation analysis between GFR and (A) Ccl5, (B) Ccr2, (C) Cybb, (D) Fcgr2b, (E) Il18, (F) Myd88; Correlation analysis between serum creatinine and (G) Ccl5, (H) Ccr2, (I) Cybb, (J) Fcgr2b, (K) Il18, (L) Myd88.

Figure 9 ROC curve analysis. ROC curves for (A) Ccl5, (B) Il18, (C) Cybb, (D) Fcgr2b, (E) Myd88, and (F) Ccr2 in the GSE30122 dataset. ROC curve of (G) Ccl5, (H) Il18, (I) Cybb, (J) Fcgr2b, (K) Myd88, and (L) Ccr2 in the GSE30528 dataset. ROC curves of (M) Ccl5, (N) Il18, (O) Cybb, (P) Fcgr2b, (Q) Myd88, and (R) Ccr2 in the GSE30529 dataset. Abbreviations: ROC, Receiver operating characteristic curve; AUC, area under the curve.

Evaluation of Immune Cell Infiltration

The GSE30122 expression profile data were uploaded to the CIBERSORTx website, where the algorithm was performed to determine the relationship between the 22 different types of immune cells and the expression profile data of the samples from different groups (normal/DKD). The composition of the 22 types of immune cells in GSE30122 for each sample is shown as a histogram (Figure 10A); the color corresponds to the proportion of various immune cells in each sample. We compared the abundance of different immune cells between the normal and DKD groups (Figure 10B). The results indicated that M2 macrophages, resting mast cells, and gamma-delta T cells were the most significant infiltrating immune cells (P <0.001). In addition, M1 macrophages (P <0.01), activated mast cells (P <0.01), plasma cells (P <0.01), M0 macrophages (P <0.05), resting NK cells (P <0.05), T follicular helper cells (P <0.05), and regulatory T cells (P <0.05) differed in their infiltration abundance. Correlation analysis was conducted to determine the relationship between FIRDEGs signature and differential immune cell infiltration abundance. In Figure 10C, FIRDEGs are correlated with the infiltration abundance of various immune cells. Among them, Fcgr2b, Cybb, Ccr2, and Ccl5 had a strong correlation with the infiltration abundance of M2 macrophages and gamma-delta T cells.

Figure 10 Immune cell-infiltration analysis of GSE30122 dataset. (A) Bar chart of immune cell infiltration in the dataset. (B) Comparison of immune cell infiltration between normal and NC groups. (C) Correlation heatmap of FIRDEGs and immune cells in GSE30122 dataset. Vs Normal group; *P < 0.05; **P < 0.01 and ***P < 0.001. Abbreviation: ns, no significant.

Validation of the Hub Genes in Mouse Kidney Samples

At 20-weeks of age, db/db mice showed significantly higher blood glucose, kidney/body weight ratio, 24-hour urinary albumin excretion (UAE), serum creatinine, and serum urea nitrogen levels than db/m mice (Figure 11A–E). Pathological staining of mouse kidney cortical area is shown in Figures 11F and G. Compared to the control group, characteristic pathological changes, such as increased glomerular size, cast formation, and tubular dilatation, were observed in the kidneys of the model mice. The Paller score (for evaluating tubular damage) and mesangial area/glomerular area were significantly higher in db/db mice than in db/m mice. These results show that db/db mice have pathologically significant kidney damage compared with db/m mice. Gpx4 and Acsl4 are key proteins involved in the physiological process of ferroptosis. The protein expression of Gpx4 and Acsl4 was significantly different in the kidneys cortical area of db/db mice compared to that in db/m mice (P < 0.05; Figure 12A–D). To confirm the expression levels of the six FIRDEGs, we performed RT-PCR on the kidney tissue of a mouse model. The mRNA expression levels of Ccl5, Il18, Cybb, Fcgr2b, Myd88, and Ccr2 were significantly higher in the db/db group compared to that in the db/m (P < 0.01; Figure 12E–J).

Figure 11 Basic information of mouse model. (A) Kidney weight/body weight ratio. (B) Blood glucose level. (C) 24 h urine albumin excretion levels. (D) Serum creatinine level. (E) Serum urea nitrogen. (F) H&E staining of the mouse kidneys. (G) PAS staining of mouse kidneys. The results are expressed as the mean ± SD. **P < 0.01, ***P < 0.001, ****P < 0.0001 vs db/m. All figures scale bar = 50 μm.

Figure 12 Validation of ferroptosis immune-related DEGs in mouse model. (A) IHC assay of Gpx4 protein expression in the mouse kidneys. (B) IHC assay of Acsl4 protein expression in the mouse kidneys. The AOD of (C) Gpx4 and (D) Acsl4 between two groups. (E–J) Relative expression of mRNA by RT-PCR (n = 6). The results are expressed as the mean ± SD. **P < 0.01, ***P < 0.001, ****P < 0.0001 vs db/m. All figures scale bar = 50 μm.

Discussion

DKD is a major cause of end-stage renal disease and its prevalence is increasing annually. The pathogenesis of DKD is complicated, and there is no effective treatment.² Abnormal iron metabolism is closely linked to impaired renal function and disease progression in patients with DKD. Iron deposition in the tubular epithelial cells of patients might cause nephrotoxicity and kidney damage.³⁷ Gpx4 is a predictor of DKD progression in type 2 diabetes mellitus.³⁸ However, the specific mechanisms underlying ferroptosis in DKD remain unclear. Accumulating evidence suggests that ferroptosis is associated with immunity. Iron is required for activated T cells. T cell proliferation could be inhibited by iron deficiency, whereas iron overload may lead to an imbalance in CD4+ and CD8+ T cells.³⁹ In macrophages, deletion of Slc7a11 promotes ferritin atrophy caused by iron overload-induced ferroptosis.⁴⁰ Therefore, we speculate that ferroptosis is strongly linked to the immune status of DKD and could trigger inflammatory factors that influence the onset of DKD. However, the mechanisms underlying the abnormal ferroptosis and immune inflammation in DKD remain unclear.

In this study, we conducted an integrated analysis to identify the ferroptosis- and immune-related genes in patients with DKD. We downloaded three datasets from the GEO database to identify DKD-related DEGs. By overlapping DEGs with ferroptosis- and immune-related genes, 80 FIRDEGs were identified. The GO results revealed that FIRDEGs were primarily involved in cytokine production, inflammation, and immunity. KEGG enrichment analysis revealed that FIRDEGs were enriched in the PI3K-Akt, NF-kappa B, NOD-like receptor, Toll-like receptor, and IL17 signaling pathways. GSEA indicated that the IL10, IL17, PI3K, and Notch pathways were the most enriched gene sets in DKD. Therefore, FIRDEGs are closely linked to important pathogenic pathways in DKD. Subsequently, we identified six hub genes (Ccl5, Ccr2, Cybb, Fcgr2b, Il18, and Myd88) from the 80 FIRDEGs and predicted mRNA-miRNA, mRNA-TF, and mRNA-drug pairs. To examine the diagnostic and prognostic value of FIRDEGs, we performed ROC curve analysis in DKD; Il18 and Cybb had a certain accuracy (AUC > 0.7 in three datasets), whereas Ccl5, Fcgr2b, Myd88, and Ccr2 were highly accurate (AUC > 0.8 in three datasets) in predicting DKD. In addition, we used RT-PCR to confirm the expression of the six hub genes in the kidney samples from db/db mouse model, and the results were generally consistent with those of RNA sequencing.

Urinary sediment Ccl5 mRNA could be a prospective prognostic biomarker for DKD.⁴¹ Ccr2 is the most common receptor for the C-C chemokine ligand 2 (Ccl2; formerly called monocyte chemoattractant peptide-1).⁴² Blockade of Ccl2/Ccr2 signaling ameliorates kidney injury in db/db mice.⁴³ Cybb encodes the Nox2 protein, which is involved in the oxidative stress response. Inhibition of Nox2 activation could be beneficial for the prevention or treatment of DKD.⁴⁴ Fcgr2b is associated with acute kidney injury⁴⁵ and IgA nephropathy.⁴⁶ Urinary Il18 levels are independently associated with mortality in patients with DKD;⁴⁷ in addition, the Il18 protein level in the renal tissue is considered an effective biomarker for renal pathology.⁴⁸ Myd88 is a crucial adaptor protein involved in TLR signaling. The polymorphisms in the TLR-Myd88-NF-κB signaling pathway confer genetic susceptibility to DKD.⁴⁹ By searching the Nephroseq database, we found that FIRDEGs were correlated with some renal function indicators in patients with DN, indicating that these six genes have clinical value.

DKD is an inflammatory disease driven by immune cells in addition to the traditional metabolic and hemodynamic risk factors.^50,51 In this study, we analyzed the trend of immune cell infiltration in DKD and examined the correlation between FIRDEGs and immune cells. Spearman correlation analysis revealed that the expression profiles of FIRDEGs were linked to certain immune cells, particularly M2 macrophages and gamma-delta T cells. Macrophages can be divided into M0/M1 and M2 types. M0 macrophages are differentiated into two distinct subtypes: classically activated (M1) and selectively activated (M2) macrophages.⁵² M2 macrophages are mainly activated in DKD, and the balance of M1/M2 macrophages is related to the renal microenvironment, which could affect the progression.⁵³ In patients with DKD, T cell infiltration can lead to immune cell accumulation and an increase in the expression of various inflammatory factors in the kidney, eventually aggravating the progression of DKD.⁵⁴ Gamma delta T cells are unconventional T cells that primarily recognize antigens in the absence of classical restrictions via the major histocompatibility complex (MHC).⁵⁵ Gamma-delta T cells can produce inflammatory cytokines, like Il17a, which are involved in the onset of DKD, ischemia-reperfusion injury, fibrosis, and et al^56,57 To some extent, these genes can be used as target genes to study the relationship between ferroptosis and immunity in DKD. The gene-drug network suggests potential drugs for the treatment of DKD. However, this finding must be validated in future studies.

This study has certain limitations as it only provides preliminary clues about the underlying mechanisms rather than conclusive evidence. First, it relied mainly on database mining without sufficient experimental validation. Second, only mouse models were used for validation and the sample size was small. There are several differences between mice and humans at the genomic and biological levels. Further verification with clinical specimens and in vitro experiments is needed in the future. Third, we did not conduct validation after successfully predicting transcription factors, miRNAs, and therapeutic agents. Finally, the co-expression relationship between FIRDEGs and immune cells requires further verification.

Conclusion

This study is the first to analyze DEGs associated with ferroptosis and immunity in DKD and to predict their correlation with immune infiltration. It offers new avenues for the investigation of possible drug targets for the management of DKD-related ferroptosis and immune disorders.