The analysis of a time-course transcriptome profile by systems biology approaches reveals key molecular processes in acute kidney injury.

BACKGROUND: Acute kidney injury is a common debilitating disease with no curative treatment. The recent development of big biological data is expected to expand our understanding of the disorder if appropriately analyzed to generate translational knowledge. We have here re-analyzed a time-course microarray data on mRNA expression of rat kidneys exposed to ischemia-reperfusion to identify key underlying biological processes.
MATERIALS AND METHODS: The dataset was quality controlled by principal component analysis and hierarchical clustering. Using limma R package, differentially expressed (DE) genes were detected which were then clustered according to their expression trajectories. The biological processes related to each cluster were harvested using gene ontology enrichment analysis. In addition, the interaction map of proteins encoded by the DE genes was constructed, and the functions related to network central genes were determined. Furthermore, signaling pathways related to the DE genes were harvested using pathway enrichment analysis.
RESULTS: We found 8139 DE genes that drive critical processes such as the control of blood circulation, reactive species metabolism, mitochondrial respiration, apoptosis, cell proliferation, as well as inflammatory and immunological reactions. The role of less recognized pathways such as olfactory signaling in acute kidney injury is also proposed that remains to be investigated in future studies.
CONCLUSION: Using systems biology top-down approach, we have suggested novel potential genes and pathways to be intervened toward kidney regeneration.

INTRODUCTION

Despite huge investigations on acute kidney injury, it is yet a common cause of death in hospital-admitted patients[1234] and the existing therapeutic options are just supportive.[5] The insufficiency of current therapies for this clinical condition is partly due to our limited understanding of the molecular mechanisms. The recent availability of omics data in line with the emerging hypothesis-free paradigm of systems biology paves the way for more comprehensive insights toward the molecular pathophysiology of the disorder.[6] This holistic approach is expected to accelerate the development of novel efficient therapeutics.

Although a huge amount of biological big data has been generated in the last few years, time-course expression profiling studies encompass a minority of the datasets. Considering the dynamic nature of biological processes, time-series evaluations are more appreciated as they represent the course of events more precisely. Here, we selected a dataset deposited by Speir et al. in which the mRNA expression profile of rat kidneys exposed to ischemia-reperfusion (IR) injury is assessed in several time steps following the insult.[7] Following quality control steps, we have re-analyze this dataset to identify key role player genes as well as critical signaling pathways and molecular processes related to acute kidney injury.

MATERIALS AND METHODS

Microarray data

mRNA expression profiling deposited by Speir et al.[7] with the accession number GSE58438 was retrieved from the Gene Expression Omnibus database.[89] The series matrix file of the data has been normalized by the method of Sketch-Quantile. Dataset quality was assessed by principal component analysis using ggplot2 package[10] and hierarchical clustering using pheatmap[11] package of R software[12] which are accessible through CRAN repository. The clustering method was “complete” with “Euclidean” distance measure. Differentially expressed (DE) genes were identified using limma package[13] of R software, available in Bioconductor repository. In the analysis, IR rats harvested at 3, 24, and 120 h postischemia as well as normal controls were compared by the method of F-statistic. Benjamini–Hochberg false discovery rate (FDR) was then used for P value correction. The genes with adjusted (adj.) P < 0.05 were considered as DE. Venn diagram of DE genes was depicted using Eulerr[14] package of CRAN repository.

Gene clustering

For clustering of the DE genes, Short Time-Series Expression Miner version 1.3.11[6] was used. For this method, the matrix files of DE genes containing the expression values of three samples in each group were provided. FDR was used as the correction method in model profile section.

Gene ontology enrichment analysis1

The ClueGO plugin[15] version 2.3.5 of Cytoscape[16] version 3.2.1 was applied to retrieve the gene ontology (GO) data from GO database.[17] Bonferroni step-down method was chosen for P value correction and adj. P < 0.05 was considered significant. REViGO[18] was used to find the parent terms of the retrieved GO terms and the similarity index was set at 0.4.

Protein-protein interaction network

