Comprehensive analysis of the long non-coding RNA expression profile and functional roles in a contrast-induced acute kidney injury rat model.

Long non-coding RNAs (lncRNAs) have been identified as a class of regulatory RNAs that participate in both physiological and pathological conditions, including acute kidney injury. However, the roles of lncRNA dysregulation in the pathogenesis of contrast-induced acute kidney injury (CI-AKI) are largely unknown. In the present study, the expression profiles of lncRNAs in kidney tissue were compared between rats with CI-AKI and controls using high-throughput RNA sequencing. In total, 910 differentially expressed (DE) lncRNAs (DElncRNAs), including 415 downregulated and 495 upregulated lncRNAs, were identified at 12 h after intra-arterial iodinated contrast medium injection (fold change ≥2; P<0.05). Eight DElncRNAs were further selected and validated using reverse transcription-quantitative polymerase chain reaction. A previous study defined microRNA (miRNA) and mRNA expression changes in the same CI-AKI model. In the present study, a lncRNA-mRNA co-expression network comprising 349 DElncRNAs and 202 DEmRNAs was constructed. The function of these DElncRNAs was mainly associated with oxidative stress and inflammation. Additionally, lncRNA-associated competing endogenous RNA (ceRNA) analysis revealed a network comprising 40 DElncRNA nodes, 5 DEmiRNA nodes and 59 DEmRNA nodes. Among which, the carnosine dipeptidase 1-specific and the transmembrane protein 184B-specific networks were likely to be associated with CI-AKI. The results of the present study revealed the expression profile and potential roles of lncRNAs in CI-AKI, and provide a framework for further mechanistic studies. Copyright: © Bao et al.

Introduction

Contrast-induced acute kidney injury (CI-AKI) is a frequent complication following intravascular administration of iodinated contrast medium (ICM), which is generally characterized by an increase in serum creatinine (SCr) of 0.5 mg/dl or a 50% relative elevation over baseline within 48 h following contrast medium exposure, in the absence of an alternative etiology (1). The development of this iatrogenic syndrome is associated with adverse early and long-term clinical outcomes (2,3). The introduction of safer contrast agents and optimized hydration strategies have led to a decline in the incidence of CI-AKI in the general population to 0.6-2.0% (4). However, the morbidity rate still remains as high as 20-50% in vulnerable subgroups, including those with chronic kidney disease, acute myocardial infarction or diabetes mellitus (1,5). Accumulating evidence has revealed a variety of pathways implicated in CI-AKI, which involve direct tubular cell toxicity, outer medullary ischemia, oxidative stress and inflammation (6,7). However, the precise mechanisms underlying the pathogenesis of CI-AKI remain largely unknown, leading to a lack of early diagnostic biomarkers or more efficient prevention strategies.

Long non-coding RNAs (lncRNAs) are a heterogeneous class of RNAs that lack a protein-coding capability and are >200 nucleotides in length (8). High-throughput RNA sequencing (RNA-seq) has led to the continuous discovery of lncRNAs, which are emerging as important regulators in a variety of physiological and pathological conditions, with a relatively tissue-specific expression manner (9-13). Mechanistically, lncRNAs have been proposed to function through cis or trans transcriptional regulation, organization of nuclear domains, and by acting as competing endogenous RNAs (ceRNAs) (14,15). The complicated regulatory mechanisms finally result in the formation of a large-scale regulatory network across the transcriptome, and have provided useful explanations of the pathological processes in various diseases. Over the past decades, several studies have discovered the vital roles of lncRNAs in acute kidney injury induced by a variety of etiologies (16-18). For instance, NEAT1 was identified to be involved in sepsis-induced kidney injury by targeting microRNA (miRNA) miR-204 and subsequently activating the nuclear factor-κB pathway (17). Recently, Cheng et al (19) explored the potential link between lncRNAs and CI-AKI. As reported, LNC_000343 potentially regulates the expression of the Kielin/chordin-like protein by acting as a ceRNA of rno-miR-1956-5p in CI-AKI rats (19). However, ICM was administered in an intravenous manner in this model, which was slightly different from the clinical practice of coronary angiography (20-22); moreover, samples were collected at 24 h after iohexol injection in that study (19). To the best of our knowledge, the SCr level is significantly elevated at 24 h following ICM intervention in rat models (23-25), indicating that the regulation of the damage response might take place even earlier. Thus, it is proposed that the lncRNA transcriptome changes in CI-AKI might vary between the different methods of modeling and time intervals of sample collection.

