DOI: https://doi.org/10.21203/rs.3.rs-1945651/v1
Diabetic kidney disease (DKD) is a major public health issue because of its refractory nature. Ferroptosis is a newly coined programmed cell death characterized by the accumulation of lipid reactive oxygen species (ROS). However, the prognostic and diagnostic value of ferroptosis-related genes (FRGs) and their biological mechanisms in DKD remain elusive. The gene expression profiles GSE96804, GSE30566, GSE99339 and GSE30528 were obtained and analyzed. We constructed a reliable prognostic model for DKD consisting of eight FRGs (SKIL, RASA1, YTHDC2, SON, MRPL11, HSD17B14, DUSP1 and FOS). The receiver operating characteristic (ROC) curves showed that the ferroptosis-related model had predictive power with an area under the curve (AUC) of 0.818. Gene functional enrichment analysis showed significant differences between the DKD and normal groups, and ferroptosis played an important role in DKD. Consensus clustering analysis showed four different ferroptosis types, and the risk score of type four was significantly higher than that of other groups. Immune infiltration analysis indicated that the expression of macrophages M2 increased significantly, while that of neutrophils decreased significantly in the high-risk group. Our study identified and validated the molecular mechanisms of ferroptosis in DKD. FRGs could serve as credible diagnostic biomarkers and therapeutic targets for DKD.
DKD also called diabetic nephropathy (DN), is the most common cause of end-stage renal disease (ESRD) as a refractory chronic microvascular complication of diabetes mellitus (DM). The main pathological features of DKD are tubular atrophy and tubulointerstitial fibrosis1. Currently, only angiotensin converting enzyme inhibitors (ACEI), angiotensin receptor blockers (ARB), and sodium-glucose cotransporter-2 inhibitors (SGLT2i) have been proven to provide partial renoprotection in the progression of DKD2. As current treatments for DKD are limited, it is difficult to inhibit DKD from progressing to ESRD, which has become a conundrum for nephrologists and endocrinologists3. Therefore, more sensitive biomarkers for early diagnosis, intervention, better classification, and management of DKD are needed.
Ferroptosis is a recently identified atypical form of programmed cell death induced by oxidative stress damage4. To date, few studies have used bioinformatics methods to improve the understanding of genes related to ferroptosis in DKD. Hu et al.5 obtained diabetic nephropathy and normal kidney samples from GSE96804 dataset. Six hub genes related to ferroptosis in DKD were analyzed. However, this study had small samples from a single microarray analysis, which may have resulted in a high false-positive rate. Therefore, it is necessary to identify new prognostic ferroptosis markers for the diagnosis and treatment of DKD using comprehensive bioinformatics analyses.
We integrated four microarray datasets from the Gene Expression Omnibus (GEO) database to identify the significant differentially expressed genes (DEGs). The least absolute shrinkage and selection operator (LASSO) regression was used to construct and verify the diagnostic model. To further explore the biological processes of DKD, we conducted a comprehensive enrichment analysis using gene ontology (GO) database, Kyoto encyclopedia of genes and genomes (KEGG) pathway database, gene set enrichment analysis (GSEA), gene set variation analysis (GSVA), and single sample GSEA (ssGSEA) methods. Weighted gene co-expression network analysis (WGCNA) was used to identify co-expressed gene modules and explore the core genes in the network. As the immune system plays an important role in the progression of DKD, we used CIBERSORT to conduct the immune infiltration analysis. Our study establishes a comprehensive network of FRGs related to DKD, providing useful evidence for identifying the role of ferroptosis in the diagnosis and targeted therapy of DKD.
Differentially Expressed Genes
The DKD dataset GSE96804 from the public data platform was downloaded and homogenized to further study the characteristics of DKD and explore efficient markers for early diagnosis and treatment (Fig. 1a–b). Principal component analysis (PCA) showed that the expression profile of normal kidney tissue was significantly different from that of the DKD group (Fig. 1c). DEGs were analyzed between the two groups using the limma package, including 516 upregulated and 441 downregulated genes. The results were visualized using volcano and heat maps (Fig. 1d–f).
We analyzed the enrichment scores of related data sets in the MSigDb database of all samples using GSVA and screened the differential enrichment pathways using the limma package with the filter criteria of P < 0.01 (see Supplementary Table S1online). We found significant differences in multiple signaling pathways between the two groups, which were visualized using a heat map (Fig. 1e). All of them proved statistically significant differences between the normal and DKD groups in multiple ferroptosis-related pathways. Next, we used GSEA to perform enrichment analysis on the differential gene list ranked from high to low in the logFC of the GO, KEGG, and MSigDb pathway-related datasets (see Supplementary Table S2-3 online). The results showed that the collagen-containing extracellular matrix and extracellular matrix were mainly enriched using GO enrichment analysis, while the mitochondrial inner membrane and mitochondrial matrix (Fig. 2a–b) were inhibited. In the pathway enrichment analysis, focal adhesion and microRNAs in cancer were enriched, while chemical carcinogenesis-reactive oxygen species and non-alcoholic fatty liver disease were inhibited (Fig. 2c–d). Multiple pathways were enriched in the MSigDB-related dataset enrichment analysis, ferroptosis pathway enrichment was inhibited, and the results were consistent with GSVA (Fig. 2e–f). Based on these results, we speculated that ferroptosis might be related to DKD (Fig. 3a). A total of 382 FRGs were downloaded from the FerrDb. The expression of 27 FRGs was significantly different between the normal and DKD groups (Fig. 3b). The FRGs were visualized using a heat map (Fig. 3e). We used the ssGSEA algorithm to calculate the ferroptosis gene enrichment score of each sample based on FRGs. We found that the degree of ferroptosis enrichment in the DKD group was significantly lower than that in the normal group, which was consistent with the GSEA and GSVA results. At the same time, the ROC curve proved that the ferroptosis-related gene set had a significant diagnostic value (Fig. 3c–d). The above results showed a significant difference between the DKD and normal groups and that ferroptosis played a certain role in DKD.
Construction and verification of Prognostic Model of FRGs
We used WGCNA co-expression network analysis to search for gene modules related to ferroptosis and DKD. We obtained 10 consensus modules according to the results of the WGCNA. Next, we analyzed the correlation between the gene module and phenotypic data based on DKD and the enrichment score of ferroptosis. The results showed that the blue module (3,788) was positively correlated with DKD and negatively correlated with ferroptosis, while the green (576), green-yellow (114), and yellow modules (1,237) were positively contrary (Fig. 4a–e). We obtained the key genes of each module with a correlation coefficient of 0.8 as the threshold and intersected with the differential genes. Finally, 305 ferroptosis gene sets related to DKD diagnosis were obtained (Fig. 5a). A diagnostic model was constructed by LASSO regression using the GLMNet package (Fig. 5b–c). Finally, we constructed a diagnostic model for DKD using eight FRGs (SKIL, RASA1, YTHDC2, SON, MRPL11, HSD17B14, DUSP1 and FOS). The results showed that the risk value in the DKD group was significantly higher (Fig. 5d–e). We then analyzed the functions of these eight FRGs. Friend analysis showed that the FOS gene had the highest correlation with the other seven genes (Fig. 6a). GO functional analysis indicated that multiple enzyme activities and transcriptional processes (Fig. 6b) may be involved (see Supplementary Table S4 online). KEGG functional analysis suggested that the MAPK signaling pathway, fluid shear stress, and atherosclerosis pathways were enriched (Fig. 6c–d) (see Supplementary Table S5 online). Immune cell infiltration of all samples was evaluated through 22 immune cell gene sets using the CIBERSORT algorithm (Fig. 6e). We found that in the high-risk group, the expression of macrophages M2 increased significantly, whereas the expression of neutrophils decreased significantly (Fig. 6f).
The GSE30566, GSE99339, and GSE30528 datasets were used to verify the effectiveness of the prognostic model. After the de-batch effect and normalization processing, we obtained the verification data set (normal group 26; DKD group 23) (Fig. 7a). The risk score for each sample in the dataset was calculated based on the correlation coefficient of our prognostic model. The box chart showed that the risk value of the DKD group is significantly higher. The ROC curve was further used to verify the accuracy of the model, and the results showed that the ferroptosis-related prognostic model still had a high diagnostic value (AUC = 0.818) (Fig. 7b–d). All results proved that the prognostic model we constructed was accurate and repeatable.
Analysis of Subtypes of Ferroptosis
We used previously obtained FRGs to construct a molecular subtype model to further study the molecular subtypes of ferroptosis in DKD. According to the expression matrix of FRGs, all cases were analyzed using consistent cluster analysis to determine their potential ferroptosis types using the ConsensusClusterPlus package. The results showed four different ferroptosis types after 1,000 repeated sampling (Fig. 8a–d), and the risk score of type four was significantly higher than that of the other groups. This may provide a theoretical basis for exploring personalized treatment measures for DKD in the future (Fig. 8e).
DM and its complications pose a major public health issue. The updated prevalence of diabetes in adults is 12.4% in China6. DKD, one of the main chronic complications of DM, is currently the main cause of renal replacement therapy7. As the concrete mechanism of DKD is unclear, there are no effective medicines or methods to prevent ESRD. Ferroptosis is an iron-dependent cell death mode induced by the excessive accumulation of lipid peroxidation products8,9. ROS plays a key role in ferroptosis10,11.The discovery of ferroptosis has provided a new understanding of the pathogenesis of various diseases12.
In the present study, we performed differential expression and principal component analysis. The results showed that the expression profiles were significantly different between the DKD and normal groups. Functional enrichment analyses based on DEGs showed that ferroptosis-related pathways were significantly different between the two groups. The ferroptosis enrichment degree in the DKD group was significantly lower than that in the normal group. GO enrichment analysis showed that the collagen-containing extracellular matrix and extracellular matrix were enriched, while the mitochondrial inner membrane and mitochondrial matrix were inhibited. Abnormal accumulation of extracellular matrix produced by endothelial cells and podocytes leads to thickening of the glomerular basement membrane, an early pathological change in DKD13. The prognostic performance of the ferroptosis-related gene set was verified using the ROC curve analysis, which indicated that the model had a good prognostic value. Immune cell infiltration was further analyzed between the high- and low-risk groups. The results showed that macrophage M2 increased significantly in the high-risk group. This may be due to an overactivated immune response in patients with DKD. These results suggest that immunotherapy may be a therapeutic target for DKD.
Eight key FRGs (SKIL, RASA1, YTHDC2, SON, MRPL11, HSD17B14, DUSP1and FOS) were identified, and their validity in predicting the prognosis of DKD was analyzed. SKIL (also known as SnoN) is a regulator of the transforming growth factor-β (TGF-β) signaling pathway, which acts as an antifibrotic factor in the pathological process of DKD14. SKIL gene depletion prompted epithelial- mesenchymal transited into renal tubular cells in the condition of high glucose. Loss of SKIL expression appears to exacerbate progressive renal fibrosis in DKD15. Li et al.16 found that high glucose induced the downregulation of SKIL through the TGF-β1/Smad signaling pathway in human renal tubule epithelial cells. Bone morphogenetic protein7 (BMP-7) ameliorates renal fibrosis by increasing the expression of SKIL in renal tubular epithelial cells17. However, the relationship between ferroptosis and SKIL in DKD remains unknown. RASA1 is a RasGAP signaling scaffold protein involved in various physiological processes, such as cell proliferation, differentiation, and apoptosis18. RASA1 inhibits renal tissue fibrosis by reducing myofibroblasts proliferation 19. In kidney carcinoma, RASA1 reduces miR-223-3p expression to inhibit the proliferation and differentiation of renal cell carcinoma20. YT521-B homology domain containing 2(YTHDC2) is an m6A reader that expedites messenger ribonucleic acid (mRNA) decay. YTHDC2 increases the translation efficiency of target genes and reduces their mRNA abundance21. Ma et al.22 found that YTHDC2 is a powerful inducer of ferroptosis and that increasing YTHDC2 is an alternative therapy for lung adenocarcinoma targeting ferroptosis. SON is a ubiquitously expressed and evolutionarily conserved DNA and RNA binding protein localized in nuclear speckles. SON is involved in multiple cellular processes, including transcription, RNA splicing, and gene repression, which regulate the cell cycle and preserve stem cells23, 24.Based on the above mechanism, SON plays a role in various diseases such as cancer, influenza, and hepatitis25. Overexpression of SON is involved in aberrant transcriptional initiation in leukemia26. MRPL11 contribute to protein synthesis as a mitochondrial ribosomal protein within the mitochondria. It mediates aerobic energy conversion through the oxidative phosphorylation system to affect the pathophysiological processes of various tumors27, 28. The protein-coding variants of the hydroxysteroid 17-β dehydrogenase 14 gene (HSD17B14) can prevent the progression of type 1 DM to ESRD29. Dual-specificity phosphatase 1 (DUSP1), a regulator of the MAPK family, is associated with various pathological changes in the kidney, including renal hypertrophy, renal fibrosis, and glomerular apoptosis. Y et al.30 proved that DUSP1 is involved in renal fibrosis in DKD through the miR-324-3p/DUSP1 axis. Another study demonstrated that DUSP1 could release DKD by targeting the JNK-Mff-mitochondrial fission pathways31. Researchers have found that DUSP1 can inhibit autophagy-dependent ferroptosis in human pancreatic cancer cells32. FOS was identified to play an important role in various kidney diseases such as membranous nephropathy33, immunoglobulin A nephropathy34, 35, and chronic glomerulonephritis36.All FRGs have not been fully illustrated in the development of DKD. Further experiments are required to verify the functions of the key genes in DKD.
KEGG pathway enrichment analysis revealed that the MAPK signaling pathway was significantly enriched. The MAPK signaling pathway plays a crucial role in various physiological processes such as proliferation, differentiation, and metastasis37,38. It is involved in ferroptosis as a regulator of oxidative stress, which regulates signal transduction in a ROS-induced manner 39. Poursaitidis40 showed that the inhibition of MAPK signaling protected lung cancer cells against ferroptosis. Wen-Tsan Chang41 demonstrated that the drug could induce hepatocellular carcinoma cell death through the MAPK pathway in the form of ferroptosis. In acute myeloid leukemia cells, inhibition of the MAPK pathway can render acute myeloid leukemia cells insensitive to ferroptosis42. Nevertheless, how FRGs affect the pathophysiological process of DKD through the MAPK pathway needs to be fully studied.
Our study has some limitations. First, only diabetic glomerular tissue samples were included, which may have led to one-sided results and selection bias. Therefore, in future studies, it will be necessary to improve the detection capability by integrating data from multiple tissue samples. Second, the sample size was relatively small, which may have resulted in a false-positive rate. This will facilitate an increase in sample size for further validation. Third, owing to the lack of pathological specimens of DKD in clinical settings, we were unable to assess the associations between risk indicators and pathological subtypes. In future studies, more pathological subtypes of DKD are needed to conduct further analyses. Fourth, our results are based on bioinformatic analysis, therefore, require further in vitro and in vivo verification. Fifth, our screening method can be used for gene screening in terms of its diagnostic value, phenotypic module clustering, differential expression and co-expression analysis, and clinical predictive models. It also has good value in screening ferroptosis-related molecules with diagnostic value in DKD. Based on these analyses, there were only eight genes left; thus, it is impossible to further screen key genes by protein-protein interaction (PPI). In the follow-up study, we will consider using a PPI interaction network to screen molecules further. Sixth, in our study, the effect of the ferroptosis pathway on the upstream and downstream regulation mechanisms of DKD was not considered. In a follow-up study, we will further study the multi-level regulatory mechanism of ferroptosis on DKD at epigenetic, transcriptional, and post-transcriptional levels through various molecular experiments and bioinformatics methods.
While an effective treatment for DKD has not yet emerged, we integrated comprehensive bioinformatic analyses to identify the biological functions and pathways associated with ferroptosis in the development of DKD. We also identified eight FRGs as potential diagnostic biomarkers and therapeutic targets for DKD. Our studies emphasize the role of ferroptosis in the progression and treatment of DKD.
Data Collection and Preprocessing
We extracted the gene expression profiles of DKD (GSE96804,GSE30566,GSE99339 and GSE30528) from the GEO database using the GEOquery package43. The GSE96804 dataset 44 from Homo sapiens, based on the GPL17586 platform, contains 61 samples, including 20 normal glomerular and 41 diabetic glomerular tissues. All the samples were included in this study. The GSE30566 dataset45 from Homo sapiens based on the GPL571 platform contains 26 samples, including 13 normal glomerular and 13 normal renal tubular control; 13 normal glomerular samples were included in this study. The GSE30528 dataset45 from Homo sapiens, based on the GPL571 platform, contains 26 samples, including 13 normal glomerular and 9 diabetic glomerular tissues. All the samples were included in this study. The GSE99339 dataset46 was obtained from Homo sapiens, and the data platforms were GPL19109 and GPL19184. It contained 184 samples, including 13 diabetic glomerular tissues, all of which were included in this study. The data were normalized and de-batch processed using the sva package47 and standardized using the limma package48.Subsequently, 382 FRGs were obtained from the FerrDb49. Our study is based on open-source data; therefore, there are no ethical issues or conflicts of interest.
Construction and Verification of the LASSO Model
Currently, LASSO regression is a commonly used machine learning algorithm for the construction of diagnostic models. Regularization was used to solve the occurrence of overfitting in the process of curve fitting and to improve the accuracy of the model. The model was built using the GLMnet package50 with a parameter set.seed (2), family = "binomial."
Difference Expression Analysis
We used the limma package to calculate the differential expression of genes between the normal and DKD groups in the GEO microarray data with log fold change (logFC) > 0.5 and adjusted P-value < 0.01 as the threshold. The genes were upregulated in the DKD group if logFC > 0.5, and downregulated if logFC < 0.5. The results of the differential expression analysis are shown in the heat map using the R package pheatmap and the volcano map using the ggplot2 package51.
Functional Enrichment Analysis
GO analysis was used to conduct large-scale functional enrichment, including biological process (BP), molecular function (MF), and cellular component (CC) analyses. The KEGG database stores the genomes, biological processes, diseases, and medical information. GO biological process and KEGG pathway enrichment analyses of the DEGs were performed using the R clusterProfiler package52. The critical value of FDR< 0.05 was considered statistically significant52.
To study the differences in biological processes among different groups, we used GSEA for enrichment analysis according to the logFC arrangement based on the GSE96804 profile. GSEA is a computational method used to analyze whether a particular gene set has a statistical difference between the two biological states. In our study, we used GSEA to explore the differences in pathways and biological processes of the samples in the datasets. The "msigdb.v7.0.symbols" gene set were downloaded from the MSigDB53database for GSEA analysis.
In addition, the enrichment scores of related pathways in the MSigDB database were calculated according to the gene expression matrix of each sample using the GSVA method using the R-packet GSVA54.The differences in samples were screened using the limma package48. Enrichment items with statistically significant differences are shown in the heat map. According to the gene expression matrix of each sample, the enrichment scores of the FRGs were calculated using the ssGSEA method and are displayed in boxplots.
Gene Co-expression Analysis
WGCNA55 aims to identify co-expressed gene modules, analyze core genes in the network, and explore the relationships between modules and phenotypes. First, the soft threshold was calculated using the pickSoftTreshold function, and five was the best soft threshold. We then built a scale-free network based on the soft threshold. Next, a topology matrix and hierarchical clustering were constructed. We dynamically cut and identified the gene module and calculated the eigengenes, and the number of genes in each module was at least 50. The correlation between modules was constructed according to the eigengenes of the modules, and hierarchical clustering was performed. The modules were merged again with a correlation of more than 0.7, and finally, 10 modules were obtained. The correlations between modules and clinical features were explored using Pearson’s correlation analysis.
Consistent Cluster Analysis
Consistent clustering is a method that can be used to determine the members and number of possible clusters in a dataset (microarray gene expression). In order to distinguish different subtypes of DKD, we carried out a consensus clustering analysis related to FRGs in the DKD group in the GSE96804 dataset using "ConsensusClusterPlus" R package56. In this process, the number of clusters was set between 2 and 10; 80% of the samples were taken each time and calculated 100 times, clusterAlg = "hc” and distance = "euclidean."
Immune Infiltration (CIBERSORT)
CIBERSORT57 is an algorithm for deconvolution of the transcriptome expression matrix, according to the principle of linear support vector regression, to estimate the composition and abundance of immune cells among mixed cells. The gene expression matrix data (TPM) were uploaded to CIBERSORT, combined with the LM22 gene matrix, and an immune cell infiltration matrix with filtration (P < 0.05) was obtained. A bar chart was drawn using the ggplot2 package in R language to show the distribution of 22 immune cell infiltrations in each sample.
Statistical Analyses
All data calculations and statistical analysis were conducted using R programming (https://www.r-projec t.org/, 4.0.2 version). For two groups of continuous variables, the statistical difference in normal distribution variables was estimated using an independent Student’s t-test, and non-normally distributed variables were analyzed using the Mann–Whitney U test (Wilcoxon rank sum test). A two-tailed P < 0.05 was considered statistically significant.
Acknowledgements
This study was supported by a grant from the Science and Technology Foundation of Nanjing Medical University (NMUB2020136).
Author contributions statement
D.L. wrote the manuscript. All authors contributed equally to researching data for the article. All authors have read and agreed to the published version of the manuscript.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Additional Information
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.