All animal experiments were performed in accordance with international guidelines and approved by the Institutional Animal Care and Use Committee of the PLA General Hospital. Male C57BL/6 mice aged 6–8 weeks were purchased from GemPharmatech Co., Ltd. Mice were housed in a pathogen-free, constant temperature environment with a 12 h light/dark cycle and allowed to acclimatize for a week in the animal facility before the operation.

Mice were subjected to bilateral renal I/R injury as previously described (Zhou et al., 2012). Briefly, mice were anesthetized with isoflurane on a thermostatic blanket while maintaining the temperature of the mice at approximately 37°C. Bilateral renal pedicles were dissociated through a middle abdominal incision, and the bilateral renal pedicles were clamped with atraumatic vascular clips for 18 or 30 min (n = 3). The incision was closed in two layers, and the mice were injected subcutaneously with 0.3 ml of warm saline on the back after surgery for volume supplementation. Sham-operated mice underwent the same procedure but without clamping of the renal pedicles.

The edgeR package from Bioconductor was used to explore DEGs. DEGs were analyzed using the edgeR package (Robinson et al., 2010). DEGs were finally determined using parameters of false discovery rate adjusted p value <0.01, and log2 fold change >1.2 or < −1.2 unless specified. The Mfuzz (Kumar and M) package from Bioconductor was used for the discrete time points analysis. GO and KEGG analyses were performed to enrich the biological processes and molecular pathways of different gene sets (Ashburner et al., 2000; Yamada et al., 2006). GSEA was used to identify the activation or inhibition of the screened pathways (Subramanian et al., 2005). Enrichment p-values were adjusted using the Benjamini and Hochberg method, and p-adjusted values of less than 0.05 were determined to be significantly enriched.

DNB analysis was performed as previously reported (Ideker and Sharan, 2008; Jin et al., 2008; Liu et al., 2014). DNB is a dynamic analysis conducted through the variation of genes themselves and the correlation between genes, A gene group with larger variation and less correlation with other genes means that the gene group has outlier performance and becomes unstable, thus determining that the biological process is transforming to another stage. When another stage of the disease is reached, this phenomenon disappears. Briefly, we selected genes with significantly high deviations at each time point. Next, we built a distance matrix with Pearson correlation coefficients (PCC) for each pair of genes and set a threshold of 0.9 to group genes. The composite index (CI) was calculated from the average standard deviation (SD) of DNB molecules, the average absolute PCC value inside the DNB cluster (PCC₁), and the average absolute PCC value of the DNB cluster and others (PCC₀), as follows: C I = P C C i P C C 0 S D

After 24 h of reperfusion, the kidneys of all mice were used for RNA extraction. TRIzol Reagent (Invitrogen, Carlsbad, CA, United States) was used in accordance with the manufacturer’s instructions, and reverse transcription (RT) was performed using the QuantiTect Reverse Transcription Kit. The cDNA samples were used for gene expression quantification via real-time PCR with a 2x Super SYBR Green qPCR Master Mix purchased from ESscience (ES-QP002, Shanghai, China), and RT-PCR was performed using an Applied Biosystems QuantStudio3 (Applied Biosystems). Each gene was run in three technical replicates and normalized to Ppia, and the fold-change relative to the sham group was calculated using the 2^−ΔΔCt formula. Primers for RT-PCR were designed using Primer-BLAST (NCBI).

Each analysis was performed in at least three independent experiments. All data are presented as the mean ± SEM (standard error of the mean). Differences between two groups were determined by two-tailed Student’s t-test using the Prism five software (GraphPad, La Jolla, CA, United States). Two-way analysis of variance (ANOVA) was used for analysis between the three groups. Statistical significance was set at p < 0.05.

To identify the differentially expressed genes at various times after renal I/R in mice, we performed RNA sequencing to explore the transcriptional progressive patterns from mild to severe kidney injury using different ischemia-time treatments. RNA sequencing data were deposited into the Gene Expression Omnibus database (GSE192883). Using differential gene analysis, we found that compared with the sham group, there were 1,156 differentially expressed genes (DEGs; 571 upregulated and 585 downregulated) in the 18-min ischemia group, 1,256 DEGs (635 upregulated and 621 downregulated) in the 22-min ischemia group, 1,325 DEGs (674 upregulated) in the 26-min ischemia group, and 651 DEGs (downregulated) and 2,119 DEGs (1,063 upregulated and 1,056 downregulated) in the 30-min ischemia group (Figure 1A). Among them, the top 10 DEGs with the greatest differences were more than 9-fold upregulated or downregulated (Table 1). Combined with principal component analysis (Figure 1B), it can be concluded that a large number of DEGs exist not only in severe renal I/R injury but also in mild renal I/R injury, indicating that many biological processes have been initiated in mild renal I/R injury and continue to change with the aggravation of the injury.

Among the top 50 DEGs in each group, the proportion of overlapping genes between the 18-min and 22-min ischemia groups was 67.1%; the proportion of overlapping genes among the 18-min, 22-min, and 26-min ischemia groups was 64.3%; and the proportion of overlapping genes among the four experimental groups was 42.9%. This suggests that DEGs gradually changed with the aggravation of injury but the most different genes changed unsustainably and discontinuously as time changes (Figure 1C). The heat map of overlapping DEGs among the four experimental groups shows that the expression levels of these DEGs changed regularly as the injury intensified (Figure 1D), supporting the results of the discrete time points analysis.

