Changes in long non-coding RNA transcriptomic profiles after ischemia-reperfusion injury in rat spinal cord

Spinal cord ischemia reperfusion injury model establishment

Male Sprague-Dawley rats (weight, 200–250 g) were purchased from the Animal center of the Second Military Medical University (Shanghai, China). All animal experiments were performed according to the guidelines of the Animal Ethics Committee of the Second Military Medical University (Shanghai, China). Thirty rats were randomly divided into the SCII group and sham group (nine in each). In SCII group, the SCII model was established as previously described (Li et al., 2014). After identification of the aortic arch branches, ischemia was achieved by cross-clamping the aortic arch for 14 min. Obstruction was confirmed by a 90% blood flow decrease in the tail artery using a Laser Doppler Bloodflow Imager (Perimed, Stockholm, Sweden). The clamps were released to reperfuse for 48 h. Rats in sham group that underwent the same procedure without occlusion were used as controls.

RNA extraction and high-throughput RNA sequencing

Spinal cord segments between L4 and L6 were harvested 48 h after reperfusion. In total, six samples (three from the SCII group and three from the sham group) were subjected for high-throughput RNA sequencing. For each sample, three random rats spinal cord tissues in the same group were mixed into one sample before RNA extraction. Total RNA was extracted using the mirVana miRNA Isolation Kit (Ambion, Texas, USA) following the manufacturer’s protocol. RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent, California, USA). The samples with RNA Integrity Number (RIN) ≥ 7 were subjected to the subsequent analysis. The libraries were constructed using TruSeq Stranded Total RNA with Ribo-Zero Gold (Illumina, California, USA) according to the manufacturer’s instructions. After validating the size and purity by Agilent Technologies 2100 Bioanalyze (Agilent, California, USA), these libraries were sequenced on the Illumina HiSeqTM 4000 sequencing platform (Illumina, California, USA) and 150 bp paired-end reads were generated. Libraries construction and RNA-sequencing were performed by OE Biotech (Shanghai, China).

Bioinformatics analysis

We obtained about 99 million raw RNA-sequencing reads from each sample. First, raw RNA-Seq data were flattened by removing adapter sequences and filtering out low quality reads using Trimmomatic software (Bolger, Lohse & Usadel, 2014). The clean reads were obtained by quality detection of the obtained reads through fastqc software (www.bioinformatics.babraham.ac.uk/projects/fastqc/), and then they were mapped to the rat reference genome (Rattus_norvegicus.Rnor_6.0) using hisat2 (Kim, Langmead & Salzberg, 2015). The transcript abundances of lncRNA and mRNA were quantified by normalized expression values as FPKM (Fragments Per kb Per Million Reads) using bowtie2 (Langmead & Salzberg, 2012) and eXpress (Roberts & Pachter, 2013). To identify differentially expressed lncRNAs or mRNAs between the spinal cord after ischemia reperfusion injury and the controls, differential expression analyses were performed using DESeq package. The p-value <0.05 and —fold change—≥ 2 were selected as the criteria for significantly DElncRNA or DEmRNA.

Quantitative real-time PCR (qRT-PCR)

To validate the reliability of the RNA sequencing results, We randomly selected and analyzed 8 DElncRNAs and 8 DEmRNAs by qRT-PCR analysis. L4–L6 spinal cord tissues were collected from the SCII and sham groups (n = 6 per group). Total RNA was extracted using the TRIzol reagent (Invitrogen, NY, USA) and reversed transcribed into cDNA using random primers. Quantification of lncRNA and mRNA was performed using an ABI PRISM® 7500 Sequence Detection System (Applied Biosystems, Forster, CA, USA). Relative lncRNA expression was normalized to the U6 expression and the expression of mRNA level was normalized to GAPDH using the 2-ΔΔCT method. Three independent experiments were conducted for each sample. The primer sequences are shown in Table 1.

Functional and pathway enrichment analyses

