Identifying of miRNA–mRNA Regulatory Networks Associated with Acute Kidney Injury by Weighted Gene Co-Expression Network Analysis

Introduction

Acute kidney injury (AKI) is a clinical syndrome of rapidly declining kidney function caused by multiple etiologies and manifested by rapid increase in serum creatinine, decrease in urine output, or both.¹ The incidence of AKI in Hospitalized adult patients ranged from 3.0% to 18.3% globally.² These differences may result from factors such as different backgrounds of patients and the different criteria used to define AKI. However, even the low end of the AKI incidence range makes it abundantly clear that AKI affects a very large number of people worldwide. AKI causes approximately 1.7 million deaths worldwide each year, and this approximation may even be underestimated in view of the silent nature of AKI.^3,4

AKI has become a worldwide health concern, but its diagnosis and treatment remain challenging. The current diagnosis of AKI still relies on serum creatinine and urine output.⁵ Unfortunately, owing to the lag time between the onset of injury and changes in serum creatinine and urine output, the diagnosis is often not clarified until after the initiation and extension phases of kidney injury have transitioned to the maintenance phase of injury.⁶ The available evidence suggests that even mild AKI is associated with increased mortality.^3,4 No interventions to improve established AKI outcomes have yet been developed, and its therapeutic strategy is limited to supportive care, in which renal replacement therapy (RRT) plays an important role.⁷ However, clinicians also lack effective decision-making tools on whether and when RRT should be commenced to improve the outcome of these patients with AKI, considering that AKI is usually a continuum of kidney injury rather than a single-hit, freestanding condition.^7,8

Early diagnosis and timely intervention can minimize kidney injury and promote kidney function recovery.⁹ Although several novel AKI biomarkers, including kidney injury molecule 1 (KIM-1), IL-18, cystatin C, and neutrophil gelatinase-associated lipocalin (NGAL) have been intensively studied in the last decades, there is still a lack of recognized biomarkers for clinical applications.¹⁰ Therefore, in this study, we searched for potential biomarkers of AKI with the aim of guiding clinical interventions and providing valuable insights into the complex molecular mechanisms of AKI.

Materials and Methods

Microarray Data Processing

The GSE85957 microarray dataset was concluded and downloaded from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) through the GEO query package in the R environment (Version 4.2.0). The GSE85957 dataset (platform: GPL1355; Rat230_2; Affymetrix Rat Genome 230 2.0 Array) compared gene expression profiles of kidney tissues from male Han Wistar rats (8 weeks old) with intraperitoneal injection of cisplatin (0, 1 or 3 mg/kg) for 3, 5, 8 and 26 days induced AKI model (n = 38) and controls (n = 19). The corresponding annotation file-GPL1355 matrix containing more than 31,000 probe sets and more than 28,000 annotated rat genes were obtained for conversion of probes to target genes. If two or more probes were annotated in a gene, the mean was calculated for further analysis. The matrix was preprocessed and normalized by using the bead array software package. Outlier microarray samples were identified using a sample network approach based on Euclidean distance and a cutoff value of −5 for Z.ku was calculated as ku-mean(k)/sqrt(var[k]).¹¹

Weighted Gene Co-Expression Network Analysis (WGCNA) Construction

WGCNA was performed to identify clusters with high correlation using the WGCNA package.¹¹ The soft threshold β was evaluated by the R function pickSoft Threshold for scale-free topology. A β value of 5 (R² > 0.8) was chosen to construct the gene network by applying the default WGCNA method. The topology matrix was obtained based on the β value, and the topology matrix was clustered and analyzed, and then the tree was cut into different modules using the dynamic cut method, and the modules with correlation higher than 0.8 were merged, and finally a total of 28 modules were obtained.

Identification of Significant Modules

For detection of module which were significantly associated with AKI and dose, we performed Pearson correlation analysis to assess the associations. Modules with P-value <0.01 were considered to be significantly associated with AKI. In the co-expression network, the module eigengene (ME) as the summary of the gene expression profiles was correlated with these characteristics to find the most important associations. Gene significance (GS) is used to determine the connective degree between the gene and the specific trait in the association quantification of individual genes with AKI, which was further measured by quantification of module membership (MM) in the correlation of ME and the gene expression profile.

Conduction of AKI Models and mRNA Sequencing

