Peritubular capillary rarefaction and hypoxia in the RIF
Results of Sirius red staining shown in Fig. 2A, obviously revealed a higher ratio of collagen deposition in the RIF specimens compared with the normal group. And when combined with the micro-vessel density, we found that the positive CD 34 staining definitely depressed in the RIF patients than the normal and was reversely related to fibrosis degree in RIF, indicating remarkable PTC loss in RIF (Fig. 2B). Recent studies have shown that PTC rarefaction deteriorate microenvironment hypoxia, both of which contributed to kidney fibrosis. That’s the reason why hypoxia HK-2 cell model was established and validated with RT-qPCR for further investigation in this work (Fig. 2C and 2D).
Identification Of Degs
As can been from scatterplots (Fig. 3), high correlation of samples confirmed the superior reproducibility and reliability of our RNA sequencing. A total of 572 DEGs were identified from the HK-2 hypoxia vs HK-2 control profile datasets, indicating numerous biological alterations in RIF. The most changed genes and their hierarchical clustering results were plotted in Figs. 4A and 4B. The comparisons of the DEGs at cellular level suggested that notable alterations in gene expression under hypoxia, which were partly in favor of the hypothesis that multiple genes contributed to the differential angiogenesis state among RIF in vivo.
Go Functional And Pathway Enrichment Of Degs
To further analyze the specific biological functions and enriched signaling pathways of DEGs, we tried to analyze GO and KEGG enrichments of these DEGs, which are demonstrated in Fig. 4C. From annotations of Top 5 cellular component (CC) in GO analysis, the cell, cell junction, organelle, organelle part and membranes were all included. Molecular function (MF) annotations demonstrated that binding and catalytic activity played the most critical roles. We focused on biological processes (BP) particularly, which demonstrated that cellular process, single-organism process, biological regulation and metabolic processes were enriched. Worthy to note, for angiogenesis analysis of biological process, 15 upregulated and 3 downregulated genes were identified respectively to modulate angiogenesis. Importantly, six of them were further selected as potent targets for their obvious change after hypoxia treatment, including Coagulation factor III, tissue factor (F3), Angiopoietin like 4 (ANGPTL4), Adrenomedullin (ADM), Serpin family E member 1 (SERPINE1), BTG anti-proliferation factor 1 (BTG1) and Endothelial PAS domain protein 1 (EPAS1) (Details were revealed in Table 1).
Table 1
Analysis of potential potent genes in response to hypoxia
Gene | Description | Log 2 fold change | p-adjust | Regulated |
VEGFA | Vascular endothelial growth factor A | 1.666 | 5.682e-292 | Upregulated |
SLC2A1 | Solute carrier family 2, member 1 | 2.042 | 0 | Upregulated |
P4HA1 | Prolyl 4-hydroxylase subunit alpha 1 | 1.286 | 1.109e-140 | Upregulated |
PGK1 | Phosphoglycerate kinase 1 | 1.610 | 1.211e-259 | Upregulated |
BNIP3 | BCL2 interacting protein 3 | 2.241 | 5.147e-301 | Upregulated |
PDK1 | Pyruvate dehydrogenase kinase, isozyme 1 | 1.433 | 1.084e-112 | Upregulated |
SLC2A3 | Solute carrier family 2, member 3 | 1.760 | 5.497e-199 | Upregulated |
GYS1 | Glycogen synthase 1 | 1.386 | 8.639e-182 | Upregulated |
ALDOA | Aldolase, fructose-bisphosphate A | 1.011 | 7.025e-194 | Upregulated |
HK2 | Hexokinase 2 | 2.161 | 7.579e-295 | Upregulated |
F3 | Coagulation factor III, tissue factor | 2.300 | 0 | Upregulated |
ANGPTL4 | Angiopoietin-like 4 | 2.766 | 8.531e-304 | Upregulated |
ADM | Adrenomedullin | 2.335 | 5.691e-239 | Upregulated |
SERPINE1 | Serpin family E member 1 | 1.339 | 2.571e-220 | Upregulated |
BTG1 | BTG anti-proliferation factor 1 | 1.384 | 2.887e-117 | Upregulated |
EPAS1 | Endothelial PAS domain protein 1 | 1.231 | 4.725e-96 | Upregulated |
From the results of KEGG analysis (Fig. 4D), we noted that pathways were markedly enriched in HIF-1 signaling pathway, cytokine-cytokine receptor interaction and MAPK signaling pathway according to p-value.
Ppis Network Analysis And Hub Genes Identification
The protein-protein interactions of DEGs were constructed and visualized to discover their possible functional sites, which were environmentally determined or predicted. The first 100 differential proteins among HK-2 normoxia and HK-2 hypoxia were selected to establish the PPIs network, and the single differential protein that could not form the interaction relationship or only two nodes with one edge would be removed from the network. As can be seen from Fig. 5A, the most of these interactively proteins in this study were enhanced their expression responding to hypoxia. To find the potent pathogenic mechanisms of RIF in human, it is very important to explore the key proteins and their associating subnetworks in the complex PPIs network. It was apparent that vascular endothelial growth factor A (VEGFA), solute carrier family 2 member 1(SLC2A1), prolyl 4-hydroxylase subunit alpha 1 (P4HA1), phosphoglycerate kinase 1 (PGK1), BCL2 interacting protein 3 (BNIP3), pyruvate dehydrogenase kinase isozyme 1(PDK1), solute carrier family 2 member 3 (SLC2A3), glycogen synthase 1 (GYS1), aldolase fructose-bisphosphate A (ALDOA) and hexokinase 2 (HK2) were the most interacted genes from hub genes analysis (Figs. 5B and 5C). Therefore, combined with antecedent results (mentioned in the functional analysis), we regarded all 16 genes as potential goals and details were revealed in Table 1.
Results Of Chip Sequencing
Hypoxia is a common phenomenon among fibrosis and as a hypoxic transcriptional activator, HIF-1α can bind to the hypoxia-responsive element (HRE) to modulate genes expression 15. Therefore, to explore more predominant biomarkers, we next investigated the genes could be activated by HIF-1α directly. To do this, we developed chip sequencing and analyzed those genes interacting with HIF-1α based on DNA-protein interactions under hypoxia. All 2915 genes were identified (Fig. 8A) through aligning to the reference genome (hg38). According to peak score, we dissected 10 genes consisting of GNAS complex locus (GNAS), EF-hand calcium binding domain 10 (EFCAB10), Y-box binding protein 1 (YBX1), actin related protein 2/3 complex subunit 2 (ARPC2), inhibitor of DNA binding 2 (ID2), centromere protein B (CENPB), ENY2 transcription and export complex 2 subunit (ENY2), phosphatase and tensin homolog (PTEN), glutamate dehydrogenase 1 (GLUD1) and eukaryotic translation initiation factor 5A (EIF5A), which HIF-1α could regulate gene expression through promoter regions (Fig. 6). We also paid attention to peak distributions of all genes binding with HIF-1α listed in Table 2, and it is worthy to note that types of non-coding RNA and the most superior 20 of them were demonstrated in Table 3 and Table S2 separately, paving future directions for ncRNA. In addition, to verify the reliability of chip sequencing, we also tested it under normoxic conditions, which showed that there were only a few binding sites (Table S3 and S4 of Supplement 1), consistent with the rapid degradation of HIF-1α under normoxic conditions.
Table 2
Distributions of binding sites with HIF-1α under hypoxia
Annotation | Number of peaks | Total size (bp) | Log2 Ratio (obs/exp) | LogP enrichment (+ values depleted) |
3UTR | 179 | 26833139 | 2.831 | -204.921 |
TTS | 79 | 32404629 | 1.379 | -29.914 |
Promoter | 69 | 35946139 | 1.034 | -16.75 |
Exon | 510 | 37120946 | 3.873 | -939.718 |
Intron | 1182 | 1257910936 | 0.003 | -0.757 |
5UTR | 17 | 2601483 | 2.801 | -20.676 |
Intergenic | 825 | 1684358172 | -0.937 | 406.224 |
Table 3
Annotations of non-coding RNA binding with HIF-1α under hypoxia
Annotation | Number of peaks | Total size (bp) | Log2 Ratio (obs/exp) | LogP enrichment (+ values depleted) |
scRNA | 0 | 97 | 0 | 0 |
snoRNA | 0 | 357 | 0 | 0 |
rRNA | 0 | 25562 | -0.034 | 0.024 |
miRNA | 0 | 97618 | -0.126 | 0.092 |
ncRNA | 18 | 7044070 | 1.446 | -8.627 |
Results were then analyzed by GO and KEGG enrichment and the top 30 enrichment analysis of BP, CC and MF were shown in Figs. 7A, 7B and 7C, and top 10 of all subtypes were emphasized in Fig. 7D. From the biological processes of GO enrichment analyses, the following were enriched: macromolecule metabolic process, cellular macromolecule metabolic process and regulation of metabolic process. As can be seen from CC annotations, cell part, intracellular and intracellular organelle were the top components. Moreover, corresponding annotations of molecular functions included binding, protein binding and organic cyclic compound binding. Additionally, we found that Ras signaling pathway, RNA transport and long-term depression may be intensively involved in HIF-1α induced hypoxia according to Fig. 7E and 7F, which revealed enriched KEGG pathways according to p-value.
Genes Modulated By Hif-1α
Subsequently, we compared the RNA sequencing with the chip sequencing, and 43 genes of them were regarded as common targets. Potential genes from this module were further prioritized based on their functional analysis targeting to angiogenesis, and 4 of them were regarded as candidates to angiogenesis and HIF-1α concurrently (Fig. 8A). However, only 2 genes were finally defined as novel targets including VEGFA and BTG1considering feasible binding sites (Fig. 8B and C). The details of them and their binding sites were listed in Table 4.
Table 4
Annotated analysis of candidates
Gene | Gene description | Chr | Start | End | Strand | Peak score | Annotation |
VEGFA | Vascular endothelial growth factor A | Chr6 | 43781625 | 43782536 | + | 246.54 | exon (NM_001171626, exon 6 of 7) |
BTG1 | BTG anti-proliferation factor 1 | Chr12 | 92143996 | 92145768 | + | 95.16 | intron (NM_001731, intron 1 of 1) |