Gene expression profiling of sepsis-associated acute kidney injury.

Sepsis accounts for more than 50% of all acute kidney injury (AKI) cases, and the combination of sepsis and AKI increases the risk of mortality from sepsis alone. However, to the best of our knowledge, the specific mechanism by which sepsis causes AKI has not yet been fully elucidated, and there is no targeted therapy for sepsis-associated AKI (SA-AKI). The present study investigated gene expression profiles using RNA sequencing (RNA-Seq) and bioinformatics analyses to assess the function of differentially expressed genes (DEGs) and the molecular mechanisms relevant to the prognosis of SA-AKI. From the bioinformatics analysis, 2,256 downregulated and 3,146 upregulated genes were identified (false discovery rate <0.1 and fold-change >2). Gene Ontology analysis revealed that the genes were enriched in cellular metabolic processes, cell death and apoptosis. The enriched transcription factors were v-rel reticuloendotheliosis viral oncogene homolog A and signaling transducer and activator of transcription 3. The enriched microRNAs (miRNAs or miRs) among the DEGs were miR-30e, miR-181a, miR-340, miR-466d and miR-466l. Furthermore, the enriched pathways included toll-like receptor signaling, nod-like receptor signaling and the Janus kinase/STAT signaling pathway. In conclusion, the present study identified certain prognosis-related genes, transcription factors, miRNAs and pathways by analyzing gene expression profiles of SA-AKI using RNA-Seq, which provides some basis for future experimental studies. Copyright: © Yang et al.

Introduction

Sepsis is a clinical syndrome characterized by organ dysfunction caused by a dysregulated immune response to infection (1). Acute kidney injury (AKI) is a serious clinical problem with high levels of mortality (2). Sepsis accounts for more than 50% of all AKI cases, and the combination of sepsis and AKI increases the risk of mortality compared with sepsis alone (3,4). However, to the best of our knowledge, the specific mechanism by which sepsis causes AKI has not yet been fully elucidated, and there is no targeted therapy for sepsis-associated AKI (SA-AKI).

Although the pathophysiological mechanism has not yet been elucidated, it has been suggested that the deleterious inflammatory cascade characteristic of sepsis contributes to AKI (5). Compared with patients with sepsis without AKI, those with AKI have a higher risk of mortality (6,7). However, the renal function of some patients with SA-AKI can return to a normal level (8). A previous study reported that resuscitation of lipopolysaccharide-treated mice with warmed 0.9% saline can increase survival rates to >80%, which allows to study the prognosis of SA-AKI (9). It is therefore vital to investigate the heterogeneity and molecular mechanisms related to the prognosis of SA-AKI, which may provide new therapeutic targets for SA-AKI.

RNA sequencing (RNA-Seq) has become a routine technique for genome-wide expression analysis as it can provide high-resolution sequence information (10). The present study investigated gene expression profiles using RNA-Seq and subsequent bioinformatics analyses to assess the function of differentially expressed genes (DEGs) and molecular mechanisms relevant to the prognosis of SA-AKI.

Materials and methods

Mouse cecal ligation and puncture (CLP) model

A total of 20 8-week-old male C57BL/6J mice (weight, 20-25 g) were obtained from Shanghai SLAC Laboratory Animal Co., Ltd., and housed in a room at 22±1˚C, 40-60% humidity, with a 12-h light/dark cycle (8:00 a.m.-8:00 p.m.), with free access to food and drinking water. All operations on mice were performed by one individual who was skilled in animal modeling. As described previously (9), a 1-cm midline laparotomy was performed following anesthesia with 3% inhaled isoflurane. The base of the cecum was ligated and punctured at its distal tip using a single 21-gauge needle. Subsequently, a small amount of stool content (about 1 mm length) was extruded and put back into the abdominal cavity. The abdomen was sealed and 1 ml warm sterile 0.9% saline was injected subcutaneously to resuscitate mice. The sham group (n=8) underwent laparotomy and cecum exposure but no cecal ligation or puncture.

At 24 h or 48 h post-surgery, all mice were anesthetized by inhalation of 3% isoflurane, and then 50 µl retro-orbital venous blood was collected to measure Scr and BUN levels. Following collection of retro-orbital venous blood to assess Scr and BUN levels, 4 mice were sacrificed (24 h post-CLP, SA-AKI group) and 8 mice were given an additional 10 ml/kg warmed sterile 0.9% saline, and then sacrificed at 48 h post-CLP by cervical dislocation. AKI was defined as a serum creatinine (Scr) level >20 µmol/l and blood urea nitrogen (BUN) >28.3 mg/dl. Mice were grouped into persistent injury (PI; n=3) and recovery (RE; n=5) groups according to renal function at 48 h post-CLP. The RE group consisted of mice whose Scr and BUN levels returned to normal (<20 µmol/l and <28.3 mg/dl, respectively) at 48 h. The PI group contained mice with high levels of Scr (>20 µmol/l) and BUN (>28.3 mg/dl) at 48 h. Mice in the SA-AKI group were sacrificed to collect blood and kidneys at 24 h (n=4). Mice in the PI (n=3) and RE (n=5) group were sacrificed to collect blood and kidneys at 48 h. Mice in the sham group were sacrificed via cervical dislocation to collect blood and kidneys at 24 h (n=4) and 48 h (n=4), respectively.

All protocols of this study were approved by the Animal Ethics Committee of the First Affiliated Hospital, College of Medicine, Zhejiang University, Hangzhou, Zhejiang (approval ID. 2016160; 25 February 2016), which follows the institutional guidelines.

Renal function assessment