In this study, mRNA sequencing was performed by intraperitoneal injection of cisplatin to induce AKI model to detect differentially expressed genes. All SD rats (male, 6 weeks, n = 6) were purchased from Shanghai Slac Laboratory Animal Co. Ltd and housed in the Experimental Animal Center of Tongji University. Rats were raised in a 2-hour dark-light cycle with free access to feed and water. And the room temperature was maintained in 20–24 °C with relative humidity of 45 ~ 55%. Cisplatin at 3 mg/kg was diluted in 0.9% saline and then administered intraperitoneally. One week later rat kidneys were collected for further mRNA sequencing. The total RNA was extracted using the TRIzol (Invitrogen, USA) and quantified by NanoDrop 2000 spectrophotometer (Thermo, USA) and the high-quality RNA sample (OD260/280 = 1.8~2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, 28S:18S ≥ 1.0, > 1μg) was further applied in the sequencing. RNA-seq transcriptome library was prepared using 1μg RNA by TruSeq^TM RNA sample preparation Kit from Illumina (San Diego, CA) and the expression of mRNA was profiled using HiSeq xten/NovaSeq 6000 sequencer (Illumina, Inc.).

Identification of Critical Genes

The identification of critical genes was conducted based on significant differential expression between AKI and control samples. The limma package was separately applied to identify differential expressed genes (DEGs) in the two traits with the cut-off criteria |log2 fold change (FC)| ≥ 1 and P value < 0.05. And the critical genes were subsequently visualized among the DEGs and MEs through the Venn diagrams package. Additionally, the diagnostic performance of critical genes was also verified by receiver operating characteristic (ROC) curves in this study, and the results were presented as the area under the curve (AUC).

Functional Enrichment Analysis of Critical Genes