In the present study, a proven CI-AKI rat model with significant elevation of SCr and remarkable histopathological alterations was adopted (26). In this model, ICM was administered in an intra-arterial manner, and kidney samples were harvested at 12 h post ICM injection. By performing deep RNA-seq analysis, the expression profiles of lncRNAs were compared between the CI-AKI group and the control group in kidney tissues. The potential function of these differentially expressed (DE) lncRNAs (DElncRNAs) was then analyzed using bioinformatic algorithms. Based on the expression profiles of miRNAs and mRNAs identified previously in the same CI-AKI model (GSE130796 for miRNA; and GSE130795 for mRNA) (26), lncRNA-mRNA co-expression analysis was performed and an lncRNA-associated ceRNA network was constructed to better understand the regulatory roles of the DElncRNAs. Generally, the present study might provide new insights of the dysregulated lncRNAs involved in the emergence of CI-AKI complicated by clinical arteriography.

Materials and methods

Materials and animals

Indomethacin, N-ω nitro-L-arginine methyl ester (L-NAME), and pentobarbital sodium were purchased from Sigma-Aldrich; Merck KGaA. Iopromide (Ultravist 370; 370 mg/ml iodine) was obtained from Bayer AG. Phosphate buffer (pH 8.4) was synthesized by the experimental center. Indomethacin was dissolved in phosphate buffer (5 mg/ml), and L-NAME was dissolved in 0.9% normal saline (10 mg/ml), immediately before injection. In total, 18 3-month old male Sprague-Dawley rats, weighing ~300-400 g at the start of the experiment, were obtained from Fujian Medical University (Fuzhou, China). The 18 rats were kept in individual cages under controlled conditions of light (12-h light/dark cycle), temperature (21-23˚C) and humidity (50-60%), with free access to tap water and standard rat chow for a 7-day adaptive period.

Ethics statement and establishment of the CI-AKI rat model

The protocols of the animal experiments were conducted in accordance with the Guiding Principles in the Use of Animals in Toxicology (27), and were approved by animal experiment ethics review committees of 900 Hospital of the Joint Logistics Team, Chinese People's Liberation Army (approval no. IACUC-2017-17).

