Genome-wide analysis of acute traumatic spinal cord injury-related RNA expression profiles and uncovering of a regulatory axis in spinal fibrotic scars

Abstract

Objectives: Long non-coding RNAs (lncRNAs) are critical for posttranscriptional and transcriptional regulation in eukaryotic cells. However, data on lncRNA expression in the lesion epicentres of spinal tissues after acute traumatic spinal cord injury (ATSCI) are scarce. We aimed to identify lncRNA expression profiles in such centres and predict latent regulatory networks.

Materials and methods: High-throughput RNA-sequencing was used to profile the expression and regulatory patterns of lncRNAs, microRNAs and messenger RNAs (mRNAs) in an ATSCI C57BL/6 mouse model. Chromosome distributions, open reading frames (ORFs), transcript abundances, exon numbers and lengths were compared between lncRNAs and mRNAs. Gene ontology, KEGG pathways and binding networks were analysed. The findings were validated by qRT-PCRs and luciferase assays.

Results: Intronic lncRNAs were the most common differentially expressed lncRNA. Most lncRNAs had <6 exons, and lncRNAs had shorter lengths and lesser ORFs than mRNAs. MiR-21a-5p had the most significant differential expression and bound to the differentially expressed lncRNA ENSMUST00000195880. The microRNAs and lncRNAs with significant differential expression were screened, and a lncRNA/miRNA/mRNA interaction network was predicted, constructed and verified.

Conclusions: The regulatory actions of this network may play a role in the pathophysiology of ATSCI. Our findings may lead to better understanding of potential ncRNA biomarkers and confer better therapeutic strategies for ATSCIs.

Keywords: acute traumatic spinal cord injury; interaction network; long non-coding RNAs; mRNA; miRNA.

FIGURE 1

Construction of mouse SCI model and a flowchart of the identification of differently expressed RNAs HE staining of spinal cord samples. Sham group (A) and ATSCI (B) group data are shown. Overview of the analysis pipeline(C). HE, haematoxylin‐eosin; ATSCI, acute traumatic spinal cord injury; lncRNA, long noncoding RNA; mRNA, messenger RNA; miRNA, microRNA

FIGURE 2

Long noncoding, micro and messenger RNA expression profiles LncRNAs are shown in the right‐hand column. Red pixels correspond to an increased abundance of the gene in the indicated sample whereas blue pixels indicate decreased levels; all differentially expressed lncRNAs show a fold change >1.5 and P < .05(A). Volcano plots elucidate the variance in differentially expressed lncRNA(B), mRNAs (C) and miRNAs (D) based on P‐values and fold changes. The x‐axis is the fold change (log 2), and the y‐axis is the P‐value (−log 10). Red points (fold change >2) indicate upregulated mRNAs or lncRNAs and blue points (fold change < −2) indicate downregulation in volcano plots (B, C). Red points in the volcano plot (D) indicate a fold change >2 or <−2. lncRNA, long noncoding RNA; mRNA, messenger RNA; miRNA, microRNA

FIGURE 3

Chromosome distribution and comparison between long noncoding RNAs and messenger RNAs Distribution of 6 types of lncRNAs along each chromosome. Known lncRNAs (class code =, depicted in orange), intronic lncRNAs (class code i, depicted in light green), lncRNAs that share a reference with at least 1 splice junction (class code j, depicted in dark green), lncRNA of generic exonic overlap with a reference transcript (class code o, depicted in blue), intergenic lncRNA (class code u, depicted in violet) and antisense lncRNA (class code x, depicted in pink) are presented in physical bins of 500 kb for each chromosome (A). The ORF lengths of lncRNAs and mRNAs (B, C) are shown. The transcript length of lncRNAs and mRNAs is shown (D). Exon numbers of lncRNAs and mRNAs are shown (E). Expression levels of lncRNAs and mRNAs (F) are shown. lncRNA, long noncoding RNA; mRNA, messenger RNA

FIGURE 4

Enriched GO terms and KEGG pathways of differentially expressed lncRNAs Enriched GO terms (A) and enriched KEGG pathways (B) are shown. The rich factor is the ratio of the number of different genes to the total number of genes in the KEGG database; the higher the rich factor value, the greater the enrichment degree. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; lncRNA, long noncoding RNA (C) LncRNA/miRNA/mRNA regulatory interaction network of TargetScan and MiRanda prediction. lncRNA, long noncoding RNA; mRNA, messenger RNA; miRNA, microRNA

FIGURE 5

Validation of differential ncRNAs and confirmation of their relationship Expression of ENSMUST00000195880 and miR‐21a‐5p in the spinal tissues validated by qRT‐PCR (A, B). A predicted miR‐21a‐5p target site in ENSMUST00000195880 using TargetScan analysis (C) is shown. The construction of ENSMUST00000195880 (WT) and ENSMUST00000195880 (MT) luciferase plasmids (D, E) is shown. Relative luciferase expression of WT and MT ENSMUST00000195880 vectors co‐transfected with miR‐21A‐5p vectors (F) is shown. Expression of Smad7 protein in the spinal tissues validated by Western blot (G); a predicted miR‐21a‐5p target site in ENSMUST00000195880 using TargetScan analysis (H); a predicted ENSMUST00000195880/miR‐21a‐5p/Smad7 target interaction axis after SCI (I). qRT‐PCR, quantitative real‐time polymerase chain reaction, WT, wild‐type; MT, mutant; ***P < .001

FIGURE 6

Ethological and histological analysis of miR‐21a‐5p knockout mice BMS scores indicate the motor functional index over 10 d after SCI (A). Immunohistochemistry to determine the expression of fibronectin (B, C, D) and Smad7 (E, F, G). Fibronectin and Smad7 levels in the lesion epicentre of spinal cords by immunohistochemistry (H). ***P < .001

FIGURE 7

The validation of regulatory network in vitro The identification of primary spinal cord fibroblasts (A, B). the predominant localization of lncENSMUST00000195880 in primary spinal cord fibroblasts (C). Expression of miR‐21a‐5p validated by qRT‐PCR (D) and Expression of fibronectin, collagen I, Smad7 and Smad2/3 phosphorylation‐related proteins by Western blot (E, F) after lncENSMUST00000195880 overexpression. ***P < .001

FIGURE 8

Working model of the target interaction axis after SCI After SCI, TGF‐β1 expression was increased, leading to phosphorylation of the Smad2/3 protein, thereby promoting miR‐21a‐5p upregulation. Smad7 is an inhibitor of the TGF‐β pathway that can be inhibited by miR‐21a‐5p. miR‐21a‐5p overexpression alleviates Smad7 inhibition on the TGF‐β pathway. However, the upregulated lncENSMUST00000195880, which binds with miR‐21a‐5p, suppresses this positive feedback loop