Expression of Ferroptosis-Related Genes is Correlated with Immune Microenvironment in Diabetic Kidney Disease

Introduction

With the continuous growth of the elderly population, the incidence of chronic diseases such as diabetes is increasing year by year, which is a major challenge facing the public medical system.^1,2 Diabetes is a disease with a longer course. If not treated in time, the development of the disease will lead to microvascular lesions and other complications.³ Among them, diabetic kidney disease (DKD) is a high-risk complication caused by progressive diabetes, which leads to a continuous decline in the patient’s renal function and seriously affects the patient’s health and even endangers the patient’s life.^4,5 DKD develops in about 40% of the patients with diabetic and is the first cause of most end-stage renal disease among the world.⁶ The pathogenesis of DKD is complex and multifactorial with the participation of many pathways and mediators, such as hemodynamic abnormalities, metabolic disorders, hormone synthesis, reactive oxygen species (ROS), protein kinase (PKC), mitogen-activated protein kinase (MAPKs) and so on. Each pathway causes damage via multiple regulators or interacts with other pathways. The exact pathogenic mechanism and molecular incidence of DN are still largely uncertain.

Ferroptosis is considered as cell death dependent on intracellular iron, which is characterized by the increase of ROS to lethal levels.⁷ Although the interaction between ferroptosis and acute kidney injury (AKI) had been continuously explored, studies on ferroptosis and chronic renal disease remain lacking. Due to the sensitivity of renal tubules to ROS, ferroptosis also occurs in the development of DKD by increasing oxidative stress and reducing antioxidant capacity.^8–10 Therefore, ferroptosis may provide a new explanation for the pathogenesis of DKD.

Li et al⁹ demonstrated that up-regulating NFE2-related factor 2 (NRF2) by treating with fenofibrate inhibited diabetes-related ferroptosis, which consequently delayed progression of DKD. Feng et al¹⁰ indicated that ferroptosis might enhance DN and damage renal tubules in diabetic models through HIF-1α/HO-1 pathway. Considering the relationship between ferroptosis and ROS, we speculated that ferroptosis may provide a new explanation for the pathogenesis of DKD. Previously, the renin–angiotensin–aldosterone system and oxidative stress factors were considered major causes of DKD. However, increasing evidence revealed immune and inflammation play an important role in DKD.^11–13 Immune cells participate the process of cell reprogramming. In this process, the cells themselves modify the microenvironment by secreting various biological factors, thereby giving the surrounding cells the ability to decide survival and development. In the early stage of DKD, T cells and macrophages migrate and accumulate in the glomeruli and interstitium, which release pro-inflammatory cytokines and ROS.¹³ However, targeted immunosuppressive therapy cannot effectively inhibit the progress of end-stage renal disease (ESRD) including DKD.^14,15 In addition to immune factors, the ultimate precise pathogenic mechanisms of DKD remain elusive. Significantly, various immune cells types are involved in the maintenance of iron homeostasis, which means immune microenvironment (IME) plays a vital role in controlling iron metabolism.¹⁶ Additionally, ferroptosis was found to work synergistically with immunoregulation in tumor immune microenvironment.¹⁷ The lethal ferroptosis in tumors cells can expose tumor antigens, thereby improving the immunogenicity of the microenvironment and enhancing the effectiveness of immunotherapy.¹⁸ Therefore, we hypothesized ferroptosis may be induced by IME, thereby accelerating the progression of DKD.

Here, we use data mining and data analysis techniques, and explore the correlation between ferroptosis-related genes (FRGs) and IME in DKD to provide a new clue for exploring the underlying molecular mechanisms.

Materials and Methods

Microarray Data Acquisition

The expression profile was downloaded from the GEO; By setting “diabetic kidney disease” OR “diabetic nephropathy” as the search keywords, we obtained the GSE30122 dataset, including 9 DKD glomeruli samples, 10 DKD tubuli samples and corresponding control samples, which was available on the GPL571 platform. GSE30122 dataset was composed of the following subseries: GSE30528, GSE30529 and GSE30566. The GSE30528 was selected as the training set (9 DKD and 13 control samples), and GSE30122 as test set (19 DKD and 50 control samples). In addition, we selected GSE96804, a dataset with more DKD patients, as the validation set (41 DKD and 20 control samples). The workflow is shown in Figure 1.

Figure 1 Flow chart of data preparation, processing and analysis.

Construction of Weighted Gene Co-Expression Network Analysis