The study design was previously published (26). Briefly, rats were deprived of water for 48 h and then anaesthetized using pentobarbital sodium (40 mg/kg, i.p.). Catheters were placed in the right femoral vein and the common carotid artery (24 Gx21 mm; SPECATH) before a baseline arterial blood sample was drawn (1 ml) to determine the SCr. Indomethacin was then administered (10 mg/kg, i.v.), followed by L-NAME (10 mg/kg, i.v.) after 15 min. After another 15 min, the rats were randomized to receive iopromide (CI-AKI group, n=9) or normal saline (control group, n=9) via the carotid artery cannulation (7.8 ml/kg, i.a.). The rats were then allowed to recover in individual cages with free access to tap water and standard chow. A blood sample (1 ml) was obtained ~12 h after iopromide injection from the abdominal aorta under pentobarbital sodium anesthesia (40 mg/kg, i.p.) to measure postoperative SCr. The kidney tissue for RNA-sequencing or reverse transcription-quantitative polymerase chain reaction (RT-qPCR) validation was then harvested immediately and stored in RNAsafety Reagent (cat. no. 01901-50; http://www.shbio.com/products/3034) at 4˚C overnight, before being transferred to a -20˚C refrigerator. For histological analysis, part of the kidney tissue (~2-mm thick), containing both the cortex and medulla, was fixed in 10% neutral formalin liquid for 12-24 h at room temperature. The kidney tissue was then dehydrated through an ascending series of ethanol, infiltrated with acetone and embedded in paraffin at 60˚C. The paraffin blocks were cut into 5-µm-thick sections and subsequently dewaxed in xylene and rehydrated in a descending ethanol gradient. Hematoxylin and eosin (H&E) staining of the sections was performed with hematoxylin for 5 min and eosin for 2 min, both at room temperature. The Paller scores were calculated to determine the severity of tubular injury (13). Finally, the rats were sacrificed with an overdose of pentobarbital anesthesia (200 mg/kg, i.p.).

RNA extraction, qualification and purification

According to the manufacturer's instructions, total RNA was extracted using a mirVana™ miRNA Isolation kit (cat. no. AM1561; Ambion, Thermo Fisher Scientific, Inc.). The RNA Integrity Number (RIN) number was checked to inspect RNA integrity using an Agilent Bioanalyzer 2100 (Agilent Technologies, Inc.). Only those samples with a RIN ≥7.0 and 28s/18s ≥0.7 were identified as qualified and were selected for further analysis. The qualified total RNA was then purified using an RNAClean XP kit (cat. no. A63987; Beckman Coulter, Inc.) and RNase-Free DNase set (cat. no. 79254; Qiagen GmBH).

Library construction and RNA sequencing

Library construction and RNA sequencing were performed by Shanghai Biotechnology Corporation (http://www.shbio.com). Strand-specific cDNA libraries were prepared from ribosomal RNA-depleted RNAs using a VAHTS Total RNA-seq (H/M/R) Library PrepKit for Illumina (cat. no. NR603-02; Vazyme Biotech Co., Ltd.). The RNAs were first interrupted into shot fragments. Next, first-strand cDNA synthesis was performed using Oligo(dT)12-18 primers and SuperScript™ II Reverse Transcriptase kit (cat. no. 18064-014; Invitrogen, Thermo Fisher Scientific, Inc.). The simplified steps were as follows: i) Heat the mixture of total RNA, primers and dNTPs to 65˚C for 5 min; ii) add First-Strand Buffer (5X) and incubate at 42˚C for 2 min; iii) add SuperScript™ II Reverse Transcriptase and incubate at 42˚C for 50 min; iv) inactivate the reaction by heating at 70˚C for 15 min. Double-strand cDNA synthesis was then performed and the cDNA fragments were purified by Agencourt® AMPure XP Beads (cat. no. A63881; Beckman Coulter, Inc.). After final PCR enrichment, the cDNA libraries were quantified by Qubit® 2.0 Fluorometer (Invitrogen, Thermo Fisher Scientific, Inc.) and further validated by the Agilent 2100 system (Agilent Technologies, Inc.) to calculate the library concentration. Cluster was generated by the Illumina cBot system (version 02; Illumina, Inc.) with the library diluted to a loading concentration to ~10 pM. Following cluster generation, the libraries were sequenced in paired-end read (2x150 bp) mode on the Illumina HiSeq 2500 platform (Illumina, Inc.). TruSeq PE Cluster Kit v3 (cat. no. PE-401-3001; Illumina, Inc.), TruSeq SR Cluster Kit v3 (cat. no. GD-401-3001; Illumina, Inc.) and TruSeq Rapid Duo Dample Loading Kit (cat. no. CT-402-4001; Illumina, Inc.) were used in library clustering and sequencing. High-quality reads were aligned to the Rattus norvegicus reference genome (ftp://ftp.ensembl.org/pub/release-83/fasta/rattus_norvegicus/dna/Rattus_norvegicus.Rnor_6.0.dna_rm.toplevel.fa.gz) using the spliced mapping algorithm in the Hisat2 software (v 2.0.4) (28). The unmatched reads were analyzed using the gffcompare software (v 0.9.8) to predict novel lncRNAs (29). For the gene fragment calculation, the fragments per kilobase of transcript per million mapped reads (FPKM) value was determined using Stringtie and trimmed mean of M values algorithm. Based on the FPKM data, lncRNAs with a fold change ≥2 and P<0.05 were identified as differentially expressed.

RT-qPCR verification of selective DElncRNAs

In order to verify the reliability of the RNA-seq results, several DElncRNAs from the RNA-seq data were selected for validation in an independent cohort of 6 CI-AKI rats and 6 controls using RT-qPCR. Total RNA was extracted from rat kidneys using mirVana™ miRNA Isolation Kit (Ambion; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocols, and reverse transcribed using ReverTra Ace qPCR kit (Toyobo Life Science). The reverse transcriptional conditions were as follows: 37˚C for 15 min, 98˚C for 5 min and then kept at 4˚C. Next, qPCR was performed in triplicate using the 7500 Real-Time PCR System with Power SYBR®-Green PCR Master Mix (Thermo Fisher Scientific, Inc.). The thermocycling conditions were 50˚C for 2 min, then 95˚C for 10 min, followed by 40 cycles at 95˚C for 15 sec and 60˚C for 1 min. The Gapdh gene (glyceraldehyde-3-phosphate dehydrogenase) was used as internal control. The expression levels of the lncRNAs was normalized and quantified using the 2-ΔΔCq method (30). P<0.05 was considered to indicate a statistically significant difference. The following reverse primers were used: NONRATT027338.2 (forward, 5'-CAGGACAAAGGAACCCAGC-3' and reverse, 5'-CCGAGAGAGAGCAGCAATGA-3'); NONRATT027428.2 (forward, 5'-GCTGTAAATGTAGGCTATGGTGGTT-3' and reverse, 5'-AGCTCTCATGGGTTCTGTCATCTC-3'), NONRATT000173.2 (forward, 5'-ACCAAACAAGACCACCAGCAT-3' and reverse, 5'-GGAGGGACTGATGTGTACGAAAC-3'), NONRATT005775.2 (forward, 5'-CCTCCCACCCTCTGATGTAG-3' and reverse, 5'-AGAAAGTGCTCGTGGACAGG-3'), NONRATT016226.2 (forward, 5'-GAACCAGAGGATGGCGACA-3' and reverse, 5'-GATGGCATGAAGGGATGAAT-3'), NONRATT018005.2 (forward, 5'-CCTTCTCCTTCCAGATAACTTACACA-3' and reverse, 5'-GTGACTGCCAGGGTGCTAAAC-3'), NONRATT023682.2 (forward, 5'-CAGGTGCCTCCTCTCAGTCAA-3' and reverse, 5'-CCTCACCCCCCTAGTCTTCTTAA-3'), MSTRG.22041.2 (forward, 5'-TGCACTGAGCAGGACTGAAAA-3' and reverse, 5'-TTATCCCTTTGCATTCACTCCAA-3'), Gapdh (forward, 5'-TGGCCTCCAAGGAGTAAGAAAC-3' and reverse, 5'-GGCCTCTCTCTTGCTCTCAGTATC-3').

DElncRNA target gene prediction

Both the cis and trans regulation analyses were applied to predict the target genes of lncRNAs. The coding genes located within 10 kb upstream and downstream of the lncRNAs were predicted as putative cis-targets. The trans-prediction were performed using RNAplex (version 2.4.14) (31).

Co-expression analysis of DElncRNAs and DEmRNAs

Pearson's correlation coefficient (PCC) between the newly identified DElncRNAs and previously identified DEmRNAs was calculated (P<0.05). The DEmRNAs were identified as having a fold change ≥2 and P<0.05. The co-expression network of DElncRNAs-DEmRNAs was presented using Cytoscape software (version 3.8.0) (32).

Construction of the lncRNA-associated ceRNA network

The significant mature DElncRNAs (fold change ≥2; q<0.05), DEmiRNAs; fold change ≥1.5; P<0.05), and DEmRNAs (fold change ≥2; q<0.05), were selected for ceRNA analysis. The lncRNA-miRNA interactions and miRNA-mRNA interactions were predicted using the miRanda database according to their shared miRNA-binding seed sequence sites (33). The lncRNA-miRNA pairs and miRNA-mRNA pairs were then integrated into a ceRNA network in the comprehensive analysis. After that, the PCCs of the expression values between DElncRNAs and DEmRNAs, as well as the associations among all the three kinds of RNA (DElncRNAs, DEmiRNAs and DEmRNAs), were determined according to the expression levels. The lncRNA-associated ceRNA network was further visualized using Cytoscape (32).

Functional enrichment analysis

The functional enrichment analysis of putative targets of lncRNAs were conducted using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) software (34). Gene ontology (GO) term analysis was performed in terms of functional classification, including biological process, cellular component and molecular function (35). Kyoto encyclopedia of genes and genomes (KEGG) analysis was further conducted to reveal the association of these genes with different pathways (36). To account for multiple comparisons, both Benjamini and Bonferroni corrections were used. P<0.05 was considered to indicate a statistically significant difference.

Statistical analysis

The SPSS 22.0 software (SPSS Inc.) was used for statistical analysis. Data are presented as the mean ± SD. The comparison between groups was conducted by unpaired Student's t-test. P<0.05 was considered to indicate a statistically significant difference.

Results

ICM exposure induced remarkable injury to the kidneys

The CI-AKI rats exhibited a significantly increased SCr, with an average level of 59.9±23.0%, as well as pronounced histopathological changes to the kidney, as observed using H&E staining, compared with those of the controls (26).

Characteristics of the RNA-seq data

Total RNA in the kidney tissue derived from 3 paired CI-AKI rats and 3 controls was isolated and analyzed. RNA-seq produced over 90 million raw reads, with clean read ratios ranging from 94.43 to 95.45%. Most of the clean reads could be mapped perfectly to the Rattus reference genome (Table I). StringTie (v1.3.0) and gffcompare (v0.9.8) were used to assemble and quantify the transcripts, respectively (29,37). The identification of lncRNAs was then performed using three tools (Pfam, CPC and CNCI) (38-40). Finally, a total of 21,248 lncRNAs were identified (Table SI; Fig. S1), including 12,081 (56.9%) exonic sense-overlapping lncRNAs, 4,584 (21.6%) intergenic lncRNAs, 1,900 (8.9%) intronic sense-overlapping lncRNAs, 1,281 (6.0%) bidirectional lncRNAs, 1,012 (4.8%) exonic antisense lncRNAs and 390 (1.8%) intronic antisense lncRNAs (Fig. 1A). The features of the lncRNAs were also analyzed. The lncRNAs tended to have a lower expression level (calculated by the FPKM method) and a shorter transcript length compared with those of the protein-coding transcripts (Fig. 1B and C). Moreover, the majority of lncRNAs were identified to have fewer exons compared with mRNAs (Fig. 1D). The raw data has been uploaded to the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE134824).

