Related gene expression datasets were identified using “liver transplantation” MeSH Terms AND “Homo sapiens” Organism AND Expression profiling via GEO (http://www.ncbi.nlm.nih.gov/geo) (9), a public database containing a large number of high-throughput sequencing and microarray datasets submitted by research institutions worldwide. The GSE12720 dataset contains 13 liver tissues from pre-ischemic donors and 13 liver tissues from reperfused livers after liver transplantation. The GSE14951 dataset contains five liver tissues from pre-ischemic donors and five liver tissues from reperfused livers after liver transplantation. The GSE15480 dataset contains six liver tissues from pre-ischemic donors and six liver tissues from reperfused livers after liver transplantation. Details of three datasets are shown in Table 1 . The probe names were transformed into gene symbols based on platform annotation information. Moreover, immune-related genes were obtained from the Immunology Database and Analysis Portal (IMMPORT) database (http://www.immport.org/) (10).

The three raw datasets were downloaded using the R “GEOquery” package, and then performed gene annotation, gene filtering, normalization, etc. on the three datasets. The differentially expressed genes (DEGs) were analyzed and compared between the pre-ischemic and reperfused liver tissue controls using the R “limma” package. Only those genes with adjusted P-value < 0.05 and |log2FC (fold change)| ≥ 0.5 were considered DEGs. The R “ggplot2” package was used to draw volcano plots and the R “pheatmap” package was used to draw heatmaps for data visualization.

The venn diagram was used to analyze common DEGs. Intersections of the DEGs from the three datasets, as well as the intersection between hub genes obtained from CytoHubba and immune-related genes were identified using the Venny online tool, version 2.1 (https://bioinfogp.cnb.csic.es/tools/venny/index.html).

In order to explore the functions of common DEGs between the three datasets and the involved modular genes obtained using molecular complex detection technology (MCODE), gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using the R “org.Hs.eg.db” and “ClusterProfiler” packages. P < 0.05 was considered statistically significant. The R “ggplot2” package was used to visualize the enrichment results.

The STRING database (https://string-db.org/) (11) was used to identify the relationships between DEGs that were present in all three databases, including direct binding relationships, or coexisting upstream and downstream regulatory pathways, and a PPI network with complex regulatory relationships was created. Interactions with a combined score of > 0.4 were considered statistically significant. The nodes represented genes, and the edges represented the links between the genes. Cytoscape software (https://cytoscape.org) (12) was used to develop the PPI network.

Cytoscape’s plug-in molecular complex detection technology (MCODE) (13) was used to analyze key functional modules. The function of MCODE is to select key sub-networks, i.e., modules. The selection criteria were set as: K-core = 2, degree cutoff = 2, max depth = 100, and node score cutoff = 0.2. GO and KEGG analyses of the involved modular genes were then conducted using the R “org.Hs.eg.db” and “ClusterProfiler” packages.

The hub genes were identified using the cytoHubba plug-in of Cytoscape. TheDegree algorithm was used to evaluate and select hub genes. And immune-related hub genes as candidates for subsequent analysis. A co-expression network of the immune-related hub genes was conducted using GeneMANIA (http://www.genemania.org/) (14), a reliable tool for identifying internal associations in gene sets.

Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST) (15) is a database used to predict transcriptional regulatory networks, which contain the target genes corresponding to TFs and the regulatory relationships between TFs. The TRRUST database includes 8,444 and 6,552 TFs target regulatory relationships from 800 human TFs and 828 mouse TFs, respectively. TRRUST database was used to obtain the TFs that regulate immune-related hub genes, and an adjusted P-value < 0.05 was considered significant.

Potential immune-related hub genes targeted miRNAs were predicted using the TargetScan (http://www.targetscan.org/vert_72/) (16), miRTarBase (https://mirtarbase.cuhk.edu.cn/php/index.php) (17), and miTDB (http://mirdb.org/) (18) databases. The miRNAs for each immune-related hub gene predicted by the three databases were taken to intersect. The results obtained from the intersection were further processed usingCytoscape, miRNAs targeting more than two genes were selected. The upstream LncRNA of the obtained miRNAs was predicted using StarBase v2.0 (19), and the interactions were visualized using Cytoscape

Twenty 6-8 week-old male C57BL/6J mice weighing 22–25 g were acquired from the Center of Experimental Animals at the Renmin hospital of Wuhan University, Hubei, China. The animal experiments were conducted according to the Guide for the Care and Use of Laboratory Animals, and the study protocol was approved by the Laboratory Animal Welfare and Ethics Committee of the Renmin hospital of Wuhan University. Mice were kept in an air-filtered, temperature-controlled (22–24°C), humidity-controlled (40-70%), and light-controlled room and were permitted free access to a standard diet.

The stable mouse model of partial liver warm ischemia was constructed according to the classical method, as previously described. The mice were anaesthetized with 3% pentobarbital sodium prior to surgery. Sham-operated mice underwent the same surgical procedure without clamping of the blood vessels. The 70% liver tissue ischemia model, in which the portal vein and hepatic artery of the middle and left lobes areblocked with non-invasive vascular clamps, was used during surgery. Ischemia was performed for 60 minutes, and reperfusion times varied. Animals were euthanized post-procedure, and the serum and liver samples were immediately collected for further analyses.

The concentrations of serum alanine aminotransferase (ALT) and serum aspartate aminotransferase (AST) were measured using the ADVIA 2400 Biochemistry Analyzer (Siemens, Tarrytown, NY, USA).

After fixation with 4% paraformaldehyde, liver samples were embedded in paraffin and cut into 5-μm sections. H&E staining was used to assess histopathological liver damage and quantified using Image J software. The percentage of necrotic area in the total area of the tissue section was quantified in more than five fields for each mouse.

An apoptosis assay kit (ApopTag^® Plus In Situ Apoptosis Fluorescein Detection Kit, S7111; Millipore, USA) was used to rapidly detect apoptosis induced by I/R injury. TUNEL assay was performed according to the manufacturer’s instructions. Images were obtained using fluorescent microscopy. Quantification of TUNEL-positive cells in each field of the tissue section was performed using Image J software

Immunohistochemical and immunofluorescence staining was conducted to analyze the protein expression levels of iNOS, CD11b, and CD206. All procedures were performed according to the manufacturer’s recommendations. Positive areas were compared between groups using microscopy and images were analyzed with Image J software.

The total mRNA was extracted from cultured cells using TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, United States). RNAs were reverse transcribed into cDNA using ImProm-II™ Reverse Transcription System (Promega, USA) according to the instructions of the manufacturer. Quantitative RT-PCR with cDNA as template and β-Actin as internal reference to analyze the relative mRNA expression of each gene. The primer sequences used for each gene have been listed in Supplementary Table S1 . The 2^−ΔΔCT method was used to calculated the relative mRNA expression levels which were normalized against β-Actin.

All data are expressed as the mean ± standard deviation. SPSS version 19.0 was used for statistical analysis. The student’s t-test was used to compare differences between groups. Differences were considered significant when the P value was <0.05. All experiments were performed at least three times.

To clarify the process of this research, the research flowchart is shown in Figure 1 . Original data were downloaded from the GSE12720, GSE14951 and GSE15480 datasets in the GEO database. After standardizing the microarray data, differential expression analysis was performed and DEGs were identified. In the GSE12720 dataset, the transcriptome data from 13 pre-ischemic donor liver tissues and 13 liver tissues after reperfusion were analyzed using criteria of |log2FC| ≥ 0.5 and P < 0.05, yielding 192 DEGs (179 upregulated and 13 downregulated). In the GSE14951 dataset, transcriptome data from five pre-ischemic donor liver tissues and five liver tissues after reperfusion were analyzed using criteria of |log2FC| ≥ 0.5 and P < 0.05, yielding 372 DEGs (365 upregulated and 7 downregulated). In the GSE15480 dataset, the transcriptome data from six pre-ischemic donor liver tissues and six liver tissues after reperfusion were analyzed using criteria of |log2FC| ≥ 0.5 and P < 0.05, yielding 353 DEGs (212 upregulated and 141 downregulated). The corresponding heatmaps and volcano plots are shown in Figure 2 . Details of the three datasets are presented in Table 1 .

A Venn diagram was created to determine the intersection between the DEGs from the three datasets and 71 common DEGs were identified ( Figure 3A ). To analyze the biological functions and involed pathways, GO and KEGG enrichment analysis were conducted. GO analysis showed that these genes were mostly enriched in the ERK1 and ERK2 cascade (p = 1.95E-08), response to lipopolysaccharide (p = 2.90E-08), response to interleukin-1 (p = 3.76E-08), regulation of transcription from the RNA polymerase II promoter in response to stress (p = 6.38E-07), and DNA-binding transcription activator activity (p = 3.40E-12) ( Figure 3B ). KEGG pathway enrichment analysis showed that the three significant enrichment pathways were IL-17 signaling pathway (p = 1.26E-08), TNF signaling pathway (p = 5.07E-08), NF-kappa B signaling pathway (p = 9.80E-06), Toll-like receptor signaling pathway (p = 1.3E-04), NOD-like receptor signaling pathway (p = 2.4E-04), and Human T-cell leukemia virus 1 infection (p = 6.6E-04) ( Figure 3C ). These results strongly indicated that immune response, inflammatory response and injury repair actions are involved in the development and occurrence of hepatic I/R injury.

In order to explore the mutual interaction between common DEGs, the PPI network of common DEGs with a combined score of > 0.4, which contained 52 nodes and 218 interaction pairs, was constructed using Cytoscape ( Figure 4A ). Two closely related gene modules were obtained by MCODE plug-in of Cytoscape, including 25 common DEGs and 73 pairs of interaction ( Figures 4B, C and Table 2 ). The genes from the obtained modules were then subjected to functional enrichment analysis. GO analysis showed that these genes are associated with response to lipopolysaccharide, cellular response to external stimulus, transcriptional activation of DNA, and chemokine activity ( Figure 4D ). KEGG pathway analysis showed that these genes were primarily involved in the response to IL-17 signaling pathway, TNF signaling pathway, NF-kappa B signaling pathway, and NOD−like receptor signaling pathway ( Figure 4E ).

Using the Degree algorithm from the plug-in cytoHubba of Cytoscape, the top 20 hub genes (SOCS3, JUNB, NFKBIA, JUND, FOSL2, FOSB, FOSL1, TNFAIP3, CCL4, JUN, GADD45B, ATF3, CDKN1A, KLF6, MYC, KLF4, CEBPB, ICAM1, CXCL8, IRF1) were calculated ( Figure 5A ). To explore immune-associated hub genes involved in hepatic I/R injury, immune-associated genes were obtained from the IMMPORT database. Subsequently, the Venn method was used to analyze the intersection between hub genes obtained in cytoHubba and immune-associated genes.Nine overlapping immune-associated hub genes (SOCS3, JUND, CCL4, NFKBIA, CXCL8, ICAM1, IRF1, TNFAIP3, JUN) were identified ( Figure 5B ). Based on the Gene MANIA database, the co-expression network and associated functions of these genes were assessed. The genes had a complex PPI network with a 78.61% co-expression, 8.08% physical interaction rate, 6.42% co-localization rate, 1.95% predicted rate, and 0.35% pathway rate ( Figure 5C ).

Twelve TFs were found to regulate these immune-related hub genes based on the TRRUST database ( Figure 6 ). Detailed information of the 12 TFs is shown in Table 3 . At the same time, we found that the Hub genes JUN and NFKBIA are also upstream transcription factors regulating these Hub genes. NFKBIA may be regulated by the transcription factors, NFKB1 and RELA, and may also act as a transcription factor that regulates ICAM1 and CXCL8. JUN acts as a transcription factor that regulates CXCL8.

To identify the upstream miRNAs of these hub genes and ensure the accuracy and reliability of the results, we took intersections of the miRNAs of these genes predicted by each of the three databases. A total of 161 predicted upstream miRNAs were obtained, of which no targeted miRNAs were identified for CCL4 ( Figure 7 ). The miRNAs with a high number of cross-linked genes (≥2) are shown in Table 4 . Eleven miRNAs had a higher number of cross-linked genes (≥2), of which only hsa-miR-19b-3p, hsa-miR-19a-3p, has-miR-23a-3p, has-miR-23c had the upstream molecule LncRNA and had higher reliability, all targeting TNFAIP3. The lncRNAs corresponding to hsa-miR-19b-3p, hsa-miR-19a-3p, has-miR-23a-3p, has-miR-23c were predicted using StarBase 2.0, and the highest reliability (very high stringency, >5) was selected as the criterion. After cross-linking the lncRNAs and miRNAs, three lncRNAs targeting four key miRNAs, ZNF518A, RP11-170L3.8, and XIST, were identified ( Figure 8 ).

The effect of hepatic I/R injury after clamping of blood vessels was verified, and the extent of apoptosis and associated immune cell infiltration was assessed. After experiencing 1 hour of ischemia, the left and middle lobe livers of mice were reperfused. As the early reperfusion time was lengthened, HE staining of the livers of mice revealed increasingly severe liver necrosis (12h>6h). Hepatocyte apoptosis and levels of the blood biochemical indexs, ALT and AST, also increased with longer reperfusion times ( Figures 9A–D ). Immunohistochemical staining analysis of the macrophage-associated molecule INOS (M1) and the neutrophil-associated molecule CD11b showed a significant increase in M1 macrophage and neutrophil infiltration at longer reperfusion times. The immunofluorescence staining analysis of the macrophage-associated molecule CD206 (M2) also showed a obvious increase in M2 macrophage infiltration after reperfusion ( Figures 9A, B ). Expression of these immune-related hub genes was validated in the three databases, indicating that they were all significantly upregulated after hepatic I/R injury( Figures 10A–C ). A mouse model of hepatic I/R injury was used to verify the expression levels of immune-related hub genes by qRT-PCR. As shown in Figure 10D , the expression levels of SOCS3, CCL4, NFKBIA, CXCL8, ICAM1, IRF1, and TNFAIP3 were significantly upregulated after 6 hours of I/R, consistent with the findings of GEO datasets. In contrast, the expression levels of JUN and JUND did not differ, which may be related to species, reperfusion time. A number of immune signalling downstream molecules were validated and TNF-α, IL-1β, IL-17, NF-κB-p65 and TLR-4 were significantly upregulated after 6h of hepatic ischemia-reperfusion, further confirming the involvement of innate immune signaling identified in the study in the process of hepatic ischemia-reperfusion injury ( Figure 10E ).