GO and KEGG pathway analyses were conducted to predict the potential functions of DEmRNAs and DElncRNAs based on their co-expressed DEmRNAs. GO analysis was performed to explore potential functions of genes and gene products from three aspects: biological process (BP), cellular component (CC), and molecular function (MF). KEGG pathway analysis was conducted to predict the involvement of the differentially expressed genes in the biological pathways. The top 10 enriched GO terms and top 20 enriched pathways were ranked by enrichment score (−log10 (p-value)) identified by the Database for Annotation, Visualization, and Integrated Discovery (DAVID; https://david.ncifcrf.gov/).

Co-expression of lncRNAs/mRNAs and functional prediction

The functions of lncRNAs were predicted by annotating the function of the co-expressed mRNAs. A co-expression network was built to explore the interaction between DElncRNAs and DEmRNAs in SCII pathogenesis. The Pearson’s correlation coefficient (PCC) was calculated between the expression levels of each DElncRNA-DEmRNA. We then selected co-expression DElncRNA-DEmRNA pairs with value of —PCC— ≥0.80 and p-value <0.05 for network construction (Alkan et al., 2017).

Genomic localizations of the paired lncRNAs and mRNAs were identified for cis prediction. The co-expression nearby gene, which is less than 100 kb upstream or downstream from the lncRNA, can act as the potential target regulated by the lncRNA in a cis manner, while a trans-regulator is one that does not meet this criterion. The RIsearch-2.0 software (Alkan et al., 2017) was used to identify target genes in trans, with the parameter set as the base number of direct interactions between lncRNA and mRNA ≥10 and free energy ≤−50.

High throughput RNA sequencing and genome-wide read mapping

L4 to L6 spinal cord segments were harvested for total RNA isolation and Illumina TruSeq library construction (Data available at GEO: GSE138966). 95.49 M, 96.94 M and 96.77 M high quality sequence reads were obtained from three samples of the SCII group, and 95.94 M, 95.22 M and 96.43 M clean reads were generated from three control tissues of the sham-operated group (Table 2). All clean reads were mapped to the rat reference genome using hisat2. Nearly 98% of the reads mapped to the rat genome; 75.25%–81.13% of the reads mapped to unique genomic regions among the aligned fragments (Table 2), indicating the reliability of the data. In total, 25,328 lncRNAs and 55,760 mRNAs were identified.

The expression profile of lncRNA and mRNA in the SCII model

In this study, 1,455 lncRNAs and 6,707 mRNAs were significantly differentially expressed, with —fold change—≥2.0, P < 0.05 and FDR <0.05. In total, there were 761 up-regulated lncRNAs, 694 down-regulated lncRNAs, 3,772 up-regulated mRNAs, and 2,935 down-regulated mRNAs. Scatter plots analyses showed the expression signatures (Figs. 1A and 1C). Hierarchical clustering expression showed significant differences in the spinal cord between SCII rats and control rats (Figs. 1B and 1D). Among the top 10 DElncRNAs, NONRATT007222.2 and NONRATT019127.2 were the most significantly up-regulated and down-regulated in the spinal cord after SCII (Table 3). Moreover, Csf3 and Ankrd9 were the most significantly up- and down-regulated mRNAs in spinal cord after SCII (Table 4). Further, transcripts were proved to be distributed on all chromosomes (Fig. 2).

Validation of lncRNAs and mRNAs

Eight mRNAs and 8 lncRNAs were randomly selected to verify the high throughput RNA sequencing results in six sample pairs by qRT-PCR. We found that the expressions of lncRNAs NONRATT007222.2, NONRATT013040.2, TCONS_00018593 and TCONS_00008647 were up-regulated, while NONRATT019127.2, NONRATT013069.2, TCONS_00018687 and TCONS_00020362 were down-regulated. Of the eight selected mRNAs, NM_017104.2, NM_001191577.1, XM_017592198.1 and NM_013037.1 were up-regulated, and NM_001142366.2, NM_001191578.1, XM_017592225.1 and NM_001191043.1 were down-regulated in the SCII model (Fig. 3). The qRT-PCR results were in line with the high throughput RNA sequencing. Therefore, the results validated our high throughput RNA sequencing profile as highly reliable and showed that these lncRNAs and mRNAs might be involved in SCII pathogenesis.

Gene ontology and KEGG Pathway annotation for differentially expressed genes

In GO and KEGG pathway enrichment analyses of DEmRNAs, we found 6707 mRNAs that were differentially expressed. GO enrichment analysis showed that the enriched GO biological processes for up-regulated genes in the SCII group were negative regulation of CD40 signaling pathway, positive regulation of T-helper 17 cell differentiation, cellular response to glucoside, cellular response to diterpene and cellular response to cyanide. The enriched GO cellular components for up-regulated genes in the SCII group were for the CRLF-CLCF1 complex, inner dense plaque of desmosome, outer dense plaque of desmosome, integrin alphaL-beta2 complex and integrin alpha1-beta1 complex. The enriched GO molecular functions for up-regulated genes in the SCII group were chemokine (C-C motif) ligand 5 binding and NEDD8 ligase activity (Fig. 4A). Down-regulated mRNAs were enriched in dendrite arborization, maintenance of synapse structure, neuron projection, axon, neuronal cell body, extracellularly glycine-gated ion channel activity, protein binding and calcium ion binding of GO biological process, cellular component and molecular function (Fig. 4B).

Similarly, differentially expressed genes were analyzed using KEGG. We found that up-regulated genes in the spinal cord after SCII were involved in cytokine-cytokine receptor interaction (rno04060), TNF (rno04668) and NF-κB signaling pathway (rno04064) (Fig. 5A), while down-regulated genes were involved in axon guidance (rno04360), calcium (rno04020) and cAMP signaling pathway (rno04024) (Fig. 5B).

lncRNAs/mRNAs co-expression analysis and functional prediction

Co-expression analysis was performed to understand the roles of lncRNAs according to their co-expressed mRNAs. 385 DEmRNAs and 221 DElncRNAs, which were composed of 606 nodes and 500 edges (Fig. S1) were involved. Additionally, positively co-expressed and negatively co-expressed networks were constructed to show two distinct co-expression patterns of lncRNA-mRNA associated with SCII (Figs. 6A and 6B). Among these co-expressed mRNAs, the top three up-regulated DEmRNAs were Ermap, Hmox1and Il1r1, and these mRNAs were involved in regulation of cytokine production, apoptotic process, wound healing involved in inflammatory response, regulation of inflammatory response and positive regulation of neuron migration. The down-regulated DEmRNAs were Tmem255b, Smarcal1and Cpne5. These genes are associated with many functions, such as cellular response to DNA damage stimulus, neuron projection and regulation of dendrite extension. To better predict the functions of DElncRNAs identified in our study, functional annotation of their co-expressed mRNAs were performed (Figs. 7A and 7B). The KEGG pathway analysis showed that the most significantly enriched pathway was the TNF signaling pathway and that, there are 62 mRNAs that might be associated with SCII among these co-expressed mRNAs.

Prediction of cis-and trans-regulated genes of lncRNAs

Cis- and trans-regulated target genes of the DElncRNAs were predicted to further explore how lncRNAs regulate the pathological process of SCII. Thirty-four lncRNAs were found to have at least five cis-regulated target genes, some of which have been reported to be associated with SCII (Fig. 8A). The top five lncRNAs ranked by the number of predicted nearby coding genes were TCONS_00026300, TCONS_00026304, TCONS_00026343, TCONS_00026342 and TCONS_00026425. LncRNAs and their trans-regulated genes are displayed in Fig. 8B. Each lncRNA has a different number of potential coding genes. For example, NONRATT007903.2 had a maximum of 12 trans-regulated genes, but NONRATT028952.2 and NONRATT001746.2 only had 1 trans-regulated gene.