Table I

Summary of draft reads of cDNA libraries.

Sample	Raw reads	Clean reads	Clean ratio, %	Mapped reads	Mapping ratio, %
Control_1	119,917,792	113,890,762	94.97	101,521,071	89.14
Control_2	116,702,770	110,756,482	94.90	98,219,953	88.68
Control_3	118,362,004	112,430,138	94.99	99,295,268	88.32
Case_1	113,431,152	107,112,640	94.43	95,644,335	89.29
Case_2	103,095,464	98,408,905	95.45	88,173,472	89.60
Case_3	115,941,526	110,000,513	94.88	98,694,544	89.72

Figure 1

Characteristics of the identified lncRNAs and protein-coding RNAs. (A) The number and proportion of each type of lncRNA, including exonic sense-overlapping lncRNAs, intergenic lncRNAs, intronic sense-overlapping lncRNAs, bidirectional lncRNAs, exonic antisense lncRNAs, and intronic antisense lncRNAs. (B) Comparison of expression levels between lncRNAs and mRNAs. (C) Comparison of transcript lengths between lncRNAs and mRNAs. (D) Exon number distribution of lncRNAs and mRNAs. The abscissa is the exon number, and the ordinate is the proportion of transcripts. lncRNA, long non-coding RNA; FPKM, fragments per kilobase of transcript per million mapped reads.

