Differential gene and lncRNA expression in the lower thoracic spinal cord following ischemia / reperfusion-induced acute kidney injury in rats

We used high-throughput RNA sequencing to analyze differential gene and lncRNA expression patterns in the lower thoracic spinal cord during ischemia/reperfusion (I/R)-induced acute kidney injury (AKI) in rats. We observed that of 32662 mRNAs, 4296 out were differentially expressed in the T8-12 segments of the spinal cord upon I/R-induced AKI. Among these, 62 were upregulated and 34 were downregulated in response to I/R (FDR < 0.05, |log2FC| > 1). Further, 52 differentially expressed lncRNAs (35 upregulated and 17 downregulated) were identified among 3849 lncRNA transcripts. The differentially expressed mRNAs were annotated as “biological process,” “cellular components” and “molecular functions” through gene ontology enrichment analysis. KEGG pathway enrichment analysis showed that cell cycle and renin-angiotensin pathways were upregulated in response to I/R, while protein digestion and absorption, hedgehog, neurotrophin, MAPK, and PI3K-Akt signaling were downregulated. The RNAseq data was validated by qRT-PCR and western blot analyses of select mRNAs and lncRNAs. We observed that Bax, Caspase-3 and phospho-AKT were upregulated and Bcl-2 was downregulated in the spinal cord in response to renal injury. We also found negative correlations between three lncRNAs (TCONS_00042175, TCONS_00058568 and TCONS_00047728) and the degree of renal injury. These findings provide evidence for differential expression of lncRNAs and mRNAs in the lower thoracic spinal cord following I/R-induced AKI in rats and suggest potential clinical applicability.


INTRODUCTION
Acute kidney injury (AKI) is a major renal disease with increasing incidence and mortality [1]. Clinically, AKI occurs due to renal or extra-renal surgery-induced ischemia/reperfusion (I/R), sepsis, and nephrotoxicity [2,3]. The current primary therapeutic strategies for AKI, which includes renal replacement, symptomatic relief and supportive treatment, still result in high mortality rates of 40-80% [4,5]. Some AKI patients cannot be treated effectively due to complications. Therefore, there is an urgent need to develop new treatment strategies that can significantly improve outcomes.
The efferent and afferent sympathetic nerves play a crucial physiological role in regulating renal function [6][7][8][9][10][11][12][13]. Activation of renal sympathetic nervous system has been implicated in ischemic AKI [14,15]. Larkin  Research Paper blood flow and resulted in acute renal failure [16]. Our previous study showed that the T9 spinal cord segment was primarily involved in sympathetic regulation of renal function [17]. In addition, the excitability of renal sympathetic nervous system increased during I/R-induced acute renal failure [18,19]. Therefore, blocking of renal sympathetic nervous activation or denervated kidneys before renal ischemia ameliorated post-ischemic renal injury to some extent [20][21][22]. Previous studies also demonstrated that neuromodulation therapy benefitted heart failure patients in preclinical and small-sized clinical studies [23,24]. Hence, we postulated that spinal cord played an important role in acute kidney injury.
In recent years, significant progress has been made in understanding the function of renal sympathetic nervous system in the pathophysiology of AKI. However, the details regarding spinal cord involvement in AKI are unknown. Hence, we used RNA-seq with a higher sequencing depth as a systems biology approach in a rat model of I/R-induced AKI to identify the full transcriptome of the lower thoracic segments of spinal cord. This novel molecular biological technique allows identification of new protein-coding transcripts and novel non-coding RNA transcripts that play critical roles in many biological processes [25][26][27][28]. Since full transcriptome analysis of spinal cord from I/R-induced AKI is not available, we performed whole transcriptome sequencing and subsequent bioinformatic analysis to identify changes in mRNA and lncRNA expression in the lower thoracic segments of spinal cord in response to renal injury.

