The raw datasets, which include gene expression data for AKI and Control, were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Detailed information was presented in Table 1. In particular, samples in the GSE192883 dataset were classified into several groups based on ischemia time. For this study, we classified the ischemia reperfusion injury (IRI) 28 min and IRI 30min groups as group AKI. Only cortex samples were used in the GSE217427 dataset.

The linear model for high-throughput data analysis (limma) (Ritchie et al., 2015) in Bioconductor was applied to find DEGs by comparing expression value between AKI samples and Control samples in GSE192883 and GSE153625. Differential expression was calculated using an empirical Bayes model. The criteria for the statistically significant difference of DEGs was | log2 fold change (FC)| ≥ 1 in expression and adjusted p-value (false discovery rate, FDR) < 0.05. Volcano plot of all DEGs was performed by ggplot2 package (Ginestet, 2011) in R.

WGCNA was adopted to explore the correlation between genes (Zhou et al., 2022) and identify important module genes in AKI. Firstly, the median absolute deviation (MAD) of each gene was determined, and top 5000 genes with the biggest MAD were included for next step. Secondly, the DEG expression matrix was filtered by the goodSamplesGenes function to omit unqualified genes and samples, and a scale-free co-expression network was built. Thirdly, adjacency was computed using the co-expression similarity-derived “soft” thresholding power (β). The adjacency was then converted into a topological overlap matrix (TOM), and the gene ratio and dissimilarity were determined. The fourth step was the detection of modules using hierarchical clustering and a dynamic tree cut function. Genes with identical expression profiles were classified into gene modules using average linkage hierarchical clustering, with a TOM-based dissimilarity metric and a minimum gene group size (n = 50) for the gene dendrogram. Fifthly, the dissimilarity of module eigengenes was computed, a cut line for the module dendrogram was chosen, and several modules were combined for further investigation. The eigengene network was finally visualized.

The CIBERSORT (Newman et al., 2015), a method using the principle of linear support vector regression to deconvolute the expression matrix of 22 immune cell subtypes, was used to explore the discrepancy in immune cell between AKI and Control samples (Fan et al., 2022). Subsequently, we screened out the immune cells that showed significant differences in infiltration between groups.

The “clusterprofiler” R package (Yu et al., 2012) was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses (Yu et al., 2012). In the GO enrichment analysis, the categories include the cellular component (CC), the biological process (BP), and the molecular function (MF) terms, and adjusted p < 0.05 was regarded as statistically significant differences. In the KEGG pathway enrichment analysis, enriched pathways were identified with an adjusted p < 0.05 (J. Chen et al., 2020).