Differential expression analysis of lncRNAs

According to the screening criteria, a total of 910 DElncRNAs, including 74 novel ones, were identified in the CI-AKI rats compared with the controls (Table SII). Among them, 495 DElncRNAs were upregulated and 415 were downregulated. Notably, NONRATT027876.2 was the most highly expressed, with a fold-change of 603.8, while NONRATT010468.2 was the most downregulated, whose level decreased by 95.2-fold. Some of the DElncRNAs had no expression in either the CI-AKI rats or the controls, resulting in an infinite fold change value; for instance, MSTRG.22041.2. Uniformly, these lncRNAs were considered as dysregulated with a fold change of 2.0. The DElncRNAs were further illustrated using a heatmap (Fig. 2).

Figure 2

Global profiling of lncRNAs in the contrast-induced acute kidney injury group and control group. (A) Hierarchical clustering analysis of the differentially expressed lncRNAs. Each group includes three duplicates. A total of 910 dysregulated lncRNAs were identified, including 495 upregulated and 415 downregulated ones (fold change ≥2.0; P<0.05). Colors from green to red represent the lncRNA expression abundance from poor to rich. (B) Scatter plot presenting differences in lncRNA expression. lncRNA, long non-coding RNA.

RT-qPCR validation of DElncRNAs

To validate the reliability of the RNA-seq data, eight dysregulated lncRNAs were selected for further RT-qPCR analysis in an independent cohort of six CI-AKI rats and six controls. All the candidate lncRNAs were verified to be differentially expressed (Fig. 3 and Table SIII; P<0.05), which was consistent with the RNA-seq data. These results confirmed the authenticity of the high-throughput sequencing data.

Figure 3

RT-qPCR validation analysis of lncRNAs. RT-qPCR verification of the expression profiles of eight lncRNAs was conducted in the kidneys of six CI-AKI rats and six control rats. *P<0.05 vs. control group. RT-qPCR, reverse transcription-quantitative PCR; lncRNA, long non-coding RNA; CI-AKI, contrast-induced acute kidney injury.

Trans and cis regulatory prediction and functional enrichment analysis of DElncRNAs

The associations between DElncRNAs and their potential targets were predicted according to both cis- and trans-acting patterns (Tables SIV and V). Location analysis identified 805 DElncRNAs that were relatively close to 756 protein-coding genes, indicating of a cis-regulatory manner. GO analysis was performed and 76 GO terms were significantly enriched (Fig. S2A and Table SVI; P<0.05). KEGG analysis revealed seven significant pathways (Fig. S2B and Table SVII; P<0.05). Notably, ‘p53 signaling pathway’ (rno04115), which plays a critical role in acute kidney injury (41), was proposed to be associated with CI-AKI. In the functional analysis of the trans DElncRNA targets, 30 GO terms and 5 KEGG pathways were significantly enriched (Fig. S3A and B, Tables SVIII and SIX; P<0.05). Some of these GO terms were associated with energy metabolism and apoptosis (GO:0005739, GO:0006919, GO:0045335, GO:0005777 and GO:0045730).

DElncRNA-DEmRNA co-expression analysis

A co-expression network of the DElncRNAs and DEmRNAs was constructed (Table SX). This network comprised 1,632 lncRNA-mRNA interactions, containing 349 DElncRNA nodes and 203 DEmRNA nodes (PCC>0.990; P<0.05). Finally, the DElncRNAs or DEmRNAs with significant expression levels were selected and visualized, including 46 DElncRNAs and 38 DEmRNAs (Fig. 4). GO functional enrichment analysis of the 203 differentially co-expressed genes revealed 35 GO terms between the CI-AKI and control group. Among them, 20 GO terms were biological process-associated, 9 were cellular component-associated, and 6 were molecular function-associated (Fig. 5A and Table SXI; P<0.05). Interestingly, a considerable number of significant oxidative stress-related GO terms were identified, and oxidative stress has been recognized as a hallmark of CI-AKI (6). Notably, 17 dysregulated genes were remarkably enriched in the GO term ‘oxidation-reduction process’ (GO:0055114), including Ugdh (encoding UDP-glucose 6-dehydrogenase), Cmah (encoding cytidine monophospho-N-acetylneuraminic acid hydroxylase), Glud1 (encoding glutamate dehydrogenase 1) and Cyp2f4 (encoding cytochrome P450, family 2, subfamily f, polypeptide 4). The aforementioned genes and associated DElncRNAs were expected to modulate oxidative reactions during CI-AKI. In addition, several GO terms were inflammation-associated (GO:0006955, GO:0002504, GO:0034341, GO:0002474 and GO:0042130). KEGG analysis identified 19 remarkable pathways (Fig. 5B and Table SXII; P<0.05). Similarly, most of the significant pathways were relevant to the pathological processes of inflammation. These results suggested that the dysregulated lncRNAs might play a role in CI-AKI pathogenesis by various mechanisms.