Validation of I/R-induced acute kidney injury in rats
The functional analysis and histopathological validation was performed to ensure that all I/R-induced AKI rats used for sequencing experiments showed significant acute renal failure. We observed that 45 min of ischemia resulted in increased serum creatinine (302.5 ± 49.8μM in I/R; 27.2 ± 1.6μM in sham; P<0.01; Figure 1A), serum BUN (42.2 ± 9.9mM in I/R; 5.6 ± 0.9mM in sham; P<0.01; Figure 1B). Also, the I/R rats showed damaged tubular cells in the PAS stained tissue sections and high histological injury scores compared to sham rats ( Figure 1C-1E). These results verified I/Rinduced AKI in the rats used in our study.

High throughput RNA-seq and genome-wide read mapping
High throughput RNA-seq was used to identify novel molecular players in the spinal cord that regulate renal function in the rat I/R-induced AKI model. Lllumina TrueSeq libraries were generated from total RNA isolated from the T8-12 of spinal cord tissues from 3 I/R and sham group rats, respectively. The flow chart of the sequencing strategy and analysis is shown in Supplementary Figure 1. To ensure accuracy of the subsequent bioinformatic analysis, low-quality reads and rRNA sequences were filtered prior to RNA-seq. We obtained over 85M of high quality sequence reads in each sample, which were mapped using Refseq (https://www.ncbi.nlm.nih.gov/ refseq/). The mapped reads were assembled into putative transcripts by the TopHat software analysis (http://tophat. cbcb.umd.edu/) (Supplementary Table 1). The quality of overall transcriptome including saturation, duplication and coverage was assessed as shown in Supplementary  Figure 2. Nearly 84.73%-91.47% of the reads mapped to the rat genome; about 79.43%-85.61% of the reads mapped to unique genomic regions among the aligned fragments, which were further verified for reliability. After aligning the sequences and identifying spliced junctions, we obtained 71906 novel transcripts.

The expression profiling of mRNAs in spinal cord after I/R-induced AKI
The mRNA expression was quantified using RPKM value (reads per kilobase per million mapped reads) to identify the gene expression changes in T8-12 spinal cord upon I/R induced AKI. Among the 32662 mRNAs that were measured by the RNA-seq, 2772 mRNAs were found deregulated upon I/R induced AKI by 2-fold, of which 1524 mRNAs were upregulated and 1248 mRNAs downregulated. Further analysis demonstrated 62 upregulated and 34 downregulated mRNAs in I/R (FDR<0.05, |log2FC|>1). The scatter and volcano plots analyses displaying the expression signatures of mRNAs are shown in Figure 2. The top 10 up-regulated and downregulated mRNAs are listed in Tables 1 and 2.

The expression profiling of lncRNA in spinal cord after I/R-induced AKI
To explore the potential role of lncRNAs in the T8-12 spinal cord at 24 hours after I/R-induced AKI, the expression profiles of lncRNAs were determined by high throughput RNA-seq. Our analyses revealed 3849 lncRNA transcripts in the spinal cord from I/R and control groups, among which 2253 known and 1596 novel lncRNAs were identified by filtering the transcripts with CPC (Coding Potential Calculator), CNCI (Coding-Non-Coding Index) and Pfam analyses (Supplementary Figure 3). Among these, 35 upregulated and 17 downregulated lncRNAs were identified in the spinal cord responding to I/Rinduced AKI (fold change >2, P < 0.05). Figure 3 show the expression signatures of lncRNAs by using scatter and volcano plots. The top 10 upregulated and downregulated lncRNAs were listed in Table 3.

Gene ontology annotation for differential expression genes
The differentially expressed genes identified by the RNA-seq analyses were annotated using the GO database (Gene Ontology, http://www.geneontology.org/) into three biological functional groups, namely, biological process ( Figure 4A), cellular component ( Figure 4B) and molecular function ( Figure 4C). We found that the differentially expressed mRNAs were primarily involved in the biological processes followed by cellular component GO functions.

Gene ontology and KEGG Pathway enrichment analysis
We performed GO and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses of deregulated genes to determine the molecular changes in the lower spinal cord tissue during I/R induced AKI. The top five enriched GO biological processes for upregulated genes in the I/R group included responses to stress, external biotic stimuli, external stimuli, killing of cells and disruption of cells ( Figure 5A). The most enriched GO cellular components for up-regulated genes in I/R group included hemoglobin complex, haptoglobinhemoglobin complex, extracellular space and DNA polymerase complex ( Figure 5B). The enriched GO molecular functions for up-regulated genes in I/R group were hemoglobin binding, damaged DNA binding and antioxidant activity ( Figure 5C).
The enriched GO biological processes for downregulated genes in I/R group were regulation of multicellular organismal process, response to organic substance, single-organism process and developmental process ( Figure 5D). The enriched GO cellular  components for down-regulated genes in I/R group included proteinaceous extracellular matrix, extracellular matrix, extracellular matrix part and extracellular space ( Figure 5E). The enriched GO molecular function for down-regulated genes in I/R included signaling receptor activity ( Figure 5F).
Next, we conducted KEGG pathway analysis for differentially expressed genes. The results showed that the up-regulated genes in I/R group were involved in cardiac muscle contraction, cell adhesion, cell cycle, reninangiotensin system, and serotonergic synapse ( Figure 6A). The downregulated genes in I/R group were involved in protein digestion and absorption, hedgehog signaling, neurotrophin signaling, MAPK signaling, and PI3K-Akt signaling pathways ( Figure 6B).

Real-time quantitative PCR and western blot validation
Then, we validated the high throughput RNA-seq results by performing qRT-PCR analysis of differentially expressed mRNAs and lncRNAs. Although the fold changes in mRNA and lncRNA were different between qRT-PCR and RNA-seq data, the expression patterns revealed similar conclusions ( Figure 7A, 7B). Our RNA-seq results showed the genes of Bax, Caspase-3 and Akt were upregulated and Bcl-2 was downregulated in I/R group. So we using Western blot validated these differential expression genes at protein level, which showed that P2X7R, S100A9, Bax and P-Akt were upregulated and Bcl-2 was downregulated in the lower   spinal cord following I/R-induced AKI ( Figure 7C-7E). These data indicated differential expression of the mRNAs and lncRNAs in the lower spinal cord tissue suggesting their involvement in the I/R-induced AKI. The qRT-PCR and western blot analyses were consistent with the high throughput RNA-seq data.

Differential lncRNA expression is correlated with the degree of renal injury
As indicated previously, many lncRNA transcripts in the spinal cord were either upregulated or downregulated during I/R induced AKI. Compared with control group, we found lncRNA TCONS_00018621, lncRNA TCONS_00034035, lncRNA TCONS_00034216 and lncRNA TCONS_00047107 were significantly increased (P<0.05), and lncRNA TCONS_00042175 and lncRNA TCONS_00058568 were significantly decreased in AKI group (P<0.05), shown in Figure 8A. Further, whether these differential lncRNAs expression associated with degree of renal injury, we even built two different animal models, such as kidney ischemia for 45 mins (I/R 45min) and 60 mins (I/R 60min), which were considered two different degrees of pathological injury. The results showed that expression of lncRNA TCONS_00042175, TCONS_00058568, TCONS_00047728 and TCONS_00034216 was lower at 60 mins compared to 45 mins post-I/R (P<0.05 or P<0.001), and the expression of lncRNA TCONS_00018621 was higher in I/R 60min group than the I/R 45min group ( Figure 8B). Interestingly, the expression of lncRNA TCONS_00034126 was 4-fold lower at 60 mins compared to at 45 mins (P<0.05, Figure 8B), however, compared with control group, the lncRNA TCONS_00034126 was obviously up-expressed in I/R group. Collectively, these results suggested that selected lncRNAs were either positively or negatively associated with the renal pathology injury induced by I/R.

DISCUSSION
Next generation sequencing techniques are useful for detecting novel genes, transcriptional and epigenetic networks. Our study is the first to explore the whole transcriptome profiles of spinal cord in a I/R rat AKI model using next generation RNA sequencing and bioinformatics analysis. We identified many differentially expressed genes, pathways and biological processes in the spinal cord during renal injury. By high throughput RNA-seq, we identified 32662 mRNAs in the spinal cord tissue from both the I/R and sham groups, of which 62 upregulated and 34 downregulated mRNAs (FDR<0.05,|log2FC|>1) were identified in response to renal ischemia. Similarly, 3849 lncRNAs of which, 2253 were known and 1596 were novel were also identified by the RNA-seq analysis. Among these, 35 up-regulated and 17 down-regulated mRNAs were identified upon renal ischemia. From this list, we confirmed and verified few differentially expressed mRNAs and lncRNAs by qRT-PCR and western blotting analyses.
Previous studies have suggested that the lower thoracic segment of spinal cord controls renal function [13,17,29]. However, the exact segments of the spinal cord that are involved and the mechanisms operating in the spinal cord as a result of renal ischemia/reperfusion injury are unclear. Many studies have shown that lncRNAs play  [27,30]. In this study, we found differential expression of many mRNAs and lncRNAs in the spinal cord further suggesting that the T8-12 spinal cord segment was involved in the neuronal response to renal injury. Although the functions of most lncRNAs are not fully known, our findings provide novel insights into their involvement in the molecular mechanism of renal failure.
Among the differentially expressed mRNAs, some including S100A9, RatNP-3b, SMOC2 and P2X7R have been reported earlier. Mitchell et al. reported that S100A8 and S100A9 expressing neutrophils traffic to the spinal cord during peripheral tissue inflammation, whereas, the S100A8 and S100A9 expression rapidly increased in response to carrageenan-induced inflammation of rat hind paws [31]. Sárvári et al demonstrated RatNP-      [32]. Recently, Roy et al. analyzed the transcriptome in hypothalamus and cerebral cortex by RNA-seq analysis and found that SMOC2 expression was related to domestication [33]. Rudqvist et al. showed downregulation of Slc47a2 in the thyroid tissue by oligonucleotide microarray, 24h after I-131 administration in rats [34]. Previous studies showed enhanced expression of P2X7R in the spinal cord dorsal horn in a rat model of neuropathic pain induced by chronic constriction injury [35].
We further verified that 3 lncRNAs (TCONS_00042175, TCONS_00058568, TCONS_ 00047728) were negatively associated with the degree of renal ischemia injury (Figure 8), they further demonstrated involvement of development process of renal injury induced by I/R. Recently, lncRNAs were shown to be part of gene regulatory networks responding to nerve injury, neuropathic pain, renal fibrosis and acute rejection after renal allografts [36][37][38][39]. Further, lncRNA-ATB was identified as a novel biomarker of acute kidney rejection, which could identify patients with acute rejection and predict loss of kidney function [39]. Tu et al. reported that silencing lncNONRATT021972 decreased the P2X7R levels in the cervical sympathetic ganglia and improved cardiac function after myocardial ischemia [40,41]. P2X7R expression was also strongly associated with kidney and nervous system diseases [42][43][44][45]. In this study, the P2X receptors P2X7R was upregulated in spinal cord after I/R-induced AKI group, consistent with previous studies that demonstrated increased P2X receptor expression in the stellate ganglia (SG) and superior cervical ganglia (SCG) after myocardial ischemia [46][47][48]. Moreover, several studies reported that the neuropeptide ghrelin or cholinergic agonists GTS-21 attenuated renal ischemia-reperfusion injury or sepsis-induced acute kidney injury through the vagus nerve [49][50][51]. Gigliotti et al and Inoue et al demonstrated that prior ultrasound prevented AKI by modulating the splenic neuroimmune axis, or stimulating the splenic cholinergic anti-inflammatory pathway, or vagus nerve stimulation [52][53][54]. Hence, we postulate that the list of novel lncRNAs (Table 3) found in our study represent new, yet undeciphered mechanisms in the response to renal injury.
Further, GO term enrichment analyses demonstrated that stress responses were the most enriched biological processes during I/R namely, response to external biotic stimuli and external stimuli. These results are consistent with previous studies showing activation of the afferent renal sympathetic nervous system during renal I/R [55]. The afferent renal nerves project to the ipsilateral dorsal horn in laminate I, III-V within the spinal cord [56]  and play an important role in mediating renal ischemic injury-induced reduction in renal hemodynamic and renal functions. The most enriched GO molecular functions were DNA binding, antioxidant and signaling receptor activities. We found increased Bax and P-Akt in the spinal cord tissue following renal I/R. Conversely, we also observed decreased Bcl-2, but Akt expression was unaltered. These indicated that pro-apoptotic and PI3K/Akt signaling pathways were upregulated in the spinal cord in response to renal injury. Additionally, KEGG analysis showed upregulation of cell cycle and renin-angiotensin system genes, whereas, the downregulated genes during I/R included those involved in protein digestion and absorption, neurotrophin, MAPK and PI3K/Akt signaling pathways. These findings are consistent with previous studies that demonstrated activation renin-angiotensin system [57], nerve growth factor promoted growth arrest and apoptosis in tubular renal cells [58] after renal injury. The data collectively indicated that inhibition of apoptotic and PI3K/Akt signaling pathways in the spinal cord may represent potential therapeutic strategy to treat I/R induced renal injury.
In conclusion, we identified and validated differentially expressed mRNA and lncRNA transcripts in the lower thoracic spinal cord by high-throughput RNA-seq that maybe involved in I/R induced acute renal injury. Further detailed studies of factors uncovered by this analysis may serve as novel diagnostic markers and therapeutics for I/R induced AKI in the future.

Animal care
All animal experiments were performed with adult male Sprague-Dawley (SD) rats (300 to 350g) in accordance with the guidelines according to the Huazhong University of Science and Technology Guide for the Care and Use of Laboratory Animals. All experimental animal procedures were approved by the Institutional Animal Care and Use Committee of Tongji Hospital. The rats were maintained and habituated in a standard 12h lightdark cycle with ad libitum access to food and drinking in a temperature-and humidity-controlled room (22 °C ± 0.5 °C, relative humidity 40%-60%).

Acute kidney ischemia/reperfusion experiments
Rats were randomly divided into sham (control, n=12) and I/R (AKI, n=12) groups. They were anesthetized with 40mg/kg sodium pentobarbital intraperitoneal injection and placed on a heating pad to maintain a constant temperature. Acute kidney ischemia/ reperfusion experiments were performed according to previously published protocols [59][60][61][62]. A midline abdominal incision was made and bilateral renal pedicles were exposed. Renal I/R were induced without trauma by clamping both renal pedicles for 45 minutes or 60 minutes. The ischemia was confirmed by visualizing dark color of kidneys. The blood flow was restored after clamp removal and the color of kidneys changed from dark to pink. To reduce abdominal air, 1ml warm normal saline was given intraperitoneally before abdominal closure. In sham controls, only bilateral renal pedicles were exposed. Rats were euthanized 24h after I/R, blood was drawn and their kidneys were thoroughly perfused with saline to remove any blood from the vascular beds. Kidneys and T8-12 spinal cord specimens were preserved at -80°C for further use.

Analysis of renal function and histology
Serum creatinine (SCr), serum BUN levels were determined in the Automated Blood or Urine Chemical Analyzer Vitro 350 (Orthoclinical Diagnostic Inc., Rochester, NY). For histology, kidney specimens were first fixed with 10% formalin solution and embedded with paraffin. Then, 3μm thick renal tissue sections were cut and stained with periodic acid-Schiff reagent (PAS). Histological analysis was performed by experienced pathologists blinded to the experimental groups. The kidney injuries were graded according to the semiquantitative scores developed by Paller et al [63,64] as follows: (1) tubular epithelial smoothness or tubular expansion: score 1, (2) loss of brush-like edge: score 1or 2, (3) obstruction of tubular lumen: score 1 or 2, (4) cytoplasmic vacuolization: score 1, and (5) cell necrosis: score 1. The highest achievable score was seven.

RNA isolation from spinal cord tissue
Total RNA from T8-12 segments of the spinal cord tissue from three rats was extracted for high-throughput sequencing using TRIzol® (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. RNA concentration and purity was determined with the NanoDrop 2000 Spectrophotometer (Thermo Scientific, Wilmington, DE) and agarose gel electrophoresis was used to assess RNA integrity. The A 260 /A 280 ratio was 1.94-2.2. for all samples. The RNA integrity numbers (RIN) were 9.5-9.8 as determined by RNA 6000 Nano kit in an Agilent 2100 Bionalyzer (Agilent technologies, Santa Clara, CA).

Validation of mRNAs and lncRNAs by quantitative real-time PCR (qRT-PCR)
Total RNA reverse transcription was conducted using PrimeScript™ RT reagent kit (TaKaRa) according to the manufacturer's instructions with either the oligo (dT) primers or specific RT primers. cDNA was quantitated by real-time PCR using the primers listed in Table 4 (Integrated DNA Technologies, Coralville, IA). Each sample was run in triplicate in 20μl reactions with 0.4 μM forward and reverse PCR primers and 10μl of Advanced Universal SYBR Green Supermix (Bio-Rad Laboratories, Hercules, CA) in a Bio-Rad CFX96 real-time PCR system. The PCR cycle parameters were set as follows: initial denaturation at 95°C for 1 min followed by 40 cycles of 95°C for 15s, 60°C for 15s, and 72°C for 45s. Relative expression was determined by normalization to GAPDH using the 2 −ΔΔCt method. The experiments were performed in triplicate.

Western blotting
Total protein extract was prepared from T8-12 spinal cord tissues with protein lysis buffer. Western blotting was performed as described previously [65,66]. 50μg of protein samples were separated by SDS-PAGE and transferred onto PVDF membranes (Millipore, USA). Then, the PVDF membrane was blocked with 5% skimmed milk in TBST for 1h at room temperature. Then, the membranes were incubated with the following primary antibodies: P2X7R (Boster, China), S100A9 (Boster, China), Akt (Cell Signaling Technology, USA), P-Akt (Cell Signaling Technology, USA), Bcl-2 (Affinity Biosciences, USA), Bax (Affinity Biosciences, USA), and Caspase-3 (ABclonal, USA). Further, after washes with 1X TBST, the membrane was incubated with anti-mouse or goat HRP-conjugated secondary antibodies (Abbkine, USA). This was followed by protein detection by ECL Assay Kit (Bipec Biopharma). β-actin (Abbkine, USA) was used as internal control. Band intensity was quantified using ImageJ software (1.44 P) and the expression of target proteins relative to β-actin control was determined.

Gene ontology (GO) and KEGG pathway analysis
The mRNA expression profiles from control and I/R group were screened by volcano plot filtering and differential mRNA expression was determined with Gene Ontology program (http://www.geneontology.org), which classified the data into biological processes, cellular components and molecular functions. Moreover, the key regulatory pathways in the spinal cord responding to I/Rinduced AKI were analyzed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (http:// www.genome.jp/kegg).

Statistical analysis
Data were presented as mean±SEM. Statistical analysis was performed using a t test between two cohorts. Data processing and subsequent quantitative normalization were performed using GeneSpring GX v12.1 software (Agilent Technologies). Differentially expressed mRNAs and lncRNAs between control group and I/R group were identified through fold change filtering. P value <0.05 was considered statistically significant. Hierarchical clustering and combined analysis was performed with homemade scripts (Multi-Experiment Viewer clustering).

Author contributions
QQL, HL, HBX and YML contributed substantially to the conception and design of the experiments. QQL, HL, ZGH, SJZ, BWL, WL, WHQ, QX and YML interpreted and analyzed results. BAD, QQL and HL drafted the initial version of the manuscript. HBX and YML were subsequently revised the manuscript.