This article is mentioned in:

Introduction

Chronic kidney disease (CKD) is increasingly recognized as a worldwide public health issue (1,2), and chronic hypoxia often occurs among patients with CKD. Renal hypoxia is emerging as a key player by influencing tubular epithelial cells during the process of the acute kidney injury (AKI)-to-CKD transition. Under hypoxia, these injured cells fail to redifferentiate, which results in aggravating hypoxia and gives rise to a vicious circle (3).

Some studies have shown that human renal proximal tubular cells (namely HK-2 cells) may be the primary target of a hypoxic insult (4). Although it is clear that HK-2 cells demand high concentrations of oxygen to maintain their metabolic functions (5–7), the relevant mechanisms have not yet been identified. Hypoxia can regulate the expression of a wide variety of genes that may be induced or suppressed by transcriptional or post-transcriptional mechanisms (8–10). Hypoxia response elements have been identified in the regulatory regions of a number of genes and contain consensus binding sites for the transcription factor hypoxia-inducible factor-1α (HIF-1α) (11–13). Several studies have demonstrated changes in collagen gene expression in various mesenchymal cells exposed to hypoxia (14–17).

Renal proximal tubular cells are a primary target for hypoxic injury in the kidneys and some studies have reported a distinct cohort of 48 genes with a closely shared hypoxia-dependent expression profile, identified by microarray analysis (18,19). Since the global transcriptional mechanisms underlying these events are not yet fully understood, in this study, we undertook an unbiased approach using RNA sequencing (RNA-Seq) to define the transcriptomic responses regulated by hypoxia in HK-2 cells. RNA sequencing technology provides major advantages over microarray technology, providing more precise information on absolute transcript levels and transcript variants, and currently unannotated transcribed regions can be analyzed; therefore, it is increasingly used in various biological applications (20). However, to date, and to the best of our knowledge, no RNA-seq assays demonstrating transcriptomic changes have been reported in HK-2 cells under hypoxic conditions.

In this study, using RNA-seq analysis, we examined HK-2 cells under hypoxic and normoxic conditions to clarify global changes in gene expression under such conditions. The expression levels of selected genes were further validated by reverse transcription-quantitative PCR (RT-qPCR). Although the expression levels of the vast majority of genes remained unaltered, a specific cluster of genes with a closely shared expression profile were regulated undr hypoxic conditions. By using stringent bioinformatics analysis of the data, our findings extend those of previous and related publications on gene expression patterns that may be critical for the development of CKD.

Materials and methods

Cell culture and exposure to hypoxia

The HK-2 cell line was purchased from the Cell Bank of Type Culture Collection of Chinese Academy of Sciences (Shanghai, China). The HK-2 cells were cultured in high-glucose DMEM supplemented with 10% fetal bovine serum (FBS), 100 U/ml penicillin and 100 µg/ml streptomycin (both from Thermo Fisher Scientific, Waltham, MA, USA) in humidified air containing 5% CO2 at 37°C. After the cells doubled and redoubled, the viability of the cells was at its maximum, and this phase can be identified as the logarithmic growth phase. Cells in the logarithmic growth phase were used in all of the experiments.

The cells were seeded in 6-well plates at a density of 5×105 cells/well and cultured for 24 h. To mimic hypoxic conditions, the cells in medium were then incubated under low oxygen conditions. The oxygen concentrations were maintained at 1–3% using a compact gas oxygen controller, which was held under positive pressure in an atmosphere of 94–92% N2/5% CO2/1–3% O2 for 24 h. The control cells were cultured under normoxic conditions (5% CO2 for 24 h).

Western blot analysis