Figure 4

Co-expression network of the DElncRNA-DEmRNA interactions. DElncRNAs or DEmRNAs with an average fragments per kilobase of transcript per million mapped reads value >2.0 in either the contrast-induced acute kidney injury group or the control group were selected. Triangles represent lncRNAs and squares represent mRNAs. Light red indicates upregulated and light green indicates downregulated. Pearson's correlation coefficient, >0.990. lncRNA, long non-coding RNA; DElncRNA, differentially expressed lncRNA; DEmRNA, differentially expressed mRNA.

Figure 5

Functional enrichment analysis of DEmRNAs involved in the DElncRNA-DEmRNA co-expression network. (A) Top 15 terms of the Gene Ontology enrichment analysis. Red bars represent biological process, green bars represent cell component, and blue bars represent molecular function. (B) Top five pathways of the Kyoto Encyclopedia of Genes and Genomes enrichment analysis. DEmRNAs, differentially expressed mRNAs; DElncRNA, differentially expressed long non-coding RNA.

The DEmRNAs involved in the co-expression network were further compared with predicted targets regulated by DElncRNAs in a cis or trans manner. Significantly, 25 co-expressed DEmRNAs were found to be potentially targeted by 19 cis-acting DElncRNAs (Table SXIII; P<0.05). The cis-acting DElncRNA-DEmRNA pairs potentially involved in CI-AKI were displayed (Table II). It was also found that some trans-acting DElncRNAs were predicted to target several co-expressed differential genes. However, the majority of these DElncRNAs and DEmRNAs had extremely low expression levels, making them unsuitable for further studies.

Table II

Cis-regulatory DElncRNA-DEmRNA co-expression pairs likely to be involved in CI-AKI.

DElncRNA	Fold change	Location	Strands	DEmRNA	Fold change	Location	Strands	r
NONRATT025711.2	3.72	7:29479206-2948005	+	Slc5a8	4.25	7:29435443-29477947	+	0.999
NONRATT016225.2	3.65	2:88404185-88408022	-	Lrrcc1	2.36	2:88384795-88414012	-	0.992
NONRATT016226.2	4.72	2:88410568-88413958	-					0.970
NONRATT016224.2	3.66	2:88402887-88402723	-					0.967
NONRATT016768.2	5.00	2:210803868-210805276	-	Gstm1	3.37	2:210803868-210809306	-	0.991
NONRATT014530.2	3.87	19:22621344-22632069	-	Gpt2	3.89	19:22590880-22626653	-	0.986
NONRATT007430.2	2.24	12:30160774-30165923	+	Asl	2.03	12:30165693-30178341	+	0.985
NONRATT026655.2	2.53	7:29959585-29961154	-	Gas2l3	2.41	7:29959596-29986163	-	0.971
NONRATT027381.2	-2.19	7:145147281-145149879	-	Ppp1r1a	-2.17	7:145146480-145154131	-	0.956

DEmRNA, differentially expressed mRNA; DElncRNA, differentially expressed long noncoding RNA; CI-AKI, contrast-induced acute kidney injury.

lncRNA-associated ceRNA network

LncRNAs can also function as molecular decoys for miRNAs through miRNA response elements and thus regulate gene expression as a ceRNA network (14). Thus, the lncRNA-associated ceRNA network was constructed to further illustrate the role of DElncRNAs in CI-AKI development. First, the candidate DElncRNAs and mature DEmiRNAs with significant connections were matched (P<0.05). Next, the selected DEmiRNAs were matched to their putative target DEmRNAs (P<0.05). The DElncRNA-DEmiRNA network and the DEmiRNA-DEmRNA network were then integrated into a ceRNA regulatory network (Table SXIV). Meanwhile, the associations among the expression levels of DElncRNAs, DEmiRNAs and DEmRNAs were calculated to further test the possibilities of these ceRNA pairs. The ceRNA pairs with negative associations between DElncRNAs and DEmRNAs were selected and visualized in a final ceRNA network, which covered 22 DElncRNAs, 5 DEmiRNAs and 37 DEmRNAs (Fig. 6). As shown in Fig. 7, the DEmRNAs within the network were significantly enriched in seven GO terms and seven KEGG pathways (Tables SXV and SXVI; P<0.05), and most of these processes were inflammation-associated. According to the ceRNA hypothesis, lncRNA levels generally correlate negatively with miRNA levels, and the miRNAs subsequently downregulate the expression of target genes. The present study identified that the downregulated lncRNAs, such as NONRATT021928.2 and NONRATT025462.2, were ceRNAs of the upregulated rno-miR-126a-5p, rno-miR-200a-5p and rno-miR-322-5p, which target Cndp1 (encoding carnosine dipeptidase 1). The downregulated lncRNAs NONRATT020679.2 and NONRATT023587.2 were ceRNAs of the upregulated rno-miR-126a-5p, rno-miR-200a-5p and rno-miR-322-5p, which target Tmem184b (encoding transmembrane protein 184B). The ceRNA pairs that are likely to be associated with CI-AKI are shown in Table III.

