To identify potential molecular markers associated with oxidative stress during RIRIs, in the first step, we searched in the GEO database (Gene Expression Omnibus, http://www.ncbi.nlm.nih.gov/geo) using the keywords “Renal and IRI” and selected the mRNA expression profiles from three datasets for this study, namely Data GSE148420, using the GPL14746 platform (Agilent-028282 whole mouse genome microarray), which included 4 IRI groups and 4 Sham groups; Data GSE58438, using the GPL11534 platform (Affymetrix Rat Gene 1.1 ST microarray), which included 9 IRI groups and 5 Sham groups; and Data GSE71647, using the GPL7202 platform (Agilent-014868 whole mouse genome microarray), which included 4 IRI groups and 4 Sham groups, To ensure homology of the analytical data, we perform homology gene conversion via the “Homologous Genes” toolbox in the online tool Bioinformatics (http://www.bioinformatics.com.cn). The above three datasets were normalized using the “Limma” R package. At the same time, we used GSE148420 as the training set, and GSE58438 and GSE71647 as the external validation sets. In the second step, we obtained a total of 92 sets of genes (OS) associated with oxidative stress from the Molecular Signatures Database.

Through the Limma R software package, the classical Bayesian algorithm was used to calculate the differentially-expressed genes (DEGs) in the dataset GSE148420. The absolute value of Log2Fold > 1.0 and adj.P<0.05 were used as significance indicators. The Hiplot (https://hiplot.org) visual biomedical online analysis toolbox was used to construct DEG heat maps and volcano maps (15). The intersection of DEGs and oxidative-stress-related genes (OSs) obtained in the Molecular Signatures Database is defined as oxidative-stress-related differential genes (DEOSGs) by us. To reveal the underlying biological functions and underlying mechanisms of DEOSGs, we used the Hiplot (https://hiplot.org) online analysis tool for gene ontology (GO) and Kyoto Encyclopedia of Gene and Genes (KEGG) enrichment analysis of DEOSGs (15). Hiplot is a visualization application equipped with 240+ pieces of biomedical data, covering basic statistics, multiomics, regression, clustering, dimensionality reduction, meta-analysis, survival analysis, risk modeling and other functions (15). GO enrichment analysis includes biological processes (BPs), cellular components (CCs) and molecular functions (MFs); adj. P<0.05 is the threshold for screening the main enrichment functions and pathways of DEOSGs.

To further explore the interactions among DEOSGs, we constructed a PPI network of DEOSGs or protein interactions using the online biological resource database String (http://www.string-db.org/). We used the confidence score of 0.4 as a cut-off value to visualize the PPI network of DEOSGs using Cytoscape software. We identified key genes in the PPI network of DEOSGs using three evaluation metrics from the CytoNCA plug-in (Degree Centrality (DC), Betweenness Centrality (BC) and Closeness Centrality (CC)). In addition, in order to obtain the target key genes more accurately, we used the “glmnet” package for Lasso regression algorithm to re-identify them. Subsequently, the key genes screened through the PPI protein interaction network and the LASSO regression algorithm were intersected, and the intersection gene was defined as the most valuable oxidative-stress-related Hub gene for RIRIs. Subsequently, we performed a ROC curve analysis of the Hub gene using a Hiplot (https://hiplot.org) visual biomedical online analysis toolbox to verify its accuracy, which was considered diagnostic with AUC>0.7 (16). At the same time, we also performed a statistical analysis on the expression of the Hub gene in IRI group and Sham group, where P<0.05 was considered to be statistically significant.

We enrolled a total of 10 adult male Sprague-Dawley (SD) rats from the Laboratory Animal Center of Chongqing Medical University (SYXK[YU] 2022-0016, Chongqing, China), weighing between 250-280 g, all of which were raised under the same conditions. All rats underwent right kidney removal and were randomly divided into two groups, the IRI group and Sham group. In the IRI group, the left renal pedicle was bluntly separated, and non-invasive vascular clamping was performed for 45min, and the color of their kidneys changed from bright red to purple-black, indicating that successful ischemia. After releasing the vascular clamp, it was observed that the color of the kidneys gradually returned to bright red, indicating a successful reperfusion. In contrast, the Sham group only underwent blunt isolation of the left renal pedicle without ischemia treatment. After 24h, both groups of rats died suddenly under anesthesia, and their left kidneys were collected for Western Blot detection to determine the expression of Hub genes. Furthermore, pathological detection was performed to evaluate kidney damage in the two groups. The animal experiments were approved by the Animal Ethics Committee of the Children’s Hospital of Chongqing Medical University (IACUC Issue No: CHCMU-IACUC20220429002).

To further determine the diagnostic value of the Hub genes, first of all, we performed a ROC curve analysis based on the external validation sets GSE58438 and GSE71647 to verify their accuracy, and the Hub genes with the same AUC>0.7 were considered diagnostic. At the same time, we performed a statistical analysis on the expression of the Hub genes between the IRI group and the Sham group in the external validation set, and P<0.05 was considered statistically significant. Secondly, we used the Western Blot method to detect the expression of the Hub genes between the RIRI group and the Sham group; Both GPX3 and GSTT1 antibodies were purchased from Wuhan ABclonal Biotechnology Co., Ltd. (A12856, A1287). The kidney tissues of each group of rats were added to the corresponding lysis buffer, the tissues were homogenized, centrifuged and the supernatant was collected. Protein concentration was determined according to the BCA method. Corresponding tissue protein was taken for SDS-PAGE electrophoresis, followed by membrane transfer. The primary antibodies (both at 1:2000 dilution) was incubated overnight at 4°C. On the following day, the secondary antibodywas incubated at room temperature for 1 hour then washed with TBST again. The band of interest was semi-quantitatively analyzed by calculating the ratio of the band of the protein of interest to the band of the β-actin protein. P<0.05 was considered significant.

The left kidney of the rat was fixed in paraformaldehyde solution, then HE staining and PAS staining were performed after routine dehydration, transparency, embedding and sectioning. Slides were observed under light microscopy. The tubular injury score was rated based on the Paller method (17). The evaluation method was to observe the renal cortex under 400x microscopy, randomly select 10 renal tubules in each high-magnification field, and randomly select 10 fields for a total of 100 renal tubules for each specimen. The tubular injury score was determined according to the following criteria: 1 point for vacuolar degeneration of tubular epithelial cells; Renal tubular conspiration with cell flattening scored 1 point; 1 point for brush edge damage and 2 points for shedding; Epithelial cell necrosis scored 2 points; 2 points for tube formation; Normal tubules score 0. As long as one of the above pathological changes occurred in the selected renal tubules, the corresponding score was recorded. If two or more pathological changes occur in the same renal tubule, the score was scored based on the highest injury score. The final score for each specimen was the total score of 100 tubules.

Kim-1 and NGAL are often used as markers of tubular injury to assess tubular epithelial cell damage. We performed immunohistochemical staining using the EnVision two-step method, strictly following the instruction. The antibodies were diluted 1:200. Both Kim-1 and NGAL were positively expressed in the cytoplasm of tubular epithelial cells. Image J software was used for semi-quantitative protein analysis of Kim-1 and NGAL. The specific operation is as follows: Under a 200x optical microscope, each specimen randomly selected 5 non-overlapping fields of view in the renal cortex and measured their average optical density (OD). The average OD value of the 5 fields of view was the OD value of the specimen, all under the same optical conditions.

Hypoxia is known contributing factor in the development of acute kidney injury (AKI), especially in the case of renal ischemia-reperfusion injury (RIRI). The expression of hypoxia-induced factor 1-α (HIF1-α) can clear ROS and regulate cell survival, exerting a protective effect on RIRI. Heme oxidase-1 (HO-1) is an essential antioxidant in the body and is often used to measure the reduction of ROS-induced damage during ischemia-reperfusion. We performed immunohistochemical staining using the EnVision two-step method, and the primary antibodies were diluted at 1:200, and both HIF1-α and HO-1 were positively expressed in the cytoplasm of renal tubular epithelial cells. Image J software was used for semi-quantitative analysis of HIF1-α and HO-1. The results were assessed as described in the previous paragraph.

We performed a SGEA using the “Clusterproiler” (V3.16.1) R package to elucidate significant functional and pathway differences in the Hub genes between the IRI and Sham group. The number of gene set permutations per analysis is 1000 times. The criterion for a significant gene enrichment was P <0.05.

The proportion of 22 different types of immune cells in the sample GSE148420 of the training set was calculated through Cibersort (http://cibersort.stanford.edu/) (18). The expression level of the 22 types of immune cells in the IRI and Sham group was compared with the “VioPlot” software package. A Spearman-related analysis of infiltrating immune cells with Hub genes was calculated using “Corrplot” in R. We used lollipops to visually diagnose the correlations between the Hub genes and immune cells. P<0.05 was considered statistically significant.

Trrust (https://www.grnpedia.org/trrust/) is a reliable and intuitive tool for transcriptional regulatory networks (19), which can provide key TFs for multiple genes and regulatory information on interactions. We used the TRRUST database to predict TFs that regulated the Hub genes. The Starbase V2.0 database (http://starbase.sysu.edu.cn/) contains data from five databases: TargetScan, Pictar, Pita, Miranda/MirsVR and RNA22, in which the most comprehensive range of miRNA-lncRNAa, miRNA-pseudogene, miRNA-circRNA, miRNA-mRNA and protein-RNA interaction networks is built (20). We used the Starbase database to predict miRNAs that were targeted at and bound to the Hub genes. Finally, we mapped the regulatory network of TF-Hub genes and miRNA-Hub genes through Cytoscape software.

In the Drug-Gene Interactions Database (DGIdb, www.dgidb.org), information on drug-gene interactions and possible druggability of genes from 30 different sources are integrated, synthesized and normalized (21). We used the DGIdb database to predict drugs and molecular compounds that could interact with the Hub genes. Cytoscape software was used to visualize drug-gene interaction networks.

All bioinformatics analyses were performed using R software; Prism software (GraphPad software, La Jolla, CA) was used for statistical analyses of experimental data. **** represents P<0.0001, *** represents P<0.001, ** represents P<0.01, * represents P<0.05.

A flowchart of this study is shown in Figure 1 . In the training set GSE148420, a total of 5456 DEGs (Log2Fold absolute value > 1.0 and adj.P <0.05) were measured, of which 2486 were upregulated and 2970 were downregulated in the IRI group. The volcano map shows the distribution of those DEGs ( Figure 2A ). In order to further identify DEGs related to oxidative stress among differentially-expressed genes, we obtained a total of 92 genes (OS) related to oxidative stress from the Molecular Signatures Database, and obtained 30 genes after an intersection with DEG-related intersection genes, which we defined as DEOSGs ( Figure 2B ). The heat map shows that DEOSGs can be used to distinguish between the IRI and Sham group ( Figure 2C ).

To further determine the potential functions and pathways of DEOSGs, through a GO enrichment analysis, we found that in the BP, DEOSGs were mainly enriched in glutathione metabolism, responses to oxidative stress stimulationthe as well as regulation of T cell activation and apoptosis signaling pathways ( Figure 2D ). Among the CCs, DEOSGs are mainly enriched in peroxisomes, intercellular bridges and NADPH oxidase complexes ( Figure 2D ). In MFs, DEOSGs are mainly enriched in the activation of antioxidant enzymes, glutathione transferase and peroxidase ( Figure 2D ). After a KEGG pathway analysis, we found that DEOSGs were mainly enriched in glutathione metabolism, cytochrome-P450-mediated metabolism of harmful substances, peroxidase, ferrozoticosis and T cell receptor signaling pathways ( Figure 2E ).

First of all, we input the 30 DEOSGs obtained above into the String database, screened the proteins that interacted with them ( Figure 3A ), and imported the results into Cytoscape software to build a PPI network ( Figure 3B ). Using the three evaluation indexes in the CytoNCA plug-in (Degree Centrality (DC), Betweenness Centrality (BC) and Closeness Centrality (CC)), the first 7 genes were identified as key genes, namely: GPX2, GPX3, GPX7, GPX8, GCLC, GSTT1 and GSTM1 ( Figure 3B ).

Secondly, in order to further improve the accuracy of key genes, we used Lasso regression analysis in machine learning algorithms to obtain 5 key genes when the Lambda value was set as 0.0676, namely GLS, GSTT3, NCF4, GSTT1 and GPX3 ( Figures 3C, D ). After intersecting the key genes of PPI with Lasso regression screening, the intersection genes: GSTT1 and GPX3, were finally determined as the most diagnostically-valuable Hub genes ( Figure 3E ).

Finally, in order to verify the accuracy of the Hub genes, we diagnosed the sensitivity and specificity of the predictive model through subject-based operating characteristic (ROC) curve analysis in the training set GSE148420, quantified with the area under the ROC curve (AUC), and the AUC of the Hub genes was greater than 0.7 ( Figure 3F ). The results of statistical analysis between the IRI and Sham group showed that P was less than 0.05 ( Figure 3G ). Therefore, we believe that the Hub genes GSTT1 and GPX3 may have some influence on the pathophysiological process of RIRIs.

Through the above analysis process, we identified GSTT1 and GPX3 as possible biomarkers associated with oxidative stress during RIRIs. In order to further verify their accuracy, we first performed a ROC curve analysis on the expression of Hub genes in external validation sets GSE58438 and GSE71647, and found that AUC was greater than 0.7 (The AUC values for GPX3 and GSTT1 in the GSE71647 and GSE58438 datasets were both 1) ( Figure 4A ). At the same time, the expression of Hub genes in the IRI and Sham group was statistically analyzed, it was found that the expression trend of Hub genes in the validation sets GSE58438 and GSE71647 was the same as that in the training set GSE148420, and there was a statistical difference between their expression (P<0.05) ( Figure 4B ). Secondly, in the rat RIRI model, after a Western Blot detection of kidney tissues, we found that the expression of GSTT1 and GPX3 in the Sham group was significantly higher than that in the IRI group, and the difference was statistically significant ( Figures 4C–F ), which was consistent with our previous results. GSTT1 and GPX3 were further identified as biomarkers associated with oxidative stress in RIRIs.

After evaluating the pathological damage caused by renal ischemia-reperfusion by HE staining, PAS staining and immunohistochemistry, we obeserved that the glomerular and renal tubular structures in the Sham group were normal under the light microscope, as shown in Figures 5A, B . In contrast, in the IRI group, the tubular structure was disordered, the renal tubule was significantly dilated, the tubular epithelial cells were necrotic and detached, and there were debris in the lumen. The tubular injury scores were significantly higher than in the IRI group compared to the Sham group ( Figure 5C ). Furthermore, the expression of renal tubular injury markers Kim-1 and NGAL was significantly increased in the IRI group compared to the Sham group ( Figures 5D, E ). These results are consistent with previous studies, indicating that renal histopathological damage is significantly aggravated when the kidneys are subjected to ischemia-reperfusion.

HIF1-α and HO-1 are often used as markers to assess oxidative stress levels during ischemia-reperfusion. Immunohistochemical detection of markers of oxidative stress, HIF1-α and HO-1 ( Figure 5F ), revealed that their expression was significantly increased in the IRI group compared to the Sham group, with a statistically significant difference ( Figures 5G, H ). It has been demonstrated that during renal ischemia-reperfusion injury, the oxidative stress response is activated and the secretion of antioxidants is increased to counter oxidative stress damage caused by renal ischemia-reperfusion.

We further explored the functions of the Hub genes GPX3 and GSTT1 in RIRIs using GSEA enrichment analysis. The caspase signaling pathway, FAS signaling pathway, IL-7 signaling pathway, NF-κB signaling pathway, B cell survival signaling pathway and TCR signaling pathway were negatively correlated with GPX3 gene expression ( Figure 6A ). TNFR1 signaling pathway, FAS signaling pathway, IL-7 signaling pathway, oxidative stress signaling pathway, IL2RB signaling pathway and 41BB signaling pathway were negatively correlated with GSTT1 gene expression ( Figure 6B ). Through further analyses of these pathways, we found that both IL-7 and NF-κB signaling pathway were key pathways for inducing oxidative stress in tissues. Both the caspase and FAS pathway were key pathways for apoptosis in tissue cells; the 41BB, TCR and IL2RB pathway were closely related to T cell activation and T-cell-mediated immune responses. The B cell survival signaling pathway was related to B-cell-mediated immune responses. Therefore, we believe that there is a strong relationship between the Hub genes GPX3 and GSTT1 and oxidative stress, apoptosis as well as immune responses during RIRIs.

In GO and KEGG enrichment analyses of the above DEOSGs, immune cells also appeared to play a key role in the RIRI process. At the same time, in the single-gene GSEA enrichment analysis of GPX3 and GSTT1, a close relationship was also found between GPX3 and GSTT1 and immune-related signaling pathways. To further confirm the role of immune cells, after an immunoinfiltration analysis, we found significant differences in 7 immune cells of the IRI and Sham group (P<0.05), namely: dendritic cells, macrophages M0, M1 and M2, monocytes, plasma cells, T cell follicular cells ( Figures 7A, B ). To further understand the role of Hub genes in immune infiltration, we used Spearman correlation analysis to determine whether they are associated with immune cell infiltration. A correlation analysis showed that genes GPX3 and GSTT1 were significantly positively correlated with plasma cells and macrophage M0, while negatively correlated with monocytes and macrophages M1 and M2 ( Figures 7C, D ).

We wanted to better understand the role of Hub genes in the regulation of RIRIs. We predicted a total of 81 transcription factors (TFs) capable of interacting with 2 Hub genes using the Strust Database and constructed a TF-Hub gene regulatory network ( Figure 8A ); among them, there were 68 transcription factors that could regulate GPX3 and 13 that could regulate GSTT1. After an intersection of the two, we found that 5 TFs could regulate GPX3 and GSTT1 at the same time ( Figuress 8A, B ). We also predicted miRNAs capable of interacting with 2 Hub genes through the Starbase Database, constructing a miRNA-Hub-gene interaction network with a total of 29 miRNAs interacting with GPX3 ( Figure 8C ) and 20 interacting with GSTT1 ( Figure 8D ). The 2 Hub genes we screened had multiple binding sites for TFs and miRNAs, which could guide our further studies on their mechanism.

Based on the DGIDB database, drugs or molecular compounds that might interact with the Hub genes were predicted, and a total of 13 drugs or molecular compounds that might have regulatory relationships with the Hub genes were screened, among which the number of drugs interacting with GSTT1 was the largest, and a total of 10 related drugs were identified. There were 3 drugs that interacted with GPX3 ( Figures 9A, B ). The 3D structure of the drugs is shown in Figures 9C, D , providing a theoretical basis for further research on Hub genes.