Gene Ontology (GO) and KEGG enrichment analysis were applied in the exploration of the related biological function in critical genes through Metascape (http://metascape.org/). The visualization of top 10 GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms were established by using ggplot2 package.

Construction of miRNA–mRNA Regulatory Network

TargetScan was used to predict miRNAs interacting with eight critical genes, and the relevant miRNAs were screened using score.percentile > 95 as criteria, and the miRNA–mRNA network was visualized in Cytoscape 3.5.1. The workflow of this study is shown in Figure 1.

Figure 1 The workflow of this study.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

Total RNA was extracted from tissue specimens using Trizol reagent (TaKaRa, China). cDNA was obtained by reverse transcription and qRT-PCR reactions were performed according to the previous study steps.¹² The relative expression of each gene was calculated using the 2^−ΔΔCt method. All primer sequences are shown in Table S1.

Results

Data Processing

The GSE85957 matrix data file was downloaded and 15323 gene expression data were annotated using the GPL1355 platform annotation file. The dataset contained 38 AKI samples and 19 control samples, and all samples had Z.ku higher than the standard Z.ku (−5) in the Euclidean distance-based sample network analysis. Therefore, expression data from all 57 samples were applied in the WGCNA.

WGCNA Construction

In WGCNA, the R function pickSoft Threshold was used to evaluate the soft threshold β. We chose a β value of 5 (R²>0.8) and constructed the gene network by applying the default WGCNA method (Figure 2A). Based on the β values to the topology matrix and further analysis, a total of 28 modules were finally obtained (Figure 2B). We then evaluated the association between modules eigengenes and AKI and dose using Pearson correlation analysis. Among the modules, black (AKI: P = 2e-08; Dose: P = 1e-09), royalblue (AKI: P = 0.007; Dose: P = 4e-05), steelblue (AKI: P = 0.004; Dose: P = 3e-05), darkgrey (AKI: P = 8e-05; Dose: P = 1e-10), blue (AKI: P = 9e-04; Dose: P = 5e-07), midnightblue (AKI: P = 0.003; Dose: P = 2e-06), orange (AKI: P = 0.004; Dose: P = 1e-11), salmon (AKI: P = 0.005; Dose: P = 3e-04), and paleturquoise (AKI: P = 4e-04; Dose: P = 3e-06) were strongly correlated with AKI and dose (P < 0.01) (Figure 2C). Correlation line analysis of the modules with AKI and dose revealed that the MEorange module (103 genes) had the highest positive correlation with AKI (cor = 0.97, p = 1e-73) and dose (cor = 0.98, p = 6.9e-84), and the MEblack module (30 genes) had the highest positive correlation of genes with AKI (cor = −0.98, p = 1e-24) and dose (cor = −0.97, p = 7.6e-22) had the highest negative correlations (Figure 2D). In addition, gene expression in the orange module was positively correlated with AKI and dose, while the opposite result was observed for the black module (Figure 2E).

Figure 2 WGCNA selected two modules associated with AKI and dose. (A) Analysis of scale-free fit index and mean connectivity for different soft thresholds (β). (B) Gene enrichment modules in the co-expression network with 28 modules identified. (C) Correlation between each module and AKI and dose. (D) Pearson correlation analysis of Orange and black modules for gene significance (GS) and module membership (MM). (E) Relationship between gene expression in Orange and black modules with AKI and dose. (*p<0.05; ***p<0.001; ****p<0.0001).

Module Genes and Functional Enrichment Analysis

We uploaded 103 genes from the MEorange module and 30 genes from the MEblack module to Metascape to explore gene related biological functions. Figure S1 shows the top 10 biological processes (BPs) and KEGG pathways annotated with genes in the orange and black modules. The BP annotations show that genes in the orange and black modules are significantly enriched in skeletal muscle cell differentiation, skeletal muscle tissue development, skeletal muscle organ development, muscle structure development, striated muscle tissue development, muscle tissue development, muscle organ development, rhythmic process, regulation of neurogenesis and fat cell differentiation (Table 1). The KEGG pathway enrichment analysis shows that genes in the orange and black modules are mainly involved in circadian rhythm, bladder cancer, MAPK signaling pathway, p53 signaling pathway, bell cycle, HTLV-I infection, apelin signaling pathway, hepatitis B, IL-17 signaling pathway and endocrine resistance (Table 1).

Table 1 Top 10 GO Biological Processes and KEGG Pathway for Orange and Black Modules

Identification and Validation of Critical Genes

To screen for critical genes, we analyzed gene expression between AKI and controls in the GSE85957 dataset and our own experimentally processed data (raw count) using the limma software package. The screening criteria were |log2 fold change (FC)| ≥ 1 and P value < 0.05. The DEGs in the GSE85957 dataset and raw count are shown in Figure 3A. In the GSE85957 dataset, 74 differentially expressed genes (65 up-regulated and 9 down-regulated) were obtained and raw count yielded 2596 differentially expressed genes (603 up-regulated and 1993 down-regulated). We intersected the GSE85957 dataset, raw count, orange module and black module, and the results showed that eight genes (Atf3, Egr1, Egr2, Fos, Fosb, Gdf15, Serpine1, Nr1d1) were screened as DEGs (Figure 3B). Table 2 shows the BPs annotated with the eight critical genes. These eight up-regulated genes were defined as important in the development of AKI (Figure 4A and Figure S2). In addition, the ROC curves showed that all critical genes had AUCs above 0.8 (Figure 4B), indicating that they could be used as potential biomarkers for the diagnosis of AKI. Moreover, we also predicted miRNAs interacting with eight genes in TargetScan and constructed a possible miRNA–mRNA regulatory network (Figure 5).

Table 2 GO Biological Processes for Eight Critical Genes

Figure 3 Selection of critical genes. (A) Volcano map of differential expressed genes (DEGs) in GSE85957 and our raw data. (B) Screening of eight critical genes by Venn diagram.

Figure 4 Expression profile of critical genes. (A) Expression profile of critical genes in GSE85957 dataset. (B) ROC curve of critical genes.

Figure 5 miRNA–mRNA regulatory network of critical genes.

In addition, we extracted total RNA from the kidneys of 15 AKI and 15 control rats, and verified the expression of the above eight critical genes by qRT-PCR, and found that Egr1 and Fos expression was increased in AKI (Figure 6), demonstrating that Egr1 and Fos may be potential key biomarkers of AKI.

Figure 6 Relative expression of eight critical genes in tissue samples.

Discussion

AKI is a clinical syndrome caused by multiple clinical factors and characterized by a dramatic decline in renal function. AKI patients are often accompanied by metabolite retention and water, electrolyte, and acid-based imbalances.² Unfortunately, without early treatment, AKI patients are associated with a very poor prognosis and may progress to chronic kidney disease (CKD) and end-stage renal disease, resulting in increased mortality.¹³ Therefore, exploring the genes and modules specific to AKI can help in early diagnosis and treatment, thus effectively reducing mortality and adverse prognostic responses in patients.

In the present study, according to WGCNA, we found two specific modules highly associated with AKI, containing a total of 133 hub genes. Development-related processes in skeletal muscle and inflammatory response-related processes were significantly involved in the pathology of AKI as revealed by BP and KEGG enrichment analysis. Eight of these hub genes were differentially expressed between the AKI and control groups and were identified as playing a key role in the development of AKI, with promising applications in early predicting AKI. Finally, we predicted and integrated the miRNAs interacting with eight hub genes in TargetScan.

Within these critical genes, Egr1, Egr2, Fos and Fosb are closely associated with cytoprotective or immune-related responses. Egr1 and Egr2, as members of the zinc finger transcription factor family, are characterized by a highly conserved DNA binding domain consisting of three zinc finger motifs.¹⁴ Egr1 gene is an important transcription factor, which is located in the human chromosome 5q23-31 region and is mainly involved in processes such as tissue damage, immune response and fibrosis.¹⁵ Egr1 can promote T cell activation after T cell receptor (TCR) stimulation.¹⁶ EGR2 is induced by TCR binding and plays a key role in the induction of T cell incompetence.¹⁷ After T cell activation, EGR2 is highly induced in order to control the inflammatory response and thus negatively regulate T cell activation.¹⁸ TCR signaling upregulates Fas ligand (FasL) on CD4⁺ T cells, while EGR2 acts as a potent activator of FASL expression to maintain immune homeostasis by inducing apoptosis.¹⁹ In the present study, Egf1 and Egf2 were significantly upregulated in the AKI group compared to the control group, suggesting that elevated Egf1 and Egf2 transcripts in AKI patients are associated with T lymphocytes activation and proliferation. This is consistent with previous studies finding that after AKI, T cells are dramatically elevated in renal draining lymph nodes, transferring and producing pro-inflammatory cytokines to the injured kidney, activating macrophages and promoting renal inflammation.²⁰ Thus, Egf1 and Egf2 may be a predictive biomarker for patients with AKI. Fos and Fosb are members of the activator protein-1 family (AP-1) of transcription factors. Fos and Fosb are regulators of cell proliferation, differentiation and transformation. They are characterized by a conserved basic leucine zipper domain, and a trans-activating domain that is primarily responsible for providing binding sites for other transcriptional co-regulators.²¹ However, it is unknown the role of Fos and Fosb in AKI and further research to explore the underlying mechanisms is needed.

Activating transcription factor 3 (ATF3) as a member of the ATF/cAMP response element-binding (CREB) family binded to the cyclic AMP response element (CRE) in numerous promoters.²² ATF3 is a stress-induced transcription factor, which has key roles in modulating metabolism, immunity, and oncogenesis. A variety of extracellular signals, such as endoplasmic reticulum (ER) stress, cytokines, chemokines, and LPS, have been implicated in the induction of ATF3.²³ It had been reported that ATF3 promoted Gadd45β/γ transcription and glomerular mesangial cells apoptosis induced by complement C5b-9.²⁴ Several studies verified in both mouse models and patients that urinary exosomal ATF3 (uATF3) expression was higher in the AKI group than in the normal group.^25,26 Therefore, uATF3 could a biomarker for early AKI. Growth differentiation factor 15 (GDF15), a member of the transforming growth factor beta family, was associated with immunosuppression, anti-apoptosis and anti-inflammation.²⁷ Furthermore, several studies had shown that GDF15 expression was upregulated in AKI induced by different factors.^28,29 However, the effect of GDF15 on AKI was still unclear.

Serpine1 encoding plasminogen activator inhibitor-1 (PAI-1), would suppress fibrinolysis and enhance renal interstitial fibrosis which were precisely the pathological characteristic of AKI.^30,31 Moreover, the expression of PAI-1 was significantly correlated with renal inflammatory injury,^32,33 and knocking down PAI-1 would play a protective role.³⁴ Another study demonstrated that the Serpine1 gene expression was notably upregulated in cisplatin-induced AKI kidney tissues.³⁵ Thus, Serpine1 may be a potential marker for AKI. NR1D1, as the member of nuclear receptor subfamily 1 group D, functioned as transcriptional repressors with haem as their endogenous ligand.³⁶ NR1D1 was considered as the components of the circadian clock system,³⁷ which play an essential role in regulating circadian rhythms.³⁸ NR1D1 was also involved in the regulation of many physiological processes and diseases, including inflammatory responses and cancer.³⁹ Researchers found inhibition of NR1D1 attenuated ferroptosis-mediated AKI.⁴⁰ The mechanism of NR1D1 in AKI was still unclear and needed further study.

With the advancement of sequencing technology, genomics is playing an essential role in disease diagnosis, mechanism research and treatment. Now, sequencing technology is more and more used in clinic, and key genes may play an important role in the occurrence and development of some disease. Thus, the diagnosis of AKI could be determined by the expression of key genes and the AKI can be treated by intervention of key genes, so as to improve patient prognosis and reduce the occurrence of serious complications.

Conclusion

In conclusion, by analysis the results of sequencing and validation, we identified eight critical genes in AKI, including GDF15, Serpine1, ATF3, NR1D1, Egr1, Egr2, Fos, Fosb, and verified the concordance of Egr1 and Fos in tissue samples. Their abnormal expression may activate multiple signaling pathways, ultimately leading to the development and progression of AKI. Furthermore, we constructed a related miRNA–mRNA regulatory network based on the above eight critical genes, which may provide new insights into the molecular pathogenesis of AKI.