Figure 6

lncRNA-associated competing endogenous RNA network of the DElncRNA-DEmiRNA-DEmRNA pairs. Triangles represent DElncRNAs, circles represent DEmiRNAs, and squares represent DEmRNAs. Light red indicates upregulated, and light green indicates downregulated. DEmRNAs, differentially expressed mRNAs; DElncRNA, differentially expressed long non-coding RNA; DEmiRNAs, differentially expressed microRNAs.

Figure 7

Functional enrichment analysis of differentially expressed mRNAs involved in the long non-coding RNA-associated competing endogenous RNA network. (A) Significant GO terms enriched in the GO analysis. Red bars represent biological process, green bars represent cell component, and blue bars represent molecular function. (B) Significant pathways enriched in the Kyoto Encyclopedia of Genes and Genomes enrichment analysis. GO, Gene Ontology.

Table III

Significant ceRNA pairs likely to be associated with CI-AKI.

DElncRNA	Fold change	P-value	DEmiRNA	Fold change	P-value	DEmRNA	Fold change	P-value
NONRATT025462.2	-14.27	<0.001	rno-miR-126a-5p	2.60	<0.001	Cndp1	-2.16	<0.001
NONRATT021928.2	-37.10	<0.001
NONRATT025462.2	-14.27	<0.001	rno-miR-200a-5p	2.24	<0.001
NONRATT021928.2	-37.10	<0.001
NONRATT025462.2	-14.27	<0.001	rno-miR-322-5p	2.41	<0.001
NONRATT021928.2	-37.10	<0.001
NONRATT004914.2	-12.83	<0.001
NONRATT001846.2	-55.97	<0.001
NONRATT020679.2	-8.43	<0.001	rno-miR-126a-5p	2.60	<0.001	Tmem184b	-2.08	<0.001
NONRATT023587.2	-12.2	<0.001
NONRATT020679.2	-8.43	<0.001	rno-miR-200a-5p	2.24	<0.001
NONRATT023587.2	-12.2	<0.001

ceRNA, competing endogenous RNA; CI-AKI, contrast-induced acute kidney injury; lncRNA, long non-coding RNA; DElncRNA, differentially expressed lncRNA; DEmiRNA; differentially expressed microRNA; DEmRNA; differentially expressed mRNA.

Discussion

The present study was designed to compare the lncRNA expression patterns of CI-AKI rats and controls from a global perspective. To that end, the transcriptome of the kidney tissue was generated using high throughput sequencing in the present study. A total of 910 DElncRNAs were identified, including 74 novel ones, that were dysregulated in CI-AKI rat kidney at 12 h following intra-arterial iopromide exposure. Importantly, DElncRNA-DEmRNA co-expression and functional enrichment analysis indicated that these lncRNAs were mainly involved in oxidative stress and inflammation reactions. Moreover, bioinformatic algorithms revealed the CNDP1-specific and Tmem184b-specific networks that might be associated with CI-AKI. The aforementioned results might provide new insights into the lncRNA-associated regulatory mechanisms underlying CI-AKI.

Depending on the CI-AKI model, striking differences were observed in the expression profiles of lncRNAs between CI-AKI rats and controls, as well as their functional characteristics. DElncRNA-DEmRNA co-expression and functional enrichment analyses implied that these DElncRNAs were likely to be associated with several key aspects of CI-AKI, including oxidative stress and inflammation. In particular, oxidative stress is recognized as one of the most significant mechanisms of CI-AKI (42). Accumulating evidence has shown that the injection of contrast medium might induce hypoxia and consequently augment reactive oxygen species generation in the kidney (43,44). Enhanced oxidative stress can cause damage to membrane lipids, cellular proteins and DNA. Membrane lipid peroxidation might change the cellular and mitochondrial membrane permeability and activate specific cell signaling cascades of cell death, ultimately resulting in renal parenchymal oxidative injury (42,45). Several lncRNAs, including PRINS, MIR210HG, linc-ATP13A4-8 and linc-KIAA1737-2, have been identified to mediate hypoxic kidney injury (18,46). In the present study, several significant oxidative stress-related GO terms were identified in the functional analysis of the DEmRNAs involved in the DElncRNA-DEmRNA co-expression network, including ‘oxidation-reduction process’ (GO:0055114), ‘removal of superoxide radicals’ (GO:0019430), and ‘oxidoreductase activity’ (GO:0016491). Moreover, several DElncRNA-DEmRNA co-expression pairs we found with potential cis-regulatory relationships. In particular, the NONRATT016768.2 locus was within the locus of the associated gene Gstm1 (encoding glutathione S-transferase mu 1). Gstm1 had been recognized as a target of nuclear factor erythroid-derived 2-like 2 and played a protective role in response to oxidative stress in kidney disease (47,48). In addition, lncRNA-associated ceRNA analysis revealed a significant Cndp1-specific network that was also likely to be associated with the process of antioxidation. Cndp1 encodes a serum carnosinase that is associated with carnosine degradation. A decreased level of CNDP1 might result in an increased carnosine concentration, and could lead to renoprotective effects by eliminating peroxyl and hydroxyl radicals (49,50). These results suggested a vital role of lncRNAs in oxidative stress-mediated injury in CI-AKI. The DElncRNAs implicated in oxidative stress are worthy of further in-depth functional study to help find indications regarding the important regulators of CI-AKI pathogenesis.