The cells were plated in 24-well plates and cultured in DMEM supplemented with 10% FBS for 24 h. Following exposure to hypoxia for 24 h, the cells were lysed with immunoprecipitation assay buffer containing protease inhibitors. The cell lysates were centrifuged at 14,000 × g for 10 min at 4°C followed by incubation on ice. The protein concentrations of the lysates were examined using the Bradford protein assay kit (Vazyme Biotech Inc., Nanjing, China). The cell lysates were boiled and separated by sodium dodecyl sulfate polyacrylamide gel electrophoresis and transferred onto polyvinylidene difluoride membranes via semi-dry transfer (Bio-Rad Laboratories Inc., Hercules, CA, USA). The membranes were washed in Tris-buffered saline containing 0.1% Tween-20 (TBST), blocked with 5% non-fat milk in TBST for 1 h at room temperature, and incubated with a primary rabbit monoclonal antibody against HIF-1α (Cat. no. ab179483; dilution 1:5,000; Abcam, Cambridge, UK) or glyceraldehyde-3-phosphate dehydrogenase (GAPDH) antibody (Cat. no. 10494; dilution 1:5,000; Proteintech Technology, Inc., Wuhan, China) overnight at 4°C. The membranes were washed 3 times in TBST, followed by incubation with the appropriate horseradish peroxidase-linked secondary anti-rabbit antibodies (Cat. no. SA00001-2; dilution 1:5,000; Proteintech Technology, Inc.) for 1 h at room temperature. The specific proteins on the blots were developed by enhanced chemiluminescence (ECL; Vazyme Biotech, Inc.) and visualized as the bands on the CL-XPosure Film (Thermo Fisher Scientific). The optical densities of the bands were measured on the GS710 Densitometer and analyzed using Quantity One image analysis software (Bio-Rad Laboratories Inc.).

RNA extraction and RNA-Seq analysis