Cytoscape CluePedia plugin[19] version 1.3.5 was used for the construction of protein-protein interaction (PPI) network. The activation, inhibition, binding, and posttranslational modification interactions were retrieved from STRING database,[20] and the confidence cutoff was set at 0.8. The network topology was analyzed using network analyzer tool.

Pathway enrichment analysis

Cytoscape ClueGO plugin was employed for pathway enrichment analysis. KEGG[2122] and REACTOME[23] databases were used to find related pathways. Bonferroni step-down method was applied for P value correction, and pathways with adj. P < 0.05 were selected. Ram and Wickham[24] package was applied for drawing the scatter plot.

RESULTS

In this study, we have re-analyzed the microarray dataset “GSE58438” which includes the mRNA profiles of rat kidneys at 3, 24, and 120 h following IR injury as well as normal controls. To assess the quality of the data, principal component analysis and hierarchical clustering were employed. Both these unsupervised methods revealed the segregation of the samples in four subsets matching with the time steps of the study [Figure 1a]. However, three samples that were not scattered according to the experiment group were omitted in downstream analyses to improve the quality of the dataset.

Figure 1

The quality of the microarray dataset was satisfactory as most samples were segregated according to the experimental groups in unsupervised principal component analysis and hierarchical clustering. Three samples that did not follow this trend were omitted in downstream analyses (a). The numbers of DE genes in each IR time point are depicted in the Venn diagram (b)

To identify DE genes, limma package of R software was employed which revealed 8139 genes to be variably expressed with FDR <0.05 [Supplementary Table 1; Figure 1b]. As the animal model used in this study is generated with a transient episode of ischemia, we were interested to focus on the profiles of DE genes related to this condition and plotted their expression trajectories over time [Figure 2]; Actb, the activator of Siah2, increases following IR which is associated with a transient Siah2 overexpression. In addition, Hif1an (Fih1) declines and Ard1a increases shortly after injury. These alterations can lead to Egln2 (Phd1) inactivation and Hif1a stabilization. Therefore, the expression of these genes is according to their role for protecting the cells against hypoxia. In addition, Cxcl12 (Sdf1) shows a declining trend which is in line with our previous data indicating the downregulation of this gene both at mRNA and protein levels following kidney IR.[25] Furthermore, Tgfb is observed to be overexpressed which is in agreement with previous studies [Figure 2].[2627]

Figure 2

The temporal expression profiles of some differentially expressed genes related to tissue hypoxia are shown

Next, to provide a holistic view, the DE genes were clustered based on their temporal expression patterns, and the functions attributable to each group were identified by GO term enrichment analysis. Interestingly, in each cluster, one or two processes were highlighted, indicating that the genes with similar trajectories drive similar functions. The identified processes were in accordance with the molecular pathophysiology of the disorder and included reactive species metabolism, lipid peroxidation, mitochondrial metabolism, cell proliferation, regeneration, cilium organization, and immunological reactions [Figure 3].

Figure 3

The DE genes were clustered based on their expression trajectories and gene ontology enrichment analysis was performed for each cluster. Parent gene ontology terms are demonstrated

To assess the functional relationships between the proteins encoded by the DE genes, a PPI map was constructed, and its topology was analyzed. As expected for biological networks,[28] the node degree in the constructed network followed a power low distribution with the form y = 1548.5x− 1.516, where y is the number of nodes and x is degree (correlation: 0.965). We have previously shown that central nodes in PPI are critical for cell functions.[2930] Therefore, based on degree, betweenness, and closeness centrality parameters, top 1% DE genes were selected [Figure 4a] and their ontologies were assessed [Figure 4b]. These central genes were mainly related to regulation of blood circulation, angiotensin signaling, reactive oxygen species metabolic processes, mitochondria organization, apoptosis, and cell proliferation. These processes are known to be of special importance in the pathogenesis of acute kidney injury.[313233]

Figure 4

Based on the degree, betweenness centrality, and closeness centrality, top 1% genes in the PPI network were determined (a). Gene ontology enrichment analysis identified key biological processes, cellular components, and molecular functions related to these central genes (b)