Considering that proteins rarely work alone, it is necessary to study the interactions among proteins. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) is an online biological resource database (https://cn.string-db.org/) that is commonly used to identify the interactions between known and predicted proteins. By searching the STRING database, the PPI network were selected with a score >0.7, and the PPI network was visualized by Cytoscape software. Finally, the hub genes were screened from this PPI network by Molecular Complex Detection (MCODE) (Bader and Hogue, 2003) and Cyto-Hubba apps (Chin et al., 2014; Chen et al., 2020).

FerrDb V2 is the world’s first database (http://www.zhounan.org/ferrdb/current/) that dedicates to ferroptosis regulators and ferroptosis-disease associations (Zhou et al., 2023). A total of 484 regulatory factors including drivers, suppressors and markers were downloaded from FerrDB V2 database. The commonly expressed DEGs (coDEGs) in GSE192883 and GSE153625 were intersected with the genes obtained from FerrDB V2 to obtain deFRGs. Furthermore, these 33 deFRGs were performed enrichment analysis and PPI network constructions.

Machine learning–based algorithms have been widely used in clinical decision making. Of them, the least absolute shrinkage and selection operator (Lasso) is one of the most commonly used algorithms and its clinical efficacy has been demonstrated previously (Huang et al., 2016; Reichling et al., 2020). Therefore, we choose Lasso model to select the gene signatures associated with AKI. The Lasso model is a dimensionality reduction method for evaluating high dimensional data and was fitted using the “cv.glmnet” function in the R package “glmnet” (Friedman et al., 2010). Firstly, the ROC curves and the area under the curve (AUC) were used to evaluate the diagnostic efficacy of 48 upregulated hub genes in two human datasets (GSE217427 and GSE139061). And then, Lasso model was applied to analysis the 25 hub genes with AUC over than 0.7 in either human dataset to obtain the final AKI-related hub genes.

We created a nomogram model to predict AKI using the R package “rms” (Liu et al., 2021). The diagnostic nomogram model of final 4 hub genes Mcm3, Pir, Atf3, and Lcn2 was constructive in GSE139061 human dataset. The expression of each gene has a corresponding point. The “Total Points” reflected the sum of all the above elements. The ROC curve of the diagnostic nomogram model was performed in GSE139061 human dataset.

C57BL/6 mice (20–25 g) were purchased from the Animal Center of Chinese PLA General Hospital. All animal procedures were approved by the Institutional Animal Care and Use Committee at the Chinese PLA General Hospital and Military Medical College. The 12 male mice were randomly assigned to three groups: Sham group (4 mice underwent sham surgery), bIRI-1d group (4 mice underwent bilateral renal ischemia for 30 min and reperfusion for 1 day), bIRI-7d group (4 mice underwent bilateral renal ischemia for 30 min and reperfusion for 7 days). Renal ischemia and reperfusion and renal sham surgery were performed as described previously (Chen et al., 2022), blood and kidney samples were harvested for further processing.

A quarter of the kidney was fixed in 4% formaldehyde, dehydrated, and embedded in paraffin. Tissue sections (4 mm) were stained with periodic acid–Schiff (PAS). Kidney injury was assessed by measuring the levels of serum creatinine (SCr) and blood urea nitrogen (BUN). Blood samples were collected from the vena cava at the indicated time points, and the serum was separated by centrifugation at 3,000 rpm for 15 min at 4°C and then sent to the PLA General Hospital Biochemistry Department to detect SCr and BUN.

Frozen tissue samples were lysed in TRIzol reagent (Invitrogen, Carlsbad, CA, United States), and total RNA was extracted according to the manufacturer’s instructions. The levels of transcripts were determined by qPCR using TransStartTM Top Green qPCR SuperMix (AQ131, Transgen, Beijing, China) on an Applied Biosystems 7500 system PCR cycler (Applied Biosystems, Foster City, CA, United States). The data were normalized to 18S expression and further normalized to the Control group. Primers were obtained from Genomics (BGI Tech, China). All of these primers are listed in Table 2 .

The workflow of the specific analysis is shown in Figure 1. Firstly, we reanalyzed our previously uploaded dataset GSE192883 (Supplementary Table S1). We combined IRI 28min and IRI 30min group into AKI group, and performed differential gene analysis on these 6 AKI samples and 3 Control samples. Figure 2A showed the principal component analysis (PCA) of these 9 samples, indicating a good distinction between AKI and the Control group. We got 2946 DEGs, including 1399 upregulated genes and 1547 downregulated genes (Supplementary Table S2). The volcano plot clearly presented the expression of all genes between each group, and the top 20 genes with lowest p-values were labeled (Figure 2B).

Then we performed WGCNA analysis on GSE192883 dataset to achieve key modules and genes. After determining the weighting coefficients, the disTOM of 5000 genes was obtained (Figure 2C), and 7 modules, each module represented by a different color. The heatmap displayed the relationship between module eigenvalues and AKI, each column showing the correlation coefficient and the corresponding p-value (Figure 2D). Red represented positive correlations, and blue represented negative correlations, and the darker the color, the larger the correlation coefficient. Figure 2E displayed the clusters and correlation of module eigengenes.

As we can see, turquoise and blue modules had the greatest correlation with AKI, and their correlation coefficients were 0.96 and 0.8, respectively. Therefore, we visualized the correlation between genes in these two modules and AKI ( Figures 2F, G), and it could be seen that most genes of these two modules clustered in the upper right corner, indicating that the module attributes of these genes and their correlation with AKI were high. We took the genes of these two modules as the key genes obtained by WGCNA analysis, and the total number of genes was 4217 (Supplementary Table S3). The Venn diagram in Figure 2H illustrated the genes obtained through both WGCNA and DEGs analyses, revealing a total of 2123 overlapping genes.

In order to obtain more reliable and robust results, we downloaded another dataset GSE153625 (Supplementary Table S4). Figure 3A showed the PCA diagram of samples in GSE153625, indicating a good distinction between AKI and the Control group. We got 2100 DEGs (Supplementary Table S5), including 1035 upregulated genes and 1065 downregulated genes. The volcano plot clearly presented the expression of all genes between each group, and the top 20 genes with lowest p-values were labeled (Figure 3B).

Then we performed WGCNA analysis on GSE153625 dataset to achieve key modules and genes. After determining the weighting coefficients, the disTOM of 5000 genes was obtained (Figure 3C), and 5 modules, each module represented by a different color. The heatmap displayed the relationship between module eigenvalues and AKI, each column showing the correlation coefficient and the corresponding p-value (Figure 3D). Figure 3E displayed the clusters and correlation of module eigengenes. As we can see, turquoise and blue modules had the greatest correlation with AKI, and their correlation coefficients were 0.98 and 0.96, respectively. Therefore, we visualized the correlation between genes in these two modules and AKI (Figures 3F, G), and it could be seen that most genes in these two modules clustered in the upper right corner, indicating that the module attributes of these genes and their correlation with AKI were high. We took the genes of these two modules as the key genes obtained by WGCNA analysis, and the total number of genes was 4249 (Supplementary Table S6). Figure 3H showed the veen map of genes obtained by WGCNA analysis and DEGs analysis, and there were 1952 overlapping genes.

Key genes were found in two datasets by limma and WGCNA analysis, as shown in Figure 4A, after Venn analysis, there were 885 genes expressed in all 4 gene clusters, named as common DEGs (coDEGs). Enrichment analysis were performed to reveal the role of the 885 coDEGs in AKI. KEGG pathway analysis (Figure 4B) showed that, these genes were mainly enriched in Complement and coagulation cascades, Cell cycle, and TNF signal pathway. GO analysis revealed that coDEGs in BP were mainly enriched in anion transport, and fatty acid metabolic process (Figure 4C); The CC were mainly enriched in collagen-containing extracellular matrix, and apical part of cell (Figure 4D); The MF were mainly enriched in active transmembrane transporter activity, and oxidoreductase activity (Figure 4E). The PPI network of coDEGs revealed that coDEGs interact with each other, and the most significant three modules were visualized using MCODE plug-in (Figure 4F). Figure 4G showed the gene nodes degree and MCODE score of 20 hub genes of MCODE Cluster 1. These 20 genes Cdc6, Cdk1, Cdt1, Chek1, Chtf18, Clspn, Dtl, Exo1, Gins2, Lig1, Mcm2, Mcm3, Mcm4, Mcm5, Ncaph, Pold1, Pole, Rad51, Smc2, and Wdhd1 in MOCDE 1 were hub genes.

Many reports have shown that immune changes are prominent during the occurrence and progression of AKI (Allison, 2018; do Valle Duraes et al., 2020; Melo Ferreira et al., 2021), so we paied special attention to the immune infiltration in AKI. We performed immune infiltration analysis on GSE153625 dataset, Figure 5A showed the heatmap of different immune cells expressed in each sample, Figure 5B presented the correlation of 22 kinds of immune cell type compositions. Figure 5C displayed the proportion of 22 kinds of immune cell type in each sample. Figure 5D illustrated the comparison of different kinds of immune cells between AKI and Control groups. These results demonstrated that AKI mice had a higher level of gamma delta T cells, CD4 memory resting cells T cells, resting mast cells, and M2 macrophages. The correlation of 22 types of immune cells revealed that CD4 memory resting cells T cells were positively associated with resting mast cells (r = 0.60), and that plasma cells were positively correlated with monocytes (r = 0.91). In general, various kinds of immune cells were differentially infiltrated in AKI, which could serve as the potential regulation point for AKI treatment.

In order to recognize ferroptosis-related genes in AKI, 484 unique FRGs were selected from FerrDb V2 database (Figure 6A), and they were overlapped with the DEGs and turquoise of GSE192883 and GSE153625 datasets. The results showed that there were 33 deFRGs commonly differentially expressed in both datasets (Figure 6B). Fads2, Dpep1, Cbs, and Hcar1 were downregulated in both datasets, Snca were upregulated in GSE192883 while downregulated in GSE153625, and other deFRGs were upregulated in both datasets (Figures 6C,D). The role of these 33 deFRGs was explored by functional enrichment analysis. The top 7 GO items under each classification were shown in bar charts (Figures 6E–G). The results revealed that these deFRGs were involved in lipid droplet, cellular response to oxidative stress, and enzyme inhibitor activity. The KEGG pathway enrichment analysis indicated significant enrichment of deFRGs in the terms of Human T cell leukemia virus 1 infection, AGE−RAGE signaling pathway in diabetic complications, and HIF1 signaling pathway (Figure 6H). Furthermore, PPI analysis demonstrated interactions among deFRGs. Cdkn1a, Tert, Jun and Nras were grouped into MCODE cluster 1 with Jun having the highest MCC score (Figure 6I).

We combined 20 MCODE cluster 1 genes of coDEGs and 33 deFRGs as hub genes for AKI. Among these 53 hub genes, there were 48 genes upregulated in both GSE192883 and GSE153625 datasets (Table 3). We then performed ROC analysis on these 48 upregulated hub genes in two human datasets GSE217427 and GSE139061. The results showed that there were 20 hub genes with AUC over than 0.7 in GSE217427 human dataset (Figure 7A), while there were 8 hub genes with AUC over than 0.7 in GSE139061 human dataset (Figure 7B). Among these 25 genes with AUC over than 0.7, LCN2, PLin2, and ATF3 had an AUC over than 0.7 in both two human datasets. Further machine learning analysis (Lasso analysis) was performed on these 25 hub genes in the GSE217427 (Figures 7C, D) and GSE139061 human datasets (Figures 7E, F). Through Lasso analysis, ATF3, CHEK1, ETV4, LCN2, MCM3, NRAS, PIEZO1, PIR, TERT, and WDHD1 were identified in GSE217427. In GSE139061 dataset, ATF3, CHTF18, EGR1, HILPDA, LCN2, MCM3, NQO1, PIR, TIMP1 were screened. Four hub genes (LCN2, ATF3, PIR, and MCM3) were found in both datasets. A diagnostic nomogram model was constructed using these final four hub genes from the GSE139061 human dataset to predict an individual’s risk of developing AKI. In addition, the ROC curve of the diagnostic nomogram model was depicted in the GSE139061 human dataset to assess its diagnostic ability (Figure 7G). The AUC of this model was 1 (Figure 7H), indicating that this diagnostic nomogram model of LCN2, ATF3, PIR, and MCM3, can effectively distinguish AKI samples. Therefore, we hypothesize that these genes are highly potential biomarkers for AKI.

To validate the above results and expression patterns of these four key diagnostic signature genes, we first reanalyzed our previously uploaded dataset GSE192883 and another AKI dataset GSE98622 to demonstrate the mRNA expression of these four genes at different ischemia times and different reperfusion times. We found that these four genes began to upregulate within 16–18 min of ischemia, indicating their sensitivity to mild ischemia. The expression level of severe ischemia was higher in Ischemia 28 min, indicating a positive correlation with the degree of ischemia (Figure 8A). The results of IRI-AKI with different reperfusion (Figure 8B) showed that the Atf3 gene began to express rapidly at 2 h after ischemia, which could predict the occurrence of AKI early and its expression remained elevated even after 24 h. The expression of Lcn2 gradually increased at 4 h after reperfusion, peaked at 24–72 h, and remained elevated at 7d after reperfusion. The expressions of Pir and Mcm3 began to increase at 24 h and reached a peak at 24–48 h. These results indicated that all four genes were highly responsive to AKI with varying degrees of ischemia and time of reperfusion. Secondly, we conducted our own bilateral IRI (bIRI) model, as shown in Figures 8C–G, the pathological injury was severe 1 day after bIRI and persisted for 7 days. The qPCR results (Figures 8H–K) demonstrated a high level of consistency between RNA-seq data and qPCR results. Finally, we utilized the Kidney Interactive Transcriptomics (KIT) online tools (http://humphreyslab.com/SingleCell/) to identify the expression pattern of these four genes through single cell sequencing and spatial transcriptomic analysis in GSE182939 (Dixon et al., 2022) (Figure 8L). The results (Figure 8M-P) showed that, Atf3 was mainly expressed in proximal tubule segments 3 to descending thin limp of loop of Henle, with the highest expression at 4 h after IRI; Lcn2 was mainly expressed in the descending and ascending thin limp of loop of Henle, as well as principal cells, with the highest expression at 12 h after IRI; Mcm3 was mainly expressed in proximal tubule segments 3 to the descending thin limp of loop of Henle, with the highest expression at 2d after IRI; Pir was mainly expressed in urothelium, with the highest expression at 2d after IRI.