Total RNA was isolated from cells using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. Total RNA sample quality and integrity were controlled using a Bioanalyzer 2100 (Agilent, Boeblingen, Germany). For library preparation, 5 µg total RNA were captured by Dynabeads Oligo(dT) sheared to fragements of ~200 bp, and reverse transcribed using the SuperScript III cDNA Synthesis kit, as previously described (21) (both from Life Technologies, Grand Island, NY, USA). cDNA was end-repaired, A-tailed and ligated to Illumina sequencing adapters and amplified by PCR. Library preparation were performed using the TruSeq RNA LT V2 Sample Prep kit, as previously described (22) (Illumina, San Diego, CA, USA). The sequencing library was qualified by Qubit 2.0 (Life Technologies) and Bioanalyzer 2100 (Agilent), then sequenced on a Illumina HiSeq 2000 with 2×100 bp paired-end sequencing, which were controlled by HiSeq Control software. Raw reads in FASTQ format were aligned to the human genome (GRCh37/hg19) using TopHat software 2.0.10, as previously described (23). This parameter defines -G genes.gtf -r 0–mate-std-dev 80–solexa1.3-qual. Transcript isoform assembly and abundance estimation were performed using Cufflinks (v1.3.0) (24) and combined with gene annotations from the National Center for Biotechnology Information (ftp://igenome:G3nom3s4u@ussd-ftp.illumina.com/Homo_sapiens/UCSC/hg19/Homo_sapiens_UCSC_hg19.tar.gz). The cDNA libraries were sequenced using the HiSeq™ 2000 sequencing platform (Illumina).

Validation by qPCR

qPCR were executed to validate the data obtained by RNA-seq analysis. Total RNA was isolated from the HK-2 cells using the RNeasy Mini kit and an RNase-free DNase set (both from Vazyme Biotech, Inc.). qPCR was performed to determine the expression of genes in a SYBR-Green PCR Master mix (Vazyme Biotech, Inc.) using the StepOnePlus™ real-time PCR detection system (Step One Plus 2.1 software) with universal thermal cycling parameters. The primer sequences are listed in Table I. GAPDH was used as an internal control for target genes for reaction efficiency. The ΔΔCt method was used to determine the relative amounts of product.

Table I Primers used in qPCR.

Table I

Primers used in qPCR.

Gene name	Primer sequences	Length (bp)
COL4A2	F: 5′-TGCTACCCGGAGAAAGGAG-3′	106
	R: 5′-CTTTGCGGCCCTGTAGTCC-3′
COL7A1	F: 5′-CCAGAGGTCGTTCCGGAG-3′	81
	R: 5′-GCTCTTCCCACTTCGACC-3′
THBS3	F: 5′-ACAGTTCTCCTGCGACTCCG-3′	172
	R: 5′-GCATCCTCAAATACGCCTTC-3′
LAMA5	F: 5′-ACCCAAGGACCCACCTGTAG-3′	169
	R: 5′-TCATGTGTGCGTAGCCTCTC-3′
LAMB1	F: 5′-TGGCTGAAGTGGAACAGCTCTC-3′	95
	R: 5′-TGTCTTCAACAGAATGTCTTCAGCA-3′
ITGB4	F: 5′-AGAGGGAGGAAGAGGATGGC-3′	168
	R: 5′-GCAGTAGGCGCAGTCCTTAT-3′
MUC1	F: 5′-CGTCATGGACATTGATGGTACC-3′	228
	R: 5′-GGTACCTCCTCTCACCTCCTCCAA-3′
VEGF	F: 5′-GAGTACATCTTCAAGCCATCCTG-3′	203
	R: 5′-TGCTCTATCTTTCTTTGGTCTGC-3′
GLUT1	F: 5′-AAGGTGATCGAGGAGTTCTACA-3′	119
	R: 5′-ATGCCCCCAACAGAAAAGATG-3′
PLIN2	F: 5′-TGAGATGGCAGAGAACGGTGTGAA-3′	84
	R: 5′-TTGCGGCTCTAGCTTCTGGATGAT-3′
TGF-β1	F: 5′-GAGCCTGAGGCCGACTACTA-3′	130
	R: 5′-CGGAGCTCTGATGTGTTGAA-3′
PDGFB	F: 5′-TGATGCCGAGGAACTATTCATCT-3′	178
	R: 5′-TTTCTTCTCGTGCAGTGTCAC-3′
HIG2	F: 5′-CCACAGTGCAAGACTCCATC-3′	150
	R: 5′-GCCATACTGCTGAGGAAAGC-3′
BNIP3L	F: 5′-TCGAGCCGCCGCCGCCCCTG-3′	138
	R: 5′-CATTGCCATTATCATTGCCATTG-3′
JUNB	F: 5′-GTCACCGAGGAGCAGGAGG-3′	63
	R: 5′-TCTTGTGCAGATCGTCCAGG-3′
MMP-1	F: 5′-GATGAAGTCCGGTTTTTCAAAG-3′	71
	R: 5′-GGGGTATCCGTGTAGCACCAT-3′
GLUT3	F: 5′-TTCGTCTCTAGCCTGCACTG-3′	79
	R: 5′-ACACAACTTCTCCGGGTGAC-3′
GAPDH	F: 5′-GGAAGGTGAAGGTCGGAGTCA-3′	90
	R: 5′-GCAACAATATCCACTTTACCAGAGTTAA-3′

[i] F, forward; R, reverse.

Gene Ontology (GO) and pathway analysis for differentially expressed genes (DEGs)

GO functional classifications were analyzed at the macroscopic level based on cellular component, biological process and molecular function categories. GO terms with a P<0.01 by Fisher's exact test were considered enriched. Pathway analysis was determined by the DEGs according to the KEGG, BioCarta and Reactome databases. Significant pathways were selected in accordance with Fisher's exact test followed by Benjamini Hochberg multiple testing correction, and the threshold of significance was defined as a value of P<0.05.

Identification of alternative splicing (AS) events and the discovery of novel mRNA transcripts

Mixture of isoforms (MISO) analysis (25) was used to examine differentially regulated exons across samples. MISO analysis was accomplished using paired-end reads according to the workflow given (25). The reads alignment files produced by TopHat and the pre-build human genome alternative events and were downloaded from the MISO reference manual page (http://genes.mit.edu/burgelab/miso/docs/#gff-event-annotation); the date of access for the databases was April 30, 2015.

Identification of single nucleotide polymorphisms (SNPs)

The identification of SNPs was successively performed using the mutation detection software GATK, as previously described (26,27).

Results

Expression of HIF-1α in HK-2 cells under hypoxic conditions

The HK-2 cells were exposed to the hypoxic milieu (1% O2 for 24 h) to examine the effects of oxygen deprivation on HK-2 cells, and the expression of HIF-1α was examined by western blot analysis. The results revealed that hypoxia upregulated the expression level of HIF-1α (Fig. 1A and B).

Figure 1 Expression of hypoxia-inducible factor-1α (HIF-1α) in human renal proximal tubular cells (HK-2) cells under hypoxic conditions and Volcano plots in the magnitude of gene expression ratios. (A) Western blot analysis of HIF-1α protein. (B) Ratios of the expression of HIF-1α relative to that of GAPDH. Hypoxia involved exposure of the cells to 1% O2 for 24 h. (C) Volcano plot of differentially expressed genes in HK-2 cells under normoxic and hypoxic conditions. The log2 fold change difference between the hypoxic samples and normoxic samples is represented on the x-axis, and the negative log of P-values is represented on the y-axis. The differentially expressed genes were plotted in blue for the upregulated genes and red for downregulated genes, and non-significant genes are shown as gray points. (D) Validation of selected genes based on RNA-Seq were verified using qPCR. The relative expression levels of differentially expressed genes are shown (relative value to that of normoxia, value of normoxia assigned values of 1). Results were analyzed by an unpaired Student's t-test, **P<0.01.

Illumina sequencing and overview of the sequence reads

We obtained a global overview of the transcriptome from the two cDNA libraries separately constructed from HK-2 cells under hypoxic and normoxic conditions by using the Illumina HiSeq™ 2000 sequencing platform. The two libraries were named after the hypoxic group and normoxic group, respectively. A total of 111,829,528 raw reads was obtained from the cDNA libraries. After the raw reads were filtered, a total of 93,002,921 high-quality mappable reads remained, of these, 51,182,720 reads were from cells under hypoxic conditions and 41,820,201 reads were from cells under normoxic conditions. The sequencing reads were mapped to the human genome (GRCh37/hg19); at least 86% of the sequence bases could be aligned to the genome, and >77% of the mapped read bases were localized in mRNA regions. In addition, the Q20 (those reads with an average quality score >20) was >96%, and the GC content was consistently ~53% for both hypoxic and normoxic conditions, which suggested that the sequencing was highly accurate.

An overview of the sequencing process is presented in Table II. Of the total number of clean reads, 87.06% mapped to multiple (1.69%) or unique (85.92%) genome locations for hypoxic conditions, and 86.89% mapped to multiple (1.79%) or unique (85.11%) genome locations for normoxic conditions. Thus, the read mapping evenness across the transcripts of two samples was good. Volcano plots in the magnitude of gene expression ratios are displayed in Fig. 1C and demonstrate the significance of differences in gene expression between the two samples.

Table II Summary of clean Illumina RNA-seq reads and Illumina transcriptome sequencing in HK-2 cells under hypoxic and normoxic conditions.

Table II

Summary of clean Illumina RNA-seq reads and Illumina transcriptome sequencing in HK-2 cells under hypoxic and normoxic conditions.

Reads category	Hypoxic HK-2 cells	Normoxic HK-2 cells
Total raw reads	61,284,796	50,544,732
Total clean reads	58,425,234	48,128,048
Total mapped reads	51,182,720 (87.06%)	41,820,201 (86.89%)
Multiple mapped reads	984,919 (1.69%)	859,998 (1.79%)
Uniquely mapped reads	50,197,801 (85.92%)	40,960,203 (85.11%)
Q20 (%)	96.70 (%)	96.92 (%)
GC content (%)	53.29 (%)	53.29 (%)

[i] HK-2, human renal proximal tubular cells.

Evaluation of DEGs

To examine the differences in gene expression patterns between the hypoxic group and normoxic group, the criteria for the detection of differential expression gene in this study was the absolute value of fold change >2 and a higher statistical significance value (P-value <0.05). A total of 279 genes was observed to be differentially expressed between the hypoxic group and normoxic group (data not shown). Of these genes, 201 were specifically upregulated, and 78 genes were markedly downregulated in the hypoxic group vs. the normoxic group (data not shown). These results demonstrated that the number of upregulated DEGs was considerably higher than the number of downregulated DEGs. Of the 201 upregulated genes, it is notable that a host of genes are associated with hypoxia, such as encode transforming growth factor (TGF)-β1, vascular endothelial cell growth factor (VEGF), 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 4 protein (PFKFB4) and glucose transporter protein (GLUT1). The above-mentioned data also demonstrated that the result of RNA-Seq in this study was reliable.

Validation of gene expression

A total of 17 genes was selected and quantified to verify the RNA-Seq data by qPCR. As shown in Figs. 1D and 4B, all of the 17 genes were confirmed. The results illustrated that the expression patterns of these genes were highly consistent with the RNA-Seq data, which suggested that our transcriptome analysis was reliable.

Figure 4 Extracellular matrix (ECM)-receptor interaction pathway. (A) Upregulated genes are colored in red. (B) Expressions of the COL4A2, ITGB4, LAMA5, LAMB1 and THBS3 were identified by qPCR (relative value to that of normoxia, value of normoxia assigned values of 1). Results were analyzed by unpaired Student's t-test, *P<0.05 and **P<0.01.

GO and KEGG pathway functional enrichment analysis of the DEGs

GO enrichment and pathway analysis of the DEGs were accomplished to determine the potential biological functions of the DEGs. The enriched GO terms are shown in Fig. 2 according to 3 top-level ontologies, namely biological process (BP), cellular component (CC) and molecular function (MF). With respect to biological processes, the upregulated genes were mainly enriched for extracellular matrix metabolism, extracellular structure organization and locomotion, and the downregulated genes were mostly related to ventricular cardiac muscle cell differentiation, regulation of cellular amino acid metabolic process and the positive regulation of cardiac muscle contraction. For molecular function, the upregulated genes were invovled in protein binding, glycoprotein binding and growth factor binding, and the downregulated genes were associated with diphosphotransferase activity, growth factor activity and MHC II protein complex binding. According to the cellular component annotation, the upregulated genes were mainly localized in the extracellular matrix (ECM) including proteinaceous ECM and basement membrane, while the downregulated genes were enriched in the cytoplasm, macromolecular complex and DNA-directed RNA polymerase I complex.

Figure 2 Gene Ontology (GO) analysis of differentially expressed genes (DEGs) in human renal proximal tubular cells (HK-2) cells under hypoxic and normoxic conditions. The top 10 (A) upregulated and (B) downregulated GO terms for biological processes, molecular function and cellular component were generated by a bar plot with mostly significant P-values. BP, biological process; CC, cellular component and MF, molecular function.

For KEGG pathway enrichment analysis, the up- and downregulated genes appeared to be significantly enriched in ECM-receptor interaction (Figs. 3 and 4A) and antigen processing and presentation, respectively. This result is consistent with that of the enriched biological process, and pathway-based analyses help to further understand the biological functions of genes. The 11 DEGs involved in the ECM-receptor interaction pathway included collagen type IV alpha 2 (COL4A2), heparan sulfate proteoglycan 2 (HSPG2), integrin subunit alpha 3 (ITGA3), integrin subunit beta 3 (ITGB3), integrin subunit beta 4 (ITGB4), integrin subunit beta 8 (ITGB8), laminin subunit alpha 3 (LAMA3), laminin subunit alpha 5 (LAMA5), laminin subunit beta 1 (LAMB1), thrombospondin 3 (THBS3) and agrin (AGRN); some of these genes were confirmed by qPCR (Fig. 4B).

Figure 3 KEGG pathway analysis of differentially expressed genes (DEGs) in human renal proximal tubular cells (HK-2) cells under hypoxic and normoxic conditions. (A and B) The top 10 upregulated and downregulated pathways, respectively. The data shown are the negative log2 P-values within each category.

Among the upregulated genes, platelet derived growth factor subunit B (PDGFB), solute carrier family 2 member 1 (SLC2A1), TGFA, TGF-β1 and VEGFA were enriched in the pathway of renal cell carcinoma (RCC). Furthermore, PDGFB, TGF-β1 and VEGFA were validated by qPCR. These genes dispalyed differential expression patterns in the HK-2 cells under hypoxic conditions. Of the downregulated genes, ornithine decarboxylase 1 (ODC1) and spermidine synthase (SRM) were enriched in the pathway of glutathione metabolism. In brief, the DEGs provide a better understanding of gene expression patterns in HK-2 cells under hypoxic conditions.

Identification of AS events and the discovery of novel mRNA transcripts

In this study, we obtained a complete picture of AS events between the hypoxic group and normoxic group by MISO (25). All theoretical splicing junctions were identified, and a total of 51,387 and 49,144 AS events were found in the hypoxic group and normoxic group, respectively. In this study, AS events identified contain 5 different types: alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), exon skipped (ES), intron retention (IR) and mutually exclusive exon (MXE) (data not shown). As shown in Table III, the frequency of each AS event was highly coincident between the hypoxic group and normoxic group. Of the 5 types of AS events, ES was the most prevalent, and was comprised of 27,574 (53.66%) and 26,375 (53.67%). A total of 9,394 novel transcripts coding potential was found in the two samples (data not shown).

Table III Numbers of alternative splicing (AS) events in 5 classes identified in HK-2 cells under hypoxic and normoxic conditions

Table III

Numbers of alternative splicing (AS) events in 5 classes identified in HK-2 cells under hypoxic and normoxic conditions

AS type	HK-2 cells under hypoxic conditions (%)	HK-2 cells under normoxic conditions (%)
A5SS	5,729 (11×.15)	5,495 (11.18)
A3SS	8,290 (16.13)	7,898 (16.07)
ES	27,574 (53.66)	26,375 (53.67)
IR	4,594 (8.94)	4,401 (8.96)
MXE	5,200 (10.12)	4,975 (10.12)
Total	51,387	49,144

[i] A5SS, alternative 5′ splice site; A3SS, alternative 3′ splice site; ES, exon skipped; IR, intron retention; MXE, mutually exclusive exon; HK-2, human renal proximal tubular cells.

Identification of SNPs

SNPs are valuable markers that are widely used in genetics and evolution, and RNA-Seq technology has the potential to detect SNPs rapidly and reliably. To understand the changes in SNPs in HK-2 cells under hypoxic conditions, we identified a total of 69,706 SNPs (data not shown). Among these, 36,747 SNPs were from the hypoxic group and 32,959 SNPs from the normoxic group. The types of SNPs included synonymous, missense, stopgain, exonic, ncRNA, UTR5, UTR5 and UTR3, UTR3, intronic and intergenic (Table IV). The identified SNPs, as molecular markers, will be extremely helpful for the genetic linkage mapping of HK-2 cells under hypoxic conditions. In brief, the transcriptome data will provide a valuable resource for the future development of molecular markers and functional genomics.

Table IV Number of different types of SNPs identified by RNA sequencing.

Table IV

Number of different types of SNPs identified by RNA sequencing.

Sample	Hypoxic group	Normoxic group
Total	36,747	32,959
Synonymous	5,465	5,070
Missense	3,885	3,583
Stopgain	5	7
Exonic	9,488	8,788
NcRNA	2,099	1,859
UTR5	1,762	1,642
UTR5 and UTR3	7	7
UTR3	9,564	8,728
Intronic	8,505	7,995
Intergenic	4,037	2,930

[i] SNPs, single nucleotide polymorphisms.

Discussion

Renal hypoxia is a major factor in the pathophysiology of AKI to CKD transition (3). To maintain normal metabolism, renal proximal tubular cells require a large amount of oxygen to function properly, and renal proximal tubular cells may be the primary target of a hypoxic insult (7). Although the hypoxia of HK-2 cells has been widely known in last decades, the comprehensive study of HK-2 cells under hypoxic conditions has not yet been reported, at least to the best of our knowledge. Many hypoxia-related genes have been identified in different cells, including TGF-β1, VEGF, PFKFB4 and GLUT. The underlying molecular events however, remain unknown. Leonard et al found 48 upregulated genes in HK-2 cells exposed to hypoxia for 8, 16 and 36 h with the aid of cDNA microarrays (19). In this study, the expression of HIF-1α was increased in the hypoxic HK-2 cells compared to the nomoxic cells, which confirmed that HK-2 cells were indeed cultured in a hypoxic environment. Therefore, we then performed a comprehensive analysis of the transcriptomes of HK-2 cells under hypoxic and normoxic conditions by RNA-Seq. Compared with a previous micorray study (19), we identified a total of 279 DEGs, and 201 DEGs (>2-fold change) were upregulated in the HK-2 cells exposed to hypoxic conditions for 24 h. The level and number of DEGs we obtained were not completely consistent with the previous microarry study (19). The different exposure times to hypoxia may have resulted in this discrepancy.

All DEGs were classified into GO and KEGG categories to provide an overview of their biological functions and associated biochemical pathways. Many studies have shown the expression of individual hypoxia-related genes in different cells (28–30). In the present study, our analysis indicated that the upregulated genes were highly enriched in ECM organization (GO:0030198), protein binding (GO:0005515), extra ventricular cardiac muscle cell differentiation (GO:0055012), diphosphotransferase activity (GO:0016778) and the cytoplasm (GO:0005737). Analysis for pathways revealed that ECM-receptor interaction was the most enriched upregulated pathway, including the genes related to the formation of the basement membrane. Of these genes, the COL4A2 gene encodes one of the 6 subunits of type IV collagen, the major structural component of basement membranes (31). Our data based on mRNA sequencing and qPCR confirmed the upregulated expression of subunits of type IV collagen gene under hypoxic conditions. LAMA3, LAMA5, LAMB1 belong to the laminin family for formation and function of the basement membrane (32). HSPG2 encodes the perlecan protein, which is a major component of basement membranes (33). Chronic hypoxia often occurs in the kidney tissues of many patients with CKD (18). Tubulointerstitial fibrosis leads to thickness of the tubular basement membrane and accumulation of interstitial ECM (34-37). In this study, hypoxia upregulated the expression of genes related to the basement membrane. Hypoxia upregulated the expression of matrix metalloproteinase-1 (MMP-1) and MMP-1 breaks down the interstitial collagens, types I, II and III (38). These results aid in the understanding of the molecular regulatory mechanisms of tubulointerstitial fibrosis and CKD.

Specifically, we found the upregulated pathway of RCC and the downregulated pathway of glutathione metabolism in HK-2 cells under hypoxic conditions. A previous study demonstrated that patients with CKD may have an increased risk of RCC (39), and another study showed that RCC may be associated with CKD in patients, although the mechanisms involved are unclear (40). In this study, our data suggest that hypoxia may be an important factor in promoting the development of RCC in patients with CKD. A previous study demonstrated that HIG2 protein was involved in RCC for molecular targeted therapy (41). In this study, HIG2 mRNA upregulation induced by hypoxia provides another clue for the interaction between CKD and RCC. To the best of our knowledge, this study is the first report that DEGs induced by exposure to hypoxia upregulated the pathway of RCC in HK-2 cells, as shown by RNA-seq. Glutathione may directly and indirectly participate in reactive oxygen (e.g., H2O2) through enzymatic reactions (42). Glutathione deficiency contributes to oxidative stress, which plays a key role in the pathogenesis of many diseases, including cancer, inflammation, kwashiorkor, seizure, Alzheimer's disease, Parkinson's disease, sickle cell anemia, liver disease, cystic fibrosis, attack, stroke and diabetes (43). The present study is the first to describe that hypoxia downregulated the pathway of glutathione metabolism in HK-2 cells by RNA-seq. Some studies have shown that the THBS3 gene, as a stimulator of tumor progression in osteosarcoma (44–46). This study confirmed that the expression of THBS3 may be associated with hypoxia in HK-2 cells under hypoxia by RNA-seq.

This study of the transcriptome analysis enhances our knowledge of the existing gene annotation and changes in expression patterns under hypoxic conditions, and supplements AS events, novel transcripts and SNP identification. AS is an efficient way for genomes to encode additional transcripts, and previous studies have shown that AS plays a role in source of protein diversity (47,48). Novel transcripts and identification of AS patterns will contribute to a better understanding of the mechanisms of transcriptional regulation. As regards SNPs in HK-2 cells under hypoxic and normoxic conditions, their greatest importance in biomedical research is for comparing regions of the genome between cohorts in genome-wide association studies. SNPs without an observable impact on the phenotype (so-called silent mutations) are still useful as genetic markers in genome-wide association studies (49).

In conclusion, this study provides a comprehensive analysis of transcriptome in HK-2 cells under hypoxic and normoxic conditions. Comparison of the transcriptomes between two groups revealed a number of DEGs, and potential candidates involved in CKD, and further functional analysis of these genes may help us to elucidate the mechanisms responsible for the development of CKD. This study establishes a solid foundation for future genetic and functional genomics studies and expression analysis of genes in CKD.

Acknowledgments

This study was financially supported by the National Nature Science Foundation of the People's Republic of China (no. 81370868), the Fundamental Research Funds for the Central Universities and Jiangsu Province Scientific Research Innovation Project for Graduate Students (no. KYLX_0197).

Abbreviations:

References