To further provide clues regarding to the biological functions in kidney ischemia, pathway enrichment analysis was performed with the DE genes and the signaling pathways with FDR <0.05 were harvested [Figure 5]. In agreement with previous studies,[3134353637] this analysis highlights the importance of processes such as mitochondrial respiration, apoptosis, cell cycle, DNA repair, and nuclear factor-kappa B activation in acute kidney injury. In addition, some pathways such as olfactory signaling are enriched that can be proposed as novel role players in this disorder. Although olfactory receptor family encompasses a large amount of genes, the employment of P value correction in the enrichment analyses precludes the chance that the enrichment of this pathway is a false-positive record.

Figure 5

Pathway enrichment analysis was performed for all differentially expressed genes and signaling pathways with false discovery rate <0.05 are demonstrated

DISCUSSION

The 20th century was a time for great achievements in medicine such as the control of communicable diseases and the development of surgical procedures and diagnostic techniques. However, the treatment options for most complex disorders are yet supportive rather than curative. The recent emergence of systems biology has raised hopes to address this limitation through more comprehensive investigations that may lead to novel therapeutics. Although a huge amount of big data has been generated in the last years, they are not commonly appropriately analyzed to generate translational knowledge. Indeed, the speed of data generation has surpassed that of data analysis and so it has been advised to focus on the interpretation of available datasets rather than the generation of new ones.[3839] In this study, we have re-analyzed a time-course microarray mRNA profiling of acute kidney injury and generated a comprehensive outlook of molecular interactions, key processes, and signaling pathways to address the scarcity of holistic investigations of this disorder.

Although the current dataset is originally developed as a part of a study on the effect of valproic acid on kidney ischemia,[7] we were here interested on the natural course of acute kidney ischemia and focused on the samples derived from rats exposed to IR but not receiving the drug. We have shown that the selected samples are acceptable for analysis according to unsupervised quality control measures. Another advantage is that the samples are harvested in different time points following the insult and so can represent the dynamic nature of the biological events. Indeed, single time-point expression profiling can be misleading due to complex trajectories and noise in genes expression.[40] However, time-course microarray experiments are not common due to practical limitations. In addition, the analysis of such data is challenging and requires specific tools.[41] The use of limma in the current work is based on a previous study indicating that this R package has an acceptable performance for the identification of DE genes in time-course microarray datasets and outperforms other similar algorithms.[42]

Following the detection of DE genes by limma, we focused on the expression patterns of special genes related to tissue hypoxia among which Cxcl12 is interesting. This chemokine is positively regulated by Hif1 and is responsible to attract stem cells and leukocytes to injured tissues.[4344] Therefore, it is expected to be upregulated in response to ischemia,[45] but it showed a declining trend in this dataset. We previously had a similar observation on Cxcl12 downregulation following kidney ischemia and using a stochastic Petri net mathematical model, proposing that this unexpected observation is due to the inhibitory effect of histone deacetylase on this gene which dominates Hif1-induced overexpression.[25] This prediction is recently experimentally corroborated (unpublished data).

To explore the underlying mechanisms of this disorder, the DE genes were clustered based on their temporal expression patterns and the functions attributable to each cluster were detected. Furthermore, the interaction map of proteins encoded by the DE genes was constructed, and the ontology of the central network genes was assessed. Furthermore, enrichment analysis was performed to identify the signaling pathways in which the DE genes are involved. These investigations highlight the importance of essential processes such as free radical metabolism, cellular respiration, immunological responses, cell proliferation, and apoptosis. In addition, the involvement of olfactory receptors in acute kidney injury is here proposed which to the best of our knowledge has not been previously indicated. Interestingly, olfactory receptors are shown to be involved in normal renal function,[4647] as well as polycystic kidney disease[48] and nondiabetic chronic kidney disease.[49] Future experimental studies are required to investigate if these receptors also contribute in the pathogenesis of acute kidney injury.

Financial support and sponsorship

This study was supported by Isfahan University of Medical Sciences (grant number 395115) and Iran National Science Foundation (grant number 9600608).

Conflicts of interest

There are no conflicts of interest.