During the time-series analysis, all genes were automated and clustered into eight groups. Cluster1 indicated that the change in expression level was constantly accompanied by injury, which was a gene set related to I/R injury (Figure 1E). Cluster2 and Cluster3 showed gene sets that changed only in severe injury and represented the gene sets related to severe I/R injury (Figure 1F). Cluster4 to Cluster8 gene sets relate to the progression from mild to severe I/R kidney injury and showed a trend of continuous upregulation (Cluster4 and Cluster5), upregulation followed by downregulation (Cluster6 and Cluster7), and continuous downregulation (Cluster8). Moreover, the 18-min and 30-min ischemia groups were the two time points with more obvious changes, suggesting progressive patterns from mild to severe renal I/R injury (Figure 1G). Phenotypically, the injury initially occurred in the IR30 group but not the IR18, IR22 or IR26 group suggests that the transcriptome in the IR18 group has undergone changes, while they are not reflected in the pathological sections (Supplementary Figure S4). The results of the Cluster1–8 Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) functional enrichment analyses are presented in Table 2.

The core pathway of Cluster4 was identified as “Ribonucleoprotein complex biogenesis” through the overlapping pathways and pathway network analysis (Figure 2A) of the 18-min and 30-min ischemia groups; the GSEA result (Figure 2B) showed that the activation level of this pathway in the 30-min ischemia group was higher than that in the 18-min ischemia group. This suggests that the biogenesis function of the ribonucleoprotein complex was activated during the transition from mild to severe I/R kidney injury. The heat map of core enrichment genes, which contributed the most to the enrichment results, showed the same transitional pattern (Figure 2C). The core pathways of Cluster5 and Cluster8 were identified as “Cell-substrate junction organization” and “Organic acid catabolic process” through the same method as Cluster4 (Figures 2D,G); GSEA results (Figures 2E,H) showed that the activation level of “Cell-substrate junction organization” in the 30-min ischemia group was higher than that in the 18-min ischemia group and the inhibition level of “Organic acid catabolic process” in the 30-min ischemia group was higher than that in the 18-min ischemia group. This indicates that the adhesive function between the cell and substrate was activated and that the organic acid catabolic process was inhibited during the transition from mild to severe I/R kidney injury.

The heat map of the core enrichment genes showed the same transitional pattern (Figures 2F,I). The core pathways of Cluster6 and Cluster7 were identified as “Chromosome segregation” and “DNA repair” through the same method (Figures 3A,D), and the GSEA results (Figures 3B,E) showed that there were no significant differences between the activation level of these two pathways in the 18-min and the 30-min ischemia groups; however, in the leading edge subset, the pathway activation level of the 18-min group was significantly higher than that of the 30-min group. The heat map shows a similar transitional pattern (Figures 3C,F). The other representative pathways of the five clusters with significant differences, as well as the representative pathways of other clusters, are also shown (Figure 3G, Supplementary Figures S1, S2). Genes with a limit fold-change (LFC) greater than 0.7 in the IR18 group, LFC greater than one in the IR30 group or LFC greater than one in the IR18 group, and LFC greater than 0.7 in the IR30 group were screened out in the leading edge subset of each core pathway as key genes (Supplementary Figure S3), and quantitative polymerase chain reaction (qPCR) verification was conducted. The results showed that S100a10 was continuously upregulated in the IR18 and IR30 groups, which is consistent with the sequencing results. Pcna was upregulated in the IR18 group and downregulated in the IR30 group. Abat, Kmo, Acadm, and Adhfe1 were continuously downregulated (Figure 4B), which is consistent with the sequencing results for these molecules (Figure 4A). However, there were no significant changes between IR18, IR22 and IR26 groups.

Specifically, we used a time-course (16–30 min) of mice for the I/R injury with a renal pedicle clip. Using our sequencing data, we obtained gene subnetworks (DNB) with fluctuations in different ischemic time groups according to DNB analysis methods. The CI value of the DNB subnetwork was the highest in the 22-min ischemia group, suggesting that the transcriptome expression level fluctuated sharply in this group and that 22 min of ischemia may be the key time point in the transition from mild to severe I/R kidney injury. When the ischemia time was prolonged and the injury was aggravated, the CI value of DNB began to decline, confirming that DNB only specifically recognizes the state of “pre-disease”, as previously reported (Liu et al., 2014) (Figure 5B). The molecular network diagram shows that except for the 22-min ischemia group, DNB in the other ischemic groups had strong connections with other genes, indicating that the internal network resonance of DNB in this group was strong, while its relationship with other genes was significantly weakened (Figure 5A). Detailed calculation results are shown in Table 3. However, the standard deviation and internal correlation of DNB in the 22-min ischemia group were not the highest, suggesting that the key contribution to the transition from mild to severe injury during renal I/R lies in the weakened relationship between DNB and other genes. Briefly, during the transition from mild to severe I/R kidney injury, renal transcriptome expression levels fluctuated strongly when the ischemia time reached 22 min, and the injury began to transition from mild to severe (Figure 5C). Therefore, 22 min of ischemia can be considered a warning point for severe irreversible I/R kidney injury.