All genes in the DKD and control samples were sequenced according to median absolute deviation (MAD) values, and the top 50% of the genes were selected for weighted gene co-expression network analysis (WGCNA). At first, the Pearson’s correlation matrices and average linkage method were both performed for all pairwise genes. Then, a weighted adjacency matrix was constructed using a power function. The specific operation parameters and methods referred to our previously published literature.¹⁹ FRGs (including 220 ferroptosis driver genes, 162 ferroptosis-related suppressors and 123 ferroptosis-related markers genes.) were obtained from the ferroptosis-related database FerrDb (http://www.zhounan.org/ferrdb/index.html).

Gene Enrichment Analysis

Gene ontology (GO) enrichment analysis including cell components, molecular functions, and biological processes was based on David database (https://david. ncifcrf.gov/summary.jsp). The potential pathways can be analyzed though the Kyoto Encyclopedia Gene and Genome (KEGG).

Identification of Molecular Subtypes and IME Evaluation

Consistent clustering provided quantitative evidence for determining the number and members of possible clusters, which was performed using ConsensusClusterPlus.²⁰ Then, we determined the clusters with a specific cluster count (k). The ESTIMATE algorithm was performed to calculate stromal score, immune score and ESTIMATE score. TIMER immune infiltrating analysis was used to estimate the abundance of immune infiltrating cells. The single-sample gene set enrichment analysis (ssGSEA) was conducted to evaluate immune signature in DKD samples.

Differentially Expressed Genes (DEGs) Analysis

According to the level of FRGs in DKD samples, molecular subtypes were divided by consistent clustering. The DEGs were identified by linear models for microarray data (Limma) between the clusters. The threshold values were |log₂ FC| >1.5 and P<0.05.

Construction of the PPI Network and Module Analysis

STRING online database (https://string-db.org/) was searched to analyze and evaluate the interaction correlations among DEGs. After removing the unlinked nodes, the remaining genes were constructed into PPI network by Cytoscape (version 3.5.1). Molecular Complex Detection (MCODE), a plug-in of Cytoscape, was widely used to evaluate core networks and molecules in the PPI network. Under statistical parameters of MCODE (scores ≥3 and node ≥4), the molecules in the network may play an important role in regulating or being regulated in ferroptosis process.

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) was used to evaluate the distribution trend of genes in a predefined gene set in the gene table sorted by phenotypic correlation, so as to judge their contribution to phenotype. In order to further explore the function, GSEA was used to explore the functional signaling pathways with the gene set ranging 5 to 5000, and 1000 resampling.

Results

WGCNA and Key Modules Identification

First, the gene expression matrix of GSE30528 was obtained after data preprocessing. Subsequently, we selected genes with variances greater than 50% median absolute deviation (MAD). We used the Pearson’s correlation coefficient to cluster the samples, and no outlier samples were deleted (Figure 2A). When the soft threshold was equal to 9, the connectivity between genes met the scale-free network distribution to construct the WGCNA network (Figure 2B). To classify genes with similar expression profiles into gene modules, average linkage hierarchical clustering was conducted according to the TOM-based dissimilarity measure with a minimum size (gene group) of 50 for the gene dendrogram. Finally, 18 modules were identified with a distance of less than 0.25 (Figure 2C and D). It is worth noting that the darkorange module (r = 0.75, p = 5.5e-5), orange module (r = −0.71, p = 2.4e-4), sienna3 module (r = −0.67, p = 6.6e-4), turquoise module (r = −0.93, p = 2.4e-10) and white module (r = −0.62, p = 2.1e-3) were significantly correlated with disease status, suggesting that these genes in above modules may play a protective role in DKD (Figure 2E). Furthermore, we also confirmed there was a highly significant correlation between Gene Significance (GS) for weight vs Module Membership (MM) in these modules (orange module: r = 0.56, Figure 3A; turquoise module: r = 0.93, Figure 3B; white module: r = 0.68, Figure 3C; darkorange module: r = 0.77, Figure 3D and sienna3 module: r = 0.37, Figure 3E).

Figure 2 Preprocessing of WGCNA. (A) Hierarchical clustering of module hub genes that summarize the modules yielded in the clustering analysis. (B) Analysis of the scale-free fit index for various soft-thresholding powers, and the mean connectivity for various soft-thresholding powers. (C) The cluster dendrogram of genes in GSE30528. (D) Interaction relationship analysis of co-expressed genes. (E) Heatmap plot of the adjacencies in the hub gene network.

Figure 3 Ferroptosis-related Genes (FRGs) in hub gene modules. (A–E) Gene Significance (GS) for weight vs Module Membership (MM) in orange, turquoise, white, dark orange and sienna3 module. (F) Venn diagram to find the FRGs in 5 hub modules.

Functional Enrichment Analysis of FRGs

We overlapped the FRGs in the above five modules. Seven FRGs (LINC00472, NNMT, OSBPL9, VEGFA, UBC, GPX2 and ALOX5) were closely related to the disease status of DKD (Figure 3F). The results of GO showed that FRGs were mainly enriched in intracellular vesicle-related functions (Figure 4A), and KEGG enriched in metabolic related pathways (Figure 4B).

Figure 4 The results of enrichment analysis for ferroptosis-related genes (FRGs). (A) The GO analysis. (B) The KEGG enrichment analysis.

Identification of Two Molecular Subtypes Based on FRGs

GSE30122 was selected as the test set. Cluster analysis was performed to divide the DKD patients into subgroups based on seven FRGs. Based on the evaluation of the cumulative distribution function (CDF) curve and relative change in area under CDF curve, when the k value increases, the offline area of the CDF curve gradually increases. Before the cluster analysis, the validation set (GSE96804) was used to judge the clinical diagnostic value. Except for ALOX5, the expression levels of other FRGs were significantly different between DKD patients and controls (Supplementary Figure S1A). The ROC curves suggested all FRGs have clinical diagnostic value (Supplementary Figure S1B-H). Therefore, seven FRGs including ALOX5 were still performed for subsequent cluster and functional analysis.

Based on the assessment of the downward trend of CDF value, two molecular subtypes were identified based on the expression of FRGs (Figure 5A and B). In addition, according to the evaluation of the average consistency in the cluster group, the number of clusters with the highest average consistency was k = 2 (Figure 5C). Therefore, the optimal clustering stability was identified when k = 2.

Figure 5 Consensus cluster based on ferroptosis-related genes (FRGs). (A-C) k = 2 was identified the optimal value for consensus clustering.

IME Between the Molecular Subtypes

When DKD patients were divided into two clusters according to FRGs, the heat map for cluster is shown in Figure 6A, the expression level of FRGs between clusters is presented in Figure 6B. Next, we performed immune analyses to explore IME between the two clusters. ESTIMATE algorithm revealed that there were significant differences in ESTIMATE score between these two clusters of DKD patients, with no significant difference in stromal and immune score (P < 0.05, Figure 6C). Then, TIMER algorithm indicated there was a significant difference in the degree of T cell infiltration (P < 0.05, Figure 6D). The ssGSEA algorithm was further performed. Interestingly, immunity was mainly concentrated in thick ascending limb and distal convoluted tubule in adult kidney (Figure 6E).

Figure 6 Immune analyses in the 2 clustered molecular subtypes. (A) The expression profile of DKD samples were divided into 2 subtypes; (B) The heatmap visualizing the expression of ferroptosis-related genes (FRGs) in the 2 subtypes; (C-E) ESTIMATE, TIMER and ssGSEA algorithm was further performed to analyze their cell type signature.

DEGs and Functional Analyses

DEGs between the two clusters were identified, and functional analyses were performed to explore the underlying signaling mechanisms. A total of 679 DEGs including 297 genes were upregulated and 400 genes were downregulated (|log₂ FC| >1.5 and P<0.05, Figure 7A and B). GO enrichment analysis revealed that the DEGs were enriched in T cell proliferation and T cell activation related immune function (Figure 7C), which was consistent with the results of TIMER algorithm.

Figure 7 The differentially expressed genes (DEGs) analysis between the clustered subtypes. (A) Volcano map; (B) The heatmap; (C) The visualizing the biological processes enriched by GO analysis.

Then, PPI network (Figure 8A) was constructed and five submodels (Figure 8B-F) were identified by MCODE plug-in, all of which were closely associated with metabolism and immunity, especially, such as organic acid catabolic process and carboxylic acid catabolic process in MCODE 1, immune system process in MCODE 2, T cell activation in MCODE 3 and lipid metabolic process in MCODE 5 (Figure 9A). To identify the potential functions, GSEA was conducted to identify KEGG pathways enriched in the DKD samples between clusters, which following the minimum gene set was 5, the maximum gene set was 5000, and 1000 resampling. Consistently, GSEA revealed that P53 signaling pathway and starch and sucrose metabolism were activated in cluster 1 (Figures 9B).

Figure 8 PPI network was constructed and five submodels were identified by MCODE plug-in. (A) PPI network of all DEGs; (B-F) Five submodels identified by MCODE.

Figure 9 The results of enrichment analysis for differentially expressed genes (DEGs) between the clustered subtypes. (A) GO analysis in 5 MCODE; (B) GSEA.

Discussion

Ferroptosis is a special type of regulatory cell death, which is different from necrosis and apoptosis.²¹ The iron metabolism is regulated by many signals, such as lysosomes activity, excessive activation of heme oxygenase-1 (HMOX1) and the suppression of nitrogen fixation 1 (NFS1).²² Abnormal iron metabolism may cause iron overload, deposition of lipid peroxidation products, ROS overproduction and weakened antioxidant capacity.²³ The accumulation of ROS induced by oxidative stress activates various signaling pathways, leading to DNA damage and cell apoptosis is considered to be one of the important mechanisms.^24,25 In our current study, we firstly found the correlation between FRGs and IME in DKD, which might provide a new clue for exploring the underlying molecular mechanisms.

Initially, we selected the microarray data (GSE30528) as a training set, including 9 DKD glomeruli samples and 13 glomeruli controls. Two molecular subtypes were classified based on the expression of FRGs in test set (GSE30122). Then, TIMER algorithm indicated differences in the degree of T cell infiltration were in distinct clusters associated with FGR. Therefore, we speculated that T cell infiltration may be a factor that induced ferroptosis. The immune disorders of DKD patients were interfered by a variety of factors, especially end-stage patients receiving regular hemodialysis (HD), those who were susceptible to infections.^26,27 It has also been reported that the decrease of T cells in HD patients may be related to the increased tendency of T cell apoptosis.²⁸ DKD is a common complication of type 1 and type 2 diabetes. Type 1 diabetes is an immune system disease, and their CD4+ T cell ratio is increased.²⁹ T cell proliferation, cytokine production and renal T cell infiltration were also observed in DKD animal model, which means T cell activation may play a major role in the initiation of DKD.^30–33 The equilibrium between proinflammatory features and stabilizing regulatory T cells seems to regulate the faint balance between systemic and local requirements for inflammatory stimuli.^34,35 However, the mechanism of T cell infiltration regulates ferroptosis in DKD has not been clarified.

Therefore, in order to further explore the mechanism affecting the IME, we analyzed the expression profiles of DKD samples with different FRGs clusters. A total of 679 DEGs including 297 genes were upregulated and 400 genes were downregulated in cluster 2, which were enriched in extracellular related biological processes and metabolic related biological processes. Similarly, KEGG enrichment analysis also identified some signaling pathways associated with immune and metabolism. PPI network was constructed and five submodels were identified, all of which were closely associated with metabolism and immunity. Consistently, GSEA revealed that P53 signaling pathway and starch and sucrose metabolism were activated in cluster 1. These phenomena indicated that metabolism and immune may be associated with the contribution of ferroptosis in DKD patients. Indeed, multiple factors are involved in metabolism implicated by ferroptosis, such as amino acid and lipid metabolism.³⁶ The energy metabolism pathway of DKD has undergone a series of reorganization and adjustment, which is reflected by starch and sucrose metabolism in GSEA. Ferroptosis is an iron-dependent form of cell death resulting from lipid-based reactive oxygen species, along with increased metabolism of arachidonic acid compromised by conglutination and infiltration of immune cells.^37,38 Numerous metabolic pathways involving iron, lipids, and amino acids regulated the process of ferroptosis in DKD.^39–42 The importance of ferroptosis in regulating IME had been implicated in cancers.^17,43,44 Tumor cells need more iron than normal cells to promote their rapid proliferation, which make them more vulnerable to iron overload and ROS accumulation. The possible mechanism for the correlation between ferroptosis and IME was that ferroptotic cells release HMGB1 (a key factor in immunogenicity), which could attract immune cells to their location.^45–47 Besides, ferroptotic cells could release a potential signal – amino acid oxidation products – to attract immune cells to participate in immune regulation.⁴⁸ Taken together, metabolism disorder initiated in the ferroptosis, which was directly or indirectly affected by IME, leading to cell damage in DKD. There were several limitations to the present study. Firstly, the sample size was small; thus, there may be some deviation in the final results. Secondly, without their own data to verify, the statistical results may be overestimated.

In summary, two molecular subtypes were identified based on FRGs in DKD via consensus clustering. Immune analysis and functional analyses revealed that ferroptosis may be induced by dysregulation of IME, thereby accelerating the progression of DKD, and metabolism-related pathways may play an important role. Our work could shed a novel light on the mechanism of DKD cell death and provides potential target for therapeutic drug for DKD patients.