Previous studies have also indicated a role of inflammation in the pathogenesis of CI-AKI (51-55). Andreucci et al (55) showed that nuclear factor-κB (NF-κB) is activated together with an increase in interleukin-8 in sodium diatrizoate-treated cells. An in vivo study also suggested that NF-κB induced a pro-inflammatory response, along with renal tubular damage, and inhibition of NF-κB using sodium butyrate could protect the kidney from contrast-induced injury (51). In addition, the nucleotide-binding oligomerization domain-like pyrin domain containing protein 3 inflammasome activated the pro-inflammatory cytokines and apoptotic pathway, and mediated kidney injury following contrast medium exposure, both in vitro and in vivo (54). Moreover, the inflammatory process of CI-AKI also involves the regulatory effects of macrophage (7). In the present study, functional analysis revealed that a fairly large number of DElncRNAs were involved in the inflammation process. For instance, previous studies had revealed that Gpnmb (encoding glycoprotein nmb) was expressed in macrophages, negatively regulated inflammation, and played a role in acute kidney injury (56,57). The present study showed that a group of DElncRNAs (such as MSTRG.11448.4, MSTRG.2420.1, MSTRG.6245.1 and NONRATT023367.2) were linked to Gpnmb in the co-expression analysis. The inflammation-associated DElncRNA-DEmRNA pairs further supported the view that the CI-AKI pathological process is accompanied by inflammation.

The mitogen-activated protein kinase (MAPK) pathway regulates inflammation, apoptosis, and oxidative stress, and is involved in kidney injury (58,59). A previous study showed that Tmem184b might activate the MAPK signaling pathway (60). The present study identified a Tmem184b-specific ceRNA network involving downregulated Tmem184b, two downregulated DElncRNAs (NONRATT020679.2 and NONRATT023587.2) and two upregulated DEmiRNAs (rno-miR-126a-5p and rno-miR-200a-5p). These aforementioned dysregulated RNAs possessed remarkable expression levels and significant fold changes; the functional roles of these RNAs in CI-AKI should be further explored.

It should also be highlighted that the dysregulation of several DElncRNAs and their putative targets likely resulted in a reno-protective effect. Indeed, the activation of the reparative response occurred at an early stage following injury. The recovery process, especially tubular recovery, is vital to protect from progression of chronic kidney disease after injury (61). Therein, the cross-talk between tissue-reparative macrophages and tubular cells is an important mechanism underlying tubular recovery after kidney injury. In addition, Gpnmb is essential to switch macrophages towards a wound-healing phenotype (57). The upregulation of Gpnmb, together with the co-expressed DElncRNAs, was potentially implicated in CI-AKI by modulating macrophage polarization and tubular recovery. Besides, Gstm1 and Cndp1 in the bioinformatic networks were also associated with a protective reaction. These results suggest that more attention should be paid to the recovery process in research on CI-AKI.

Despite the significant findings in the present study, several limitations should be emphasized. Firstly, the RNAseq analysis was performed using total RNA of kidneys, not on specific cells. The specific roles of the distinct cells in kidney tissue still remain unclear. Secondly, subcellular location is a key feature for understanding a lncRNA's function; however, this was not analyzed in the present study. Thirdly, these dysregulated RNAs were identified in rats, the importance of these RNAs in clinical samples is yet to be elucidated. Finally, the levels of mRNAs do not necessarily associate with levels of their translation products; however, the translation products of the critical DEmRNAs had not been analyzed in animal models or clinical samples. Intensive studies should be conducted to confirm the preliminary results in the present study and further reveal the functional characteristics of lncRNAs in CI-AKI development.

In conclusion, the present study identified the lncRNA expression profile in vulnerable rats following intra-arterial administration of contrast medium. The DElncRNAs and their targets were likely to be associated with several key aspects of CI-AKI, including oxidative stress and inflammation. The present study provides an insight into the roles of lncRNAs in the pathogenesis of CI-AKI.