Blood samples were centrifuged at 13,000 x g for 10 min at 4˚C to collect the supernatant. Scr and BUN levels were measured using the DRI-CHEM dry biochemical analyzer (FUJIFILM Wako Pure Chemical Corporation).

Hematoxylin-eosin staining

The kidney tissues of mice were fixed in 4% paraformaldehyde at room temperature overnight, embedded in paraffin wax, sliced into 4-mm-thick sections, and stained with hematoxylin and eosin at 37˚C for 24 h. Slides were prepared in triplicate, and sections were evaluated under a light microscope (magnification, x400). The pathological injury score of all the samples was assessed twice in a blinded fashion by independent investigators, and the mean score from the two assessments was recorded as the composite score. The grading criteria for renal injury were as follows: Edema of tubular epithelial cells, loss of normal morphology, and cytoplasmic vacuoles. The grading standard for renal tubular involvement were as follows: Observing the aforementioned characteristics of the samples under the same microscope at x400 magnification. Renal tissue pathological injury score was assessed with a 0-5 scoring system: 0, no renal tubular involvement; 1, <25% renal tubular involvement; 2, 25-<50% renal tubular involvement; 3, 50-75% renal tubular involvement; 4, >75% renal tubular involvement.

RNA-Seq and gene expression analysis

Total RNA of 3 kidneys in the PI group and 3 kidneys in the RE group were extracted using TRlzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. RNA purity was assessed using the kaiaoK5500® Spectrophotometer (Kaiao). RNA integrity and concentration were detected using the RNA Nano 6000 assay kit of the Bioanalyzer 2100 system (Agilent Technologies, Inc.). A total amount of 2 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were built using the NEBNext® Ultra™ RNA Library Prep kit for Illumina® (cat. no. E7530L; New England Biolabs, Inc.) following the manufacturer's instructions, and index codes were added to attribute sequences to each sample.

Based on the manufacturer's instructions, the index-coded samples were clustered on a cBot cluster generation system using HiSeq PE Cluster kit v4-cBot-HS (Illumina, Inc.). After the cluster was generated, the libraries were sorted on an Illumina platform to generate 150 bp paired-end reads. Genes with a false discovery rate (FDR) <0.1 and fold-change (FC) >2 were identified as DEGs.

Data analysis and pathway enrichment analysis were performed using integrated Differential Expression and Pathway analysis (iDEP, version 0.90; http://ge-lab.org/idep) (11). iDEP is an online application that integrates many Bioconductor packages and annotation databases to enable users to perform intensive bioinformatics analysis. In the present study, DESeq2 package, hierarchical clustering and κ-means clustering analysis, principal component analysis, Gene Ontology (GO) analysis, and transcription factor binding motifs and microRNA enrichment analysis were used in iDEP application.

Gene set enrichment analysis (GSEA)

GSEA was conducted using the software GSEA v2.2.2 (www.broadinstitute.org/gsea). The enrichment degree and statistical significance were quantified by normalized enrichment score (NES), nominal P-value and FDR (12). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were conducted by using GSEA software.

Statistical analysis

All data are presented as mean ± standard deviation. The results were analyzed by one-way analysis of variance and Tukey's post hoc test using SPSS 22.0 software (IBM Corp.). P<0.05 was considered to indicate a statistically significant difference.

Results

Sepsis-associated AKI mice

CLP was performed in 8-week-old male C57BL/6J mice. At 24 h after the CLP procedure, kidney function biomarkers BUN and Scr were significantly elevated. Following resuscitation, 5 mice exhibited recovery to normal kidney function by 48 h post-CLP operation, however 3 mice had persistent AKI. The levels of kidney function biomarkers are presented in Fig. 1A. It was identified that Scr and BUN levels did not significantly change in sham group at each time point. Renal tissue pathological staining and injury score were used to evaluate renal injury in the present study. In the sham group, the kidney tissue was demonstrated to be normal and no signs of damage were present. To further quantify the degree of renal damage in each group of mice, the pathological grading of renal pathological sections was evaluated. The renal tissue pathological injury score in the SA-AKI group (3.50±0.58) was significantly higher compared with that in the sham group (0.50±0.58) and RE group (1.60±0.55; P<0.05). There was no significant difference in renal tissue pathological injury score between the SA-AKI and PI groups (3.50± 0.58 vs. 3.33± 0.58; P=0.72; Fig. 1B). To investigate this heterogeneity, expression profiling was performed with the kidney-derived RNA obtained from the PI and RE groups 48 h post-CLP.

Figure 1

Renal function and renal structures in SA-AKI mice. (A) A total of 24 h post-CLP, mice received 10 ml/kg saline, and by 48 h post-CLP, a number of mice had recovered from AKI. (B) Hematoxylin-eosin staining revealed micro-vacuolization in proximal tubular cells in SA-AKI mice. The renal tissue pathological injury score in the SA-AKI group was significantly higher compared with that in the sham and recovery groups. There was no significant difference in the renal tissue pathological injury score between the SA-AKI and persistent injury groups. Magnification, x400. All experiments were performed a minimum of three times, and the results are presented as the mean ± standard deviation. *P<0.05 vs. Sham group. #P<0.05 vs. SA-AKI group. NSP<0.05 vs. SA-AKI group. SA-AKI, sepsis-associated acute kidney injury; CLP, cecal ligation and puncture.

Expression profiling analyses of kidneys from mice with SA-AKI

iDEP was used to analyze the RNA-Seq dataset. The total read counts per library exhibited a small amount of variation in size (Fig. 2A). Regularized log transformation was performed using the DESeq2 package, and the distribution of the transformed data is presented in Fig. 2B and C, D presented the variation between the RE 1 and RE 2 samples to represent the variation among each replicate. Variation among each replicate was small.

Figure 2

Diagnostic plots for read-counts data. (A) Total read-counts in the RE and PI samples. (B) Distribution of transformed data as presented using a density plot. (C) Boxplot of transformed data. (D) Scatter plot of the first two samples: RE_1 and RE_2. RE, recovery; PI, persistent injury.

Hierarchical clustering and κ-means clustering analysis

Hierarchical clustering analysis evaluated the top 2,000 significant DEGs among all DEGs, including the 2256 downregulated and 3146 upregulated DEGs, which were ranked by their standard deviation across all samples. The results of hierarchical clustering analysis results are presented in Fig. 3A, in which each column represents the expression pattern of one sample, the red lines represent the upregulated genes, and the green lines represent the downregulated genes (Fig. 3A). A scatter plot was generated to depict the gene expression distribution (Fig. 3B), and a volcano map was created to visualize significant DEGs (Fig. 3C). A principal component analysis (PCA) plot using the first and second principal components is presented in Fig. 3D. There was a significant difference between the RE and PI groups along the first principal component, which explained the 62% variance.

Figure 3

Bioinformatics analysis of mRNA expression patterns in the PI and RE groups by RNA-Seq. (A) Hierarchical clustering was used to evaluate the top 2,000 genes with significant differential expression in the PI and RE group samples. Each column represents one samples expression pattern. The red lines represent the upregulated genes, and the green lines represent the downregulated genes. (B) Scatter plot analysis was conducted to demonstrate the mRNA expression distribution. (C) The volcano plot was generated to visualize the significant differentially expressed genes. (D) Principle component analysis demonstrated that there was a significant difference between the RE and PI group samples, along the first principal component that explained 62% of the variance. *P<0.05 vs. RE group; RE, recovery; PI, persistent injury; FDR, false discovery rate.

Subsequently, κ-means clustering was used to divide the top 2,000 DEGs into groups. Based on the within-group sum of squares plot as a reference, κ=4 was selected and indicated the number of clusters used for further analysis (Fig. 4A). The four gene clusters underwent GO analysis and the results are presented in Fig. 4B and Table I.

Figure 4

κ-means clustering and GO enrichment analysis of the RE and PI group samples. (A) Based on the within-group sum of squares plot as a reference, a larger κ=4 was selected. (B) The clusters A, B, C and D were predominantly enriched in cell development, metabolic process, regulation of apoptosis and immune response, and defense response, respectively. RE, recovery; PI, persistent injury; GO, Gene Ontology.

Table I

Enrichment analysis of k-means clusters using various gene sets.

Cluster	Adjusted P-value	Genes (n)	Pathway
A	2.03x10-5	42	GO:0009605 response to external stimulus
A	0.000133	40	GO:0048468 cell development
A	0.000311	17	GO:0006935 chemotaxis
A	0.000311	31	GO:0040011 locomotion
A	0.000555	32	GO:0022610 biological adhesion
A	0.000555	8	GO:0098754 detoxification
A	0.000613	4	GO:0015671 oxygen transport
A	0.000613	27	GO:0032989 cellular component morphogenesis
A	0.0009	31	GO:0007155 cell adhesion
A	0.001458	31	GO:0006928 movement of cell or subcellular component
A	0.001516	25	GO:0016477 cell migration
A	0.001626	24	GO:0048646 anatomical structure formation involved in morphogenesis
A	0.001668	19	GO:0040012 regulation of locomotion
B	0.000751	46	GO:0006082 organic acid metabolic process
B	0.000751	26	GO:0006631 fatty acid metabolic process
B	0.000751	33	GO:0032787 monocarboxylic acid metabolic process
B	0.000991	41	GO:0019752 carboxylic acid metabolic process
B	0.001931	7	GO:0060740 prostate gland epithelium morphogenesis
B	0.004579	8	GO:0030850 prostate gland development
B	0.004829	6	GO:0055081 anion homeostasis
B	0.00809	52	GO:0006629 lipid metabolic process
B	0.00809	5	GO:0060572 morphogenesis of an epithelial bud
C	2.48x10-46	109	GO:0034097 response to cytokine
C	5.84x10-42	174	GO:0009605 response to external stimulus
C	2.22x10-41	135	GO:0006952 defense response
C	1.04x10-40	174	GO:0002376 immune system process
C	1.71x10-38	168	GO:0071310 cellular response to organic substance
C	2.96x10-36	104	GO:0009607 response to biotic stimulus
C	2.96x10-36	101	GO:0043207 response to external biotic stimulus
C	2.96x10-36	88	GO:0071345 cellular response to cytokine stimulus
C	3.22x10-34	153	GO:0008219 cell death
C	3.57x10-34	118	GO:0006955 immune response
C	6.62x10-34	64	GO:0032496 response to lipopolysaccharide
C	7.23x10-34	147	GO:0012501 programmed cell death
C	1.98x10-33	145	GO:0006915 apoptotic process
C	3.00x10-33	148	GO:0048584 positive regulation of response to stimulus
D	6.27x10-8	25	GO:0034097 response to cytokine
D	1.24x10-8	8	GO:0035458 cellular response to interferon-beta
D	4.1x10-6	7	GO:0006953 acute-phase response
D	5.58x10-5	27	GO:0006952 defense response
D	5.77x10-5	18	GO:0071345 cellular response to cytokine stimulus
D	7.70x10-5	35	GO:0009605 response to external stimulus
D	0.00017	8	GO:0034341 response to interferon-gamma
D	0.000215	14	GO:0051346 negative regulation of hydrolase activity
D	0.000269	6	GO:0050994 regulation of lipid catabolic process
D	0.000509	13	GO:0052548 regulation of endopeptidase activity
D	0.000515	27	GO:1901700 response to oxygen-containing compound
D	0.00062	3	GO:0010727 negative regulation of hydrogen peroxide metabolic process
D	0.00062	4	GO:0035634 response to stilbenoid

GO, Gene Ontology.

Identification of DEGs

Compared with the PI group, a total of 2,256 downregulated and 3,146 upregulated DEGs were identified in the RE group using a threshold of FDR <0.1 and FC >2 (Fig. 5A). The upregulated and downregulated DEGs then underwent GO analysis for Biological Process terms. In Fig. 5B, the upregulated DEGs were mainly involved in ‘metabolic processes’, the downregulated DEGs were enriched in ‘regulation of response to stress’, ‘cytokine production’, ‘apoptotic process’ and ‘cell death’. The details of the enriched GO Biological Process terms are presented in Table Ⅱ. For GO Cellular Component terms, it was identified that the upregulated DEGs in the RE group were enriched in the ‘mitochondrion’, ‘cytoskeleton’ and ‘centrosome’; the downregulated DEGs in the RE group were enriched in the ‘cell surface’, ‘cell junction’ and ‘adherens junction’ (Fig. 5C and Table III). For GO Molecular Function terms, it was identified the downregulated DEGs in the RE group were enriched in ‘enzyme binding’ and ‘receptor binding’. Furthermore, the upregulated DEGs were enriched in ‘transporter activity and active transmembrane transporter activity’ (Fig. 5D and Table IV).

Figure 5

DEGs and enrichment analysis. (A) A total of 2,256 downregulated and 3,146 upregulated DEGs were identified in the RE group compared with the PI group using a threshold of false discovery rate <0.1 and fold-change >2. (B) The significantly enriched GO Biological Process terms for the upregulated and downregulated DEGs. (C) The significantly enriched GO Cellular Component terms for the upregulated and downregulated DEGs. (D) The significantly enriched GO Molecular Function terms for the upregulated and downregulated DEGs. DEGs, differentially expressed genes; GO, Gene Ontology; RE, recovery; PI, persistent injury.

Table III

Enriched GO Cellular Component terms for the upregulated and downregulated genes.

A, Downregulated genes
Adjusted P-value	Genes (n)	Cellular component
<0.001	159	GO:0005730 nucleolus
<0.001	115	GO:0005912 adherens junction
<0.001	116	GO:0070161 anchoring junction
<0.001	71	GO:0005925 focal adhesion
<0.001	72	GO:0030055 cell-substrate junction
<0.001	24	GO:0030684 preribosome
<0.001	180	GO:0030054 cell junction
<0.001	59	GO:0005913 cell-cell adherens junction
<0.001	98	GO:0005911 cell-cell junction
<0.001	108	GO:0009986 cell surface
<0.001	36	GO:0044452 nucleolar part
<0.001	13	GO:0032040 small-subunit processome
<0.001	69	GO:0098552 side of membrane
<0.001	178	GO:0005615 extracellular space
B, Upregulated genes
Adjusted P-value	Genes (n)	Cellular component
<0.001	104	GO:0005929 cilium
<0.001	276	GO:0042995 cell projection
<0.001	72	GO:0044441 ciliary part
<0.001	35	GO:0036064 ciliary basal body
<0.001	157	GO:0044463 cell projection part
<0.001	176	GO:0015630 microtubule cytoskeleton
<0.001	116	GO:0005815 microtubule organizing center
<0.001	89	GO:0005813 centrosome
<0.001	275	GO:0005856 cytoskeleton
<0.001	46	GO:0016323 basolateral plasma membrane
<0.001	243	GO:0005739 mitochondrion
<0.001	141	GO:0098590 plasma membrane region
<0.001	28	GO:0005814 centriole
<0.001	13	GO:0030990 intraciliary transport particle
<0.001	34	GO:0044450 microtubule organizing center part

GO, Gene Ontology.

Table IV

Enriched GO Molecular Function terms for the upregulated and downregulated genes.

A, Downregulated genes
Adjusted P-value	Genes (n)	Molecular function
<0.001	266	GO:0019899 enzyme binding
<0.001	213	GO:0005102 receptor binding
<0.001	204	GO:0042802 identical protein binding
<0.001	100	GO:0008134 transcription factor binding
<0.001	109	GO:0019900 kinase binding
<0.001	133	GO:0044212 transcription regulatory region DNA binding
<0.001	87	GO:0050839 cell adhesion molecule binding
<0.001	256	GO:0032553 ribonucleotide binding
<0.001	99	GO:0019901 protein kinase binding
<0.001	254	GO:0017076 purine nucleotide binding
<0.001	252	GO:0032555 purine ribonucleotide binding
<0.001	247	GO:0001883 purine nucleoside binding
B, Upregulated genes
Adjusted P-value	Genes (n)	Molecular function
<0.001	204	GO:0005215 transporter activity
<0.001	56	GO:0015291 secondary active transmembrane transporter activity
<0.001	78	GO:0022804 active transmembrane transporter activity
<0.001	165	GO:0022857 transmembrane transporter activity
<0.001	145	GO:0022891 substrate-specific transmembrane transporter activity
<0.001	166	GO:0022892 substrate-specific transporter activity
<0.001	13	GO:1901618 organic hydroxy compound transmembrane transporter activity
<0.001	33	GO:0015293 symporter activity
<0.001	18	GO:0016712 oxidoreductase activity acting on paired donors with incorporation or reduction of molecular oxygen reduced flavin or flavoprotein as one donor and incorporation of one atom of oxygen
<0.001	124	GO:0016491 oxidoreductase activity
0.001	27	GO:0005342 organic acid transmembrane transporter activity
0.001	118	GO:0015075 ion transmembrane transporter activity
0.001	24	GO:0015294 solute:cation symporter activity
0.002	33	GO:0004497 monooxygenase activity
0.002	37	GO:0016829 lyase activity

GO, Gene Ontology.

According to the enriched genes in the upregulated genes of the RE group, it was identified that solute carrier family 7 member 7 (Slc7a7) was highly enriched in the top six pathways (Table SI).

Transcription factor (TF) binding motifs and microRNA (miRNA or miR) target gene sets

As presented in Table Ⅴ, using the TF target gene sets in the enrichment analysis, it was detected that target genes of v-rel reticuloendotheliosis viral oncogene homolog A (RELA; FDR <1.92x10-8) and signaling transducer and activator of transcription 3 (STAT3; FDR <4.57x10-7) were enriched for the downregulated genes. After analysis, the upregulated genes did not enrich any TF. The present study then analyzed the enriched TF binding motifs in the promoters of DEGs. DEG promoters were revealed to contain a number of G-rich motifs, and would be bound by activator protein-2 (AP-2), and other factors such as basic helix-loop-helix protein (bHLH) and Cys2-His2 zinc finger (C2H2 ZF; Table VI).

Table VI

TF motifs enriched in gene promoters (300 bp) of the upregulated or downregulated genes.

A, Downregulated genes
Motif	TF	TF family	FDR
GCCTCAGG	Tcfap2a	AP-2	<0.001
TCGCCTCAGG	Tcfap2b	AP-2	<0.001
GCCCGAGGC	Tcfap2c	AP-2	<0.001
GCCTGAGG	Tcfap2e	AP-2	<0.001
CACGCG	Hes1	bHLH	<0.001
CGCGTG	Hes7	bHLH	<0.001
CGTGC	Sohlh2	bHLH	<0.001
CACGTG	Tcfl5	bHLH	<0.001
GGGGGCGG	Sp1	C2H2 ZF	<0.001
GGGGGGTC	Glis2	C2H2 ZF	<0.001
GGGGC	Plagl1	C2H2 ZF	<0.001
GGGGGCGG	Sp4	C2H2 ZF	<0.001
GGGGGCGG	Sp4	C2H2 ZF	<0.001
GGCC	Zfp711	C2H2 ZF	<0.001
GGGG	Zfp202	C2H2 ZF	<0.001
GGGCGTG	Klf7	C2H2 ZF	<0.001
GGGGGC	Zbtb7b	C2H2 ZF	<0.001
CACAGCGGG	Zic1	C2H2 ZF	<0.001
TGCGGG	Zbtb1	C2H2 ZF	<0.001
CGTGGGCG	Egr3	C2H2 ZF	<0.001
B, Upregulated genes
Motif	TF	TF family	FDR
GGGGGCGG	Sp1	C2H2 ZF	<0.001
GGGGGCGG	Sp4	C2H2 ZF	<0.001
GGGGGCGG	Sp4	C2H2 ZF	<0.001
GGGCGTG	Klf7	C2H2 ZF	<0.001
TGCGGG	Zbtb1	C2H2 ZF	<0.001
GGGCG	Klf8	C2H2 ZF	<0.001
GGGGGGG	Zfp740	C2H2 ZF	<0.001
GGGGGG	Zfp740	C2H2 ZF	<0.001
CGCGC	Zfp161	C2H2 ZF	<0.001
CG	Cxxc1	CxxC	<0.001
CG	Kdm2b	CxxC	<0.001
GCGC	E2f3	E2F	<0.001
GGCGC	E2f2	E2F	<0.001
GTGGGGGCGGGAG	E2f3	E2F	<0.001
GGGGGCGGGGC	Sp2	C2H2 ZF	<0.001
GGGCGGGGC	Klf5	C2H2 ZF	<0.001
GGGGGGGGGCC	Patz1	C2H2 ZF	<0.001
CACAGCGGGGGGTC	Zic4	C2H2 ZF	<0.001
GTGGGGGGG	Zfp740	C2H2 ZF	<0.001
CACAGCGGGGGGTC	Zic3	C2H2 ZF	<0.001

TF, transcription factor; FDR, false discovery rate.

The predicted target genes for numerous miRNAs were then analyzed. As presented in Table Ⅶ, it was identified that miR-30e-5p, miR-181a-5p and miR-340-5p were significantly enriched target genes for the downregulated genes. In addition, miR-466d-5p and miR-466l-5p were significantly enriched target genes for the upregulated DEGs.

GSEA

To detect KEGG pathways that were differentially enriched in the PI and RE groups, GSEA was conducted between the PI and RE groups. The most significantly enriched signaling pathways were selected based on their NES (Fig. 6 and Table VIII). As presented in Fig. 6, ‘toll like receptor (TLR) signaling’, ‘nod like receptor (NLR) signaling’, ‘pathogenic Escherichia coli infection’, ‘cytokine-cytokine receptor interaction’ and ‘Janus kinase (JAK)/STAT signaling’ were differentially enriched in the PI group. Furthermore, ‘aminoacyl tRNA biosynthesis’, ‘valine leucine and isoleucine degradation’, ‘one carbon pool by folate’ and ‘oxidative phosphorylation’ were differentially enriched in the RE group.

Figure 6

Enrichment plots from GSEA. GSEA results showing that (A) ‘Toll-like receptor signaling’, (B) ‘NOD-like receptor signaling’, (C) ‘Pathogenic Escherichia coli infection’ were differentially enriched in the PI group. Enrichment plots from GSEA. GSEA results showing that (D) ‘Cytokine-cytokine receptor interaction’ and (E) ‘JAK/STAT signaling’ were differentially enriched in the PI group. Enrichment plots from GSEA. GSEA results showing that (F) ‘Aminoacyl tRNA biosynthesis’, (G) ‘Valine leucine and isoleucine degradation’ were differentially enriched in the RE group. Enrichment plots from GSEA. GSEA results showing that (H) ‘One carbon pool by folate’ and (I) ‘Oxidative phosphorylation’ were differentially enriched in the RE group. GSEA, gene set enrichment analysis; RE, recovery; PI, persistent injury.

Table VIII

Differentially expressed pathways enriched in the PI and RE mice according to gene set enrichment analysis.

A, PI group
Pathway name	NES	Nominal P-value	FDR value
Toll_Like_Receptor_Signaling	2.90	<0.001	<0.001
Nod_Like_Receptor_Signaling	2.81	<0.001	<0.001
Pathogenic_Escherichia_Coli_Infection	2.74	<0.001	<0.001
Cytokine_Cytokine_Receptor_Interaction	2.60	<0.001	<0.001
JAK_STAT_Signaling_Pathway	2.59	<0.001	<0.001
RIG_I_Like_Receptor_Signaling	2.53	<0.001	<0.001
Leishmania_Infection	2.51	<0.001	<0.001
Small_Cell_Lung_Cancer	2.45	<0.001	<0.001
Cytosolic_Dna_Sensing	2.42	<0.001	<0.001
Chronic_Myeloid_Leukemia	2.33	<0.001	<0.001
B, RE group
Pathway name	NES	Nominal P-value	FDR value
Aminoacyl_tRNA_Biosynthesis	-1.75	0.001	0.095
Valine_Leucine_And_Isoleucine_Degradation	-1.73	<0.001	0.065
One_Carbon_Pool_By_Folate	-1.70	0.006	0.065
Oxidative_Phosphorylation	-1.69	<0.001	0.053

PI, persistent injury; RE, recovery; NES, normalized enrichment score.

Discussion

Sepsis is a clinically common and refractory critical illness. As a serious complication of sepsis, SA-AKI is the most common type of AKI in the clinic, with a high mortality rate (8). At present, owing to the lack of specific and effective means of SA-AKI prevention and treatment, it is of critical importance to study the pathophysiological mechanism and investigate new therapeutic targets in SA-AKI. Therefore, RNA-Seq was used in the present study to investigate gene expression profiles related to mice with SA-AKI that exhibited persistent injury or recovery. The current study identified some prognosis-related genes, transcription factors, miRNAs and pathways associated with the differential gene expression profiles, including metabolic process, the transcription factor RELA and miR-30e-5p.

κ-means clustering was used to gain insight into the molecular pathways underlying different gene expression patterns. This divided the top 2,000 DEGs into four clusters. Genes in cluster B were upregulated in the RE group, and this cluster was strongly enriched in genes related to cellular metabolic processes, such as lipid metabolism and fatty acid metabolism. It has been reported that lipoproteins play an important role in preventing infection and inflammation (13). All lipoproteins can bind and neutralize toxic bacterial substances to modulate cytokine production during the inflammatory response, which can weaken the host response (14-16). Therefore, low lipoprotein levels may damage the innate immune response to endotoxins, leading to a deterioration of the systemic inflammatory cascade (15,17-20). Svahn et al (21) reported that dietary omega-3 fatty acids can increase sepsis survival and stimulate the immune system. Omega-3 fatty acids can be metabolized into a novel substance group termed resolvins (Rv) (22). It has been suggested that treatment with RvD1 and RvD2 can increase the survival of patients with sepsis (22-24). One putative mechanism underlying this process may be the combination of an increased frequency of neutrophils and precursor cells in bone marrow and restored phagocytotic capacity (21).

Genes in cluster C were highly upregulated in the PI group, and this cluster was associated with cytokines, cell death and apoptosis. Multiple pathways related to cell death are involved in sepsis due to the pro-inflammatory, anti-inflammatory, coagulation and complement systems. Apoptosis is a natural process of cell death, which is genetically programmed and plays a crucial role in normal physiology and in the pathophysiology of sepsis (25,26). Apoptosis can be initiated by pro-inflammatory cytokines, such as tumor necrosis factor α (TNF-α), interleukin (IL)-1 and IL-6(27). It has been suggested that prevention of cell apoptosis can improve survival in animal models of sepsis and endotoxemia (27).

With FDR <0.1 and FC >2 as the cutoff, 5,402 DEGs (2,256 downregulated and 3,146 upregulated genes) were identified, which have potential to be novel regulators and may play a role in the pathophysiological mechanism underlying SA-AKI development. To further understand the molecular pathways, enrichment analysis of the DEGs was performed. For the GO Biological Process analysis, consistent with the results of k-means clustering analysis, it was identified that the upregulated genes were predominantly involved in cellular metabolic processes. Furthermore, the downregulated genes were associated with cell response to stress, cytokine production, apoptosis and cell death.

As for the GO Cellular Component analysis, it was identified that upregulated DEGs in the RE group were enriched in the mitochondrion, cytoskeleton and centrosome. Mitochondria are cytoplasmic organelles with a double phospholipid membrane that generate energy via oxidative phosphorylation (28). Mitochondria are also associated with calcium homeostasis, intracellular reactive oxygen species (ROS) generation and cell signaling functions (29,30). In response to sepsis, the inflammatory cytokines of the innate immune response, including TNF-α, IL-1 and IL-6, can promote mitochondrial permeability transition, inhibit oxidative phosphorylation and improve ROS production (31,32). It has been demonstrated that mitochondrial dysfunction plays a critical role in SA-AKI (33). However, further studies are still needed to clarify the molecular mechanism underlying mitochondrial dysfunction in SA-AKI.

The cytoskeleton is a necessary dynamic structure for cells that plays an important role in regulating cell permeability under physiological and pathophysiological conditions (34). The cytoskeleton pathway is a known modulator of endothelial barrier function and microvascular permeability (35). It is well known that an inflammatory stimulus can lead to the development of capillary leak and tissue edema due to cytoskeletal rearrangement, which is a characteristic feature of sepsis; its development is one of the causes of organ dysfunction in sepsis (36). A meta-analysis of transcriptomic data demonstrated that the cytoskeletal pathway was upregulated in sepsis patients compared with controls (37). Previous studies have found that the breakdown in blood/endothelial barrier function plays a vital role in the pathogenesis and organ function and could be a new therapeutic target for sepsis (38,39).

Cytokine production is a milestone in the immune response during inflammation and is associated with mortality during sepsis (40). Vertii et al (41) found that centrosome integrity is critical for cytokine production. Impaired cytokine production may predispose the host to infection.

For the GO Molecular Function analysis, it was confirmed that the downregulated genes in the RE group were enriched in enzyme binding and receptor binding. Furthermore, the upregulated genes in the RE group were enriched in transporter activity and active transmembrane transporter activity. This could indicate that transporter activity and active transmembrane transporter activity may be involved in sepsis recovery. The Slc7a7 gene provides instructions for the production of a protein termed y+L amino acid transporter 1 (y+LAT-1), which is involved in transporting certain protein building blocks (amino acids), namely, lysine, arginine and ornithine. The transport of amino acids from the small intestines and kidneys to the rest of the body is necessary for the body to utilize proteins. The y+LAT-1 protein forms one part (the light subunit) of a complex called the heterodimeric cationic amino acid transporter. This subunit is responsible for binding to the amino acids that are transported (42). However, whether the Slc7a7 gene is involved in the recovery of sepsis requires further study.

The present study identified that the TF target genes RELA and STAT3 were significantly associated with the downregulated genes in the RE group. RELA is a REL-associated protein involved in NF-κB heterodimer formation, nuclear translocation and activation (43). As a nuclear transcription factor of various cells, NF-κB plays a crucial role in the coordination of innate and adaptive immune responses in sepsis (44). NF-κB binds with IκB-α and forms a p65/p60 dimer in the cytoplasm in physiological conditions (45). In stress, IκB-α isolates NF-κB and is transferred to the nucleus, where it binds proinflammatory cytokines and promotes the gene transcription of TNF-α and IL-6 to further activate NF-κB to expand the inflammatory response (46). STAT3 is a member of the STAT protein family and can mediate the expression of various genes that stimulate cells. In response to cytokines and growth factors, STAT3 can be phosphorylated by receptor-related JAK and translocate to the cell nucleus where it acts as a transcriptional activator (47). The STAT3 pathway plays an important role in inflammatory signaling cascades and its activation can play a crucial role in host-bacterial interactions (48). Cytokines, such as IL-6 and IL-10, can cause the phosphorylation of the tyrosine and serine residues of STAT3 via JAK. Activated STAT3 can in turn translocate into the nucleus where it binds to specific promotor sequences and regulate the transcription of target genes (49). Zhuo et al (24) demonstrated that RvD1 may improve survival and attenuate the degree of lung inflammation of septic mice by suppressing STAT3 and NF-κB expression through a mechanism that is partly dependent on SIRT1.

TF binding motifs (TFBMs) have conserved DNA sequence elements in their promoter regions. TFBMs can act as binding sites for transcription factors and coordinate the expression of genes in the promoter regions that they appear (50). As TFs typically bind to DNA at sites matching specific sequences motifs, knowledge of the motifs for a TF will be useful to determine a potential binding site of the TF (51). The present study detected that the DEG promoters were overrepresented with numerous G-rich motifs that are bound by AP-2 and other factors, such as bHLH and C2H2 ZF. AP-2 transcription factors (TFAP-2) constitute a family of closely related and evolutionarily conserved proteins that bind to the DNA consensus sequence GCCN3GGC and stimulate target gene transcription (52). AP-2 consists of five different proteins in human and mouse: AP-2α, AP-2β, AP-2γ, AP-2δ and AP-2ε (53). A previous study found that various inflammatory cytokines and prostaglandins can induce the expression of TFAP-2, which in turn causes aberrant activation of genes associated with hyperproliferation of mesangial cells and nephrosclerosis (54). Knocking down TFAP2A significantly decreases TF binding and the gene expression of collagen IV, which may play a critical role in diabetic nephropathy (55).

Our study demonstrated that miR-30e-5p, miR-181a-5p and miR-340-5p were significantly enriched target genes for the downregulated genes, and miR-466d-5p and miR-466l-5p were significantly enriched target genes for the upregulated DEGs. miRNAs and TFs share a common link as two vital gene regulatory molecules in multicellular organisms (56). miRNAs are a family of small non-coding RNAs that modulate gene expression in a sequence-specific manner, and as a result they regulate numerous cellular processes at the post-transcriptional level (57). miR-30e-5p is a member of the miR-30 family, and plays a critical role in renal development and maintaining renal function (58). Sun et al (59) reported that the circulating miR-30e-5p level is significantly elevated in contrast-induced AKI. This may be associated with the pathophysiology of AKI. The miR-181 family has been demonstrated to control inflammation under physiological and pathological conditions by modulating various key aspects of growth, development and activation (60). A previous study demonstrated that miR-340 can inhibit cell proliferation, tumor migration and invasion by targeting JAK via the JAK1/STAT3 signaling pathway (61). It has been indicated that miR-466l may to be a negative regulator of multiple proinflammatory cytokines, such as interferons and IL-10 (62,63). These miRNAs may be potential therapeutic targets for SA-AKI.

The innate immune system is the first line of defense against microbial invasion, relying on pattern recognition receptors to recognize external pathogenic microorganisms and then remove them (64). TLRs and NLRs are important receptors that mediate immune recognition by recognizing pathogen-associated molecular patterns derived from various microbes. These processes are bridges between innate immunity and adaptive immunity (65).

In conclusion, the present study examined gene expression profiles of a SA-AKI model that exhibited persistent renal injury or renal recovery using RNA-Seq. A number of prognosis-related genes, transcription factors, miRNAs and pathways were identified to be associated with the differential expression profile, including metabolic process, the transcription factor RELA and miR-30e-5p, which provide some basis for future experimental studies. Further validation, particularly in human tissues, may provide a more comprehensive understanding of the underlying pathophysiology during SA-AKI.

Table II

Enriched GO Biological Process terms for the upregulated and downregulated genes.

A, Downregulated genes
Adjusted P-value	Genes (n)	Biological process
<0.001	362	GO:0071310 cellular response to organic substance
<0.001	185	GO:0034097 response to cytokine
<0.001	358	GO:0002376 immune system process
<0.001	343	GO:0009605 response to external stimulus
<0.001	315	GO:0048584 positive regulation of response to stimulus
<0.001	321	GO:0008219 cell death
<0.001	155	GO:0071345 cellular response to cytokine stimulus
<0.001	243	GO:0006952 defense response
<0.001	141	GO:0001817 regulation of cytokine production
<0.001	298	GO:0006915 apoptotic process
<0.001	301	GO:0012501 programmed cell death
<0.001	180	GO:0051707 response to other organism
<0.001	149	GO:0001816 cytokine production
<0.001	214	GO:0080134 regulation of response to stress
B, Upregulated genes
Adjusted P-value	Genes (n)	Biological process
<0.001	293	GO:0044281 small molecule metabolic process
<0.001	151	GO:0006082 organic acid metabolic process
<0.001	138	GO:0019752 carboxylic acid metabolic process
<0.001	139	GO:0043436 oxoacid metabolic process
<0.001	186	GO:0044711 single-organism biosynthetic process
<0.001	99	GO:0032787 monocarboxylic acid metabolic process
<0.001	61	GO:0060271 cilium morphogenesis
<0.001	147	GO:0044255 cellular lipid metabolic process
<0.001	183	GO:0006629 lipid metabolic process
<0.001	72	GO:0006631 fatty acid metabolic process
<0.001	40	GO:0044782 cilium organization
<0.001	150	GO:0019637 organophosphate metabolic process
<0.001	163	GO:0055114 oxidation-reduction process
<0.001	76	GO:0006820 anion transport
<0.001	54	GO:0044282 small molecule catabolic process

GO, Gene Ontology.

Table V

Enriched transcription factor binding motifs for the downregulated genes.

Adjusted P-value	Genes (n)	TF binding motifs
<0.001	14	TRANSFAC: RELA
<0.001	12	TRANSFAC: STAT3
<0.001	7	TRANSFAC: REL
<0.001	14	TRANSFAC: CEBPB
<0.001	12	TRANSFAC: CREB1
<0.001	10	TRANSFAC: JUN
<0.001	8	TRANSFAC: FOS
<0.001	11	TRANSFAC: SFPI1
<0.001	6	TRANSFAC: CEBPD

0.002 3 TRANSFAC: NFKB1

Table VII

Enriched miRNA target genes of the upregulated and downregulated genes.

A, Downregulated genes
Adjusted P-value	Genes (n)	Target miRNAs
<0.001	71	MiRTarBase:mmu-miR-30e-5p
<0.001	80	MiRTarBase:mmu-miR-181a-5p
<0.001	93	MiRTarBase:mmu-miR-340-5p
<0.001	73	MiRTarBase:mmu-miR-124-3p
<0.001	21	MiRTarBase:mmu-miR-155-5p
<0.001	14	MiRTarBase:mmu-miR-338-5p
<0.001	57	MiRTarBase:mmu-miR-129-5p
<0.001	81	MiRTarBase:mmu-miR-329-3p
<0.001	62	MiRTarBase:mmu-miR-466i-3p
<0.001	53	MiRTarBase:mmu-miR-19b-3p
<0.001	33	MiRTarBase:mmu-miR-377-3p
<0.001	25	MiRTarBase:mmu-miR-342-3p
0.003	71	MiRTarBase:mmu-miR-15a-5p
0.004	52	MiRTarBase:mmu-let-7b-5p
B, Upregulated genes
Adjusted P-value	Genes (n)	Target miRNAs
<0.001	84	MiRTarBase:mmu-miR-466d-5p
<0.001	84	MiRTarBase:mmu-miR-466l-5p
<0.001	30	MiRTarBase:mmu-miR-150-5p
<0.001	67	MiRTarBase:mmu-miR-1195
<0.001	27	MiRTarBase:mmu-miR-5127
<0.001	69	MiRTarBase:mmu-miR-3089-5p
<0.001	29	MiRTarBase:mmu-miR-466n-5p
<0.001	51	MiRTarBase:mmu-miR-466m-3p
<0.001	41	MiRTarBase:mmu-miR-3064-5p
0.003	45	MiRTarBase:mmu-miR-665-3p
0.003	61	MiRTarBase:mmu-miR-129-5p