Research Paper Volume 11, Issue 8 pp 2352—2368

Identification of noncoding RNA expression profiles and regulatory interaction networks following traumatic spinal cord injury by sequence analysis

Wenzhao Wang1,2, , Yanlin Su1, , Shi Tang3, , Hongfei Li1, , Wei Xie4, , Jianan Chen1, , Lin Shen1, , Xinda Pan1, , Bin Ning1, ,

  • 1 Jinan Central Hospital Affiliated to Shandong University, Jinan, China
  • 2 Department of Orthopaedics, West China Hospital, Sichuan University, Chengdu, China
  • 3 Department of Radiology, West China Hospital, Sichuan University, Chengdu, China
  • 4 Department Emergency Medicine, Affiliated Hospital of Taishan Medical University, Taian, China

Received: March 11, 2019       Accepted: April 10, 2019       Published: April 18, 2019      

https://doi.org/10.18632/aging.101919
How to Cite

Copyright: Wang et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Aim: To systematically profile and characterize the noncoding RNA (ncRNA) expression pattern in the lesion epicenter of spinal tissues after traumatic spinal cord injury (TSCI) and predicted the structure and potential functions of the regulatory networks associated with these differentially expressed ncRNAs and mRNAs.

Results: A total of 498 circRNAs, 458 lncRNAs, 155 miRNAs and 1203 mRNAs were identified in TSCI mice models to be differentially expressed. The regulatory networks associated with these differentially expressed ncRNAs and mRNAs were constructed.

Materials and methods: We used RNA-Seq, Gene ontology (GO), KEGG pathway analysis and co-expression network analyses to profle the expression and regulation patterns of noncoding RNAs and mRNAs of mice models after TSCI. The findings were validated by quantitative real-time PCR (qRT-PCR) and Luciferase assay.

Conclusion: noncoding RNAs might play important roles via the competing endogenous RNA regulation pattern after TSCI, further findings arising from this study will not only expand the understanding of potential ncRNA biomarkers but also help guide therapeutic strategies for TSCI.

Introduction

Spinal cord injury (SCI) is a disabling neurological condition with high economic and social costs; SCI is characterized by the loss of neural tissue and consequent deficits in sensory and motor functions [1]. Each year, half a million people damage their spinal cord, and the injury is almost always life-changing [2]. SCI increases the risk of involuntary movements, bladder and gastrointestinal disorders, and depression [1]. SCI can also be caused by iatrogenic procedures, infection, vascular lesions or tumors, but the most common cause is trauma [3]. Traumatic SCI (TSCI) causes cell necrosis, the disconnection of surviving neurons, and the irreversible interruption of ascending and descending neurotransmission [4]. Unfortunately, recent studies have demonstrated that no effective treatments exist for achieving complete neurological or functional recovery after TSCI. Moreover, the key mechanisms governing the cellular response to injury are largely unknown [5]. A better understanding of the cellular and molecular mechanisms following TSCI is necessary to develop new strategies to promote axonal regeneration and functional recovery.

Noncoding ribose nucleic acids (ncRNAs), which are a class of genetic, epigenetic and translational regulators, have been found to play key roles in various physiological and pathological processes [6]. No less than 70% of the human genome is transcribed, but protein-coding transcripts account for no more than 2%, and extensive transcripts derived from most of the genome generate a large proportion of ncRNAs [7]. Theoretically, ncRNAs do not encode proteins but instead functionally regulate the translation of proteins and can be classified into two types: housekeeping ncRNAs, which consist of small nucleolar RNAs (snoRNAs), small nuclear RNAs (snRNAs), rRNAs, and tRNAs; and regulatory ncRNAs, which consist of microRNAs (miRNAs), long ncRNAs (lncRNAs) with a relatively flexible length of >200 nucleotides, and circular RNAs (circRNAs) with a closed-loop structure [8]. To date, miRNAs are the most extensively studied class of small noncoding RNAs (sncRNAs) [9]; miRNAs are present in a wide range of tissues and fluids [10, 11] and play an essential role in neurological and neurodegenerative diseases by regulating cell-to-cell communication as hormone-like molecules to influence the behaviors of different cells in a paracrine or endocrine manner [12]. Compared to sncRNAs, lncRNAs are more heterogeneous in size, often polyadenylated, longer and lack open reading frames (ORFs). In early studies, the importance of lncRNAs was vastly underestimated because of their low levels of sequence conservation and expression [13]. However, accumulating evidence indicates that lncRNAs play essential roles in the development of diseases in various organisms [14]. In addition, circRNA has recently been identified as a novel type of endogenous ncRNA that is abundant yet enigmatic in mammalian cells. Unlike linear RNAs that are terminated with 5′ caps and 3′ tails, circRNAs are characterized by a covalent closed-loop structure formed by a back-splicing event, without 5′ caps or poly-A tails. Notably, one of the most frequently studied functions of circRNA is the miRNA sponge [15]. Therefore, ncRNAs have potential as candidate diagnostic biomarkers and therapeutic targets in patients with TSCI.

The regulating functions of ncRNAs after TSCI and their underlying functional mechanisms have not yet been sufficiently and systematically described. Therefore, extensive prediction and analysis of the ncRNAs regulating the progression of TSCI is fundamental for the development of understanding the underlying mechanisms and finding effective therapeutic strategies. Our study analyzed the profiles and predicted the function of differentially expressed (DE) ncRNAs in the epicenter of spinal cord lesions in a modified Allen’s weight-drop model using RNA sequencing techniques to provide a better comprehending of the diagnostic, prognostic and therapeutic value of ncRNAs.

Results

DE ncRNAs and mRNAs

To identify the effect of TSCI on ncRNA expression in the lesion epicenter, we applied a standard Allen’s weight-drop model. SCI mice started to show improvements in locomotor function 2 days after TSCI, but during the first two days, the BMS was rated zero. The BMS of mice in the sham group showed an improvement on day 1 and returned to normal on day 3 postsurgery (Figure 1A). The spinal tissues of inbred C57 mice damaged by Allen’s impactor were sliced and stained with H&E. Staining results demonstrated severe damage to the blood-spinal cord barrier and the structural integrity of the lesion epicenter, including rupture, hemorrhage and inflammatory cell infiltration (Figure 1B1D). To further determine whether ncRNAs were involved in the TSCI, the total RNA of the lesion epicenter at the T8–10 level was analyzed by RNA sequencing techniques. P-values <0.05 were used to assess the normalized expression of genes. The dysregulated ncRNAs and mRNAs are shown in a table, cluster map, volcano plot and Venn diagram. Information on the top 40 dysregulated circRNAs, lncRNAs, miRNAs, and mRNAs is listed in order of ascending p-value (Tables 14). The cluster map, volcano plot and Venn diagram of DE circRNAs, lncRNAs, miRNAs, and mRNAs after TSCI are shown in Figure 2. According to the data, we summarized the dysregulated RNAs in the TSCI samples compared with those in the sham samples, as follows: 249 circRNAs were upregulated, and 249 circRNAs were downregulated; 356 lncRNAs were upregulated, and 93 lncRNAs were downregulated; 94 miRNAs were upregulated, and 61 miRNAs were downregulated; 1098 mRNAs were upregulated, and 105 mRNAs were downregulated (Figure 3A). DE ncRNAs could directly or indirectly target genes and regulate the expression of target mRNAs. The results of an intersectional analysis of DE circRNAs, lncRNAs, and miRNAs and their target DE mRNAs are shown in a Venn diagram (Figure 3B3D).

Establishment of SCI animal model. (A) BMS scores indicate the motor functional index 3 days after SCI. ***PB–D) H&E staining of spinal cord samples from the sham and SCI groups at days 1 and 3 postsurgery.

Figure 1. Establishment of SCI animal model. (A) BMS scores indicate the motor functional index 3 days after SCI. ***P<0.001. (BD) H&E staining of spinal cord samples from the sham and SCI groups at days 1 and 3 postsurgery.

Table 1. Top 40 differently expressed miRNAs in SCI tissues comparing with Sham tissues.

miRNAGenome IDStrandp-valuelog2_FoldChangeRegulation
mmu-miR-344e-3pchr70.0005751.33up
mmu-miR-106b-3p_R-2chr50.0015981.12up
mmu-miR-5099_L+2R-1chr12+0.0020162.03up
mmu-miR-15b-3pchr3+0.0023720.91up
mmu-miR-7688-5pchr10+0.0032692.01up
mmu-miR-1964-3pchr7+0.0045161.73up
mmu-miR-130b-3pchr160.0046661.65up
mmu-miR-155-5pchr160.0054211.96up
mmu-miR-27a-5pchr8+0.0059292.30up
mmu-miR-18a-3pchr14+0.0064403.07up
mmu-miR-18a-5pchr14+0.0072642.03up
mmu-miR-223-3p_R+1chrX+0.0072872.52up
mmu-miR-214-3pchr1+0.0081002.37up
mmu-miR-92a-1-5pchr14+0.0098573.16up
mmu-miR-28a-3pchr16+0.0106360.98up
mmu-miR-877-5p_R+4chr170.0113480.90up
mmu-miR-21a-5p_R+1chr110.0121982.46up
mmu-miR-144-3p_R-1chr11+0.0136320.87up
mmu-miR-222-3p_R+2chrX0.0137010.91up
mmu-miR-511-3pchr2+0.0137591.02up
mmu-miR-369-3pchr12+0.001349−0.68down
mmu-miR-384-3pchrX0.001676−0.87down
mmu-miR-325-5p_R-2chrX0.002973−0.49down
mmu-miR-34a-5pchr4+0.004289−0.84down
mmu-miR-383-5pchr80.005397−0.51down
mmu-miR-128-3pchr1+0.005516−0.75down
mmu-miR-30e-5p_R+2chr40.006287−0.67down
mmu-miR-411-3p_R-1chr12+0.006918−0.52down
mmu-miR-329-5p_R+2chr12+0.007447−0.83down
mmu-miR-1298-3pchrX+0.009603−1.07down
mmu-miR-1264-5pchrX+0.009893−1.20down
mmu-miR-135a-5pchr9+0.010159−0.92down
mmu-miR-218-5pchr5+0.011109−0.90down
mmu-miR-6516-5p_R+3chr11+0.011772−0.86down
mmu-miR-488-3pchr1+0.013430−0.85down
mmu-miR-7b-5p_R+1chr17+0.013986−1.00down
mmu-miR-1843a-3pchr120.014024−0.68down
mmu-miR-3069-3pchr120.015073−1.37down
mmu-miR-582-5pchr130.015134−0.90down
mmu-miR-204-5pchr19+0.015457−0.77down

Table 2. Top 40 differently expressed CircRNAs in SCI tissues comparing with Sham tissues.

CircRNA IDChromGene NameStrandp-valuelog2_FoldChangeRegulation
circRNA8075chr14Diaph30.0002883.38up
circRNA81chr2Thbs1+0.0004594.71up
circRNA172chr11Col1a10.0006184.71up
circRNA2355chr4Kif2c0.0006943.73up
circRNA9086chr5Antxr20.0008912.02up
circRNA9357chr4Pgd0.0009302.73up
circRNA7783chr15Racgap10.0011372.42up
circRNA5480chr13Nln0.0011741.10up
circRNA8879chr6Zc3hav10.0013082.09up
circRNA3369chr2Atp8b40.0013672.88up
circRNA13618chr11Psmd3+0.0017361.06up
circRNA6169chr11Ankfy1+0.0020371.66up
circRNA8690chr8Gm20388+0.0021832.11up
circRNA8894chr6Skap20.0023602.31up
circRNA7312chr16Pak20.0024201.62up
circRNA8074chr14Diaph30.0025073.69up
circRNA14292chr5Gsap+0.0026211.83up
circRNA7782chr15Racgap10.0026232.33up
circRNA8712chr7Lig1+0.0029842.66up
circRNA4569chr9Fli1+0.0032971.98up
circRNA2098chr4Anp32b+0.0033531.43up
circRNA15118chr1Dpp100.000142−1.36down
circRNA13219chr1Pld5+0.002372−1.41down
circRNA4130chr1Pld50.000752−1.51down
circRNA11130chr15Lrrc60.000863−2.00down
circRNA15238chr18Nol40.000946−1.20down
circRNA4505chr9Cntn50.001476−1.46down
circRNA4277chr18Asxl3+0.001539−1.02down
circRNA1084chr6Dync1i1+0.001770−1.13down
circRNA2986chr2Cacna1b0.002072−1.62down
circRNA402chr9Myrip+0.002090−1.54down
circRNA7067chr17L3mbtl4+0.002227−1.83down
circRNA2737chr3Pogz+0.002259−1.76down
circRNA4247chr18Greb1l+0.002611−1.69down
circRNA14575chr4Rimkla0.002659−1.64down
circRNA8251chr19Cpeb30.002934−1.13down
circRNA6448chr11Tbc1d160.003047−1.20down
circRNA16650chr14Sfmbt1+0.003097−1.12down
circRNA12235chr12Dtnb+0.003164−1.10down
circRNA16513chr14Sfmbt1+0.003400−1.14down

Table 3. Top 40 differently expressed lncRNAs in SCI tissues comparing with Sham tissues.

Gene idGene nameStatusp-valuelog2_FoldChangeRegulation
MSTRG.34122Homeznovel0.0000051.98up
MSTRG.67284Tpd52novel0.0000331.30up
MSTRG.111948Nanovel0.0000462.66up
MSTRG.74174Gm15689novel0.0000541.26up
MSTRG.124087F630028O10Rikknown0.0000851.46up
MSTRG.40161Rbfox2novel0.0000951.55up
MSTRG.113073Cfap20novel0.0001091.64up
MSTRG.95251Gm44170known0.0001091.80up
MSTRG.93249Nanovel0.0001592.29up
MSTRG.55738BE692007known0.0001892.64up
MSTRG.52721Dpysl3novel0.0002251.45up
MSTRG.60157Arhgap15novel0.0002691.33up
MSTRG.33815Nanovel0.0002841.62up
MSTRG.25765Nanovel0.0003092.17up
MSTRG.77436Gm11216known0.0003101.53up
MSTRG.107512Fgfr2novel0.0003912.08up
MSTRG.94129Nanovel0.0003951.88up
MSTRG.11903Pcbp3novel0.0003991.55up
MSTRG.123405Mamld1novel0.0004211.34up
MSTRG.105294Nanovel0.0005091.61up
MSTRG.35205Gm26908novel0.0006091.25up
MSTRG.14082Ppm1hnovel0.0006251.24up
MSTRG.33653NAnovel0.0009542.29up
MSTRG.33653Thumpd2novel0.0009871.18up
MSTRG.83712AI506816known0.0010012.80up
MSTRG.10902NAnovel0.0010141.60up
MSTRG.63860Gm14005novel0.0010421.50up
MSTRG.110801Ddx60novel0.0010811.47up
MSTRG.75978NAnovel0.0010842.21up
MSTRG.99804Kcnn4novel0.0011401.92up
MSTRG.108182Tacc1novel0.0011431.26up
MSTRG.28498Cenppnovel0.0012031.12up
MSTRG.47244NAnovel0.0012551.63up
MSTRG.24916NAnovel0.0012602.09up
MSTRG.125937Diaph2novel0.0012631.18up
MSTRG.987013300002P13Rikknown0.000457−1.34down
MSTRG.52759Pcdha6novel0.000674−1.45down
MSTRG.83528Pclonovel0.000711−1.02down
MSTRG.24945NAnovel0.000748−1.47down
MSTRG.11329Rhobtb1novel0.001009−1.41down
NA, not annotated.

Table 4. Top 40 differently expressed mRNAs in SCI tissues comparing with Sham tissues.

Gene idGene namep-valuelog2_FoldChangeRegulation
MSTRG.101002Snord340.0000131.51up
MSTRG.80254Gjb30.0000251.58up
MSTRG.104414P2ry60.0000702.41up
MSTRG.114824Mmp30.0000742.00up
MSTRG.65092H130.0000741.03up
MSTRG.27838F13a10.0000762.01up
MSTRG.13769Lyz20.0000864.44up
MSTRG.121761Ccr10.0000873.38up
MSTRG.20499Naglu0.0001292.95up
MSTRG.3656Ugt1a10.0001541.21up
MSTRG.107440F70.0001771.60up
MSTRG.105998Nupr10.0001852.18up
MSTRG.12171Sbno20.0001951.97up
MSTRG.111467Bst20.0001962.64up
MSTRG.19369Ccl70.0002123.65up
MSTRG.28823Rgs140.0002211.01up
MSTRG.108769Msr10.0002414.58up
MSTRG.117687Tagln0.0002414.68up
MSTRG.52513Cd140.0002433.52up
MSTRG.86784Cxcl10.0002442.26up
MSTRG.92027Flnc0.0002481.89up
MSTRG.6767Selp0.0002661.50up
MSTRG.31848Nid20.0002721.27up
MSTRG.19493Ccl40.0002751.37up
MSTRG.104483Folr20.0002941.84up
MSTRG.21200Cd300ld0.0003251.72up
MSTRG.104929Adm0.0003301.32up
MSTRG.1578Col5a20.0003471.87up
MSTRG.100407Fxyd30.0003522.38up
MSTRG.33361Wdfy40.0003871.01up
MSTRG.89616Lat20.0004353.08up
MSTRG.114822Mmp120.0004511.10up
MSTRG.34965Pbk0.0004582.92up
MSTRG.97389Ptpn60.0004712.19up
MSTRG.127268Tmsb4x0.0004882.82up
MSTRG.20948Milr10.0005012.23up
MSTRG.40133Ncf40.0005082.51up
MSTRG.105902Il4ra0.0005092.79up
MSTRG.78174Elavl20.000260-1.28down
Expression profiles of DE ncRNAs and mRNAs in the lesion epicenter after SCI. (A) Heat map of DE lncRNAs in the SCI group compared with the sham group. (B) Heat map of DE circRNAs. (C) Volcano plot indicating the differential expression of lncRNAs. (D) Volcano plot of circRNAs. (E) Heat map of DE miRNAs. (G) Volcano plot of miRNAs. (F) Heat map of DE mRNAs. (H) Volcano plot of mRNAs. Up-regulated and down-regulated genes are colored in red and blue, respectively.

Figure 2. Expression profiles of DE ncRNAs and mRNAs in the lesion epicenter after SCI. (A) Heat map of DE lncRNAs in the SCI group compared with the sham group. (B) Heat map of DE circRNAs. (C) Volcano plot indicating the differential expression of lncRNAs. (D) Volcano plot of circRNAs. (E) Heat map of DE miRNAs. (G) Volcano plot of miRNAs. (F) Heat map of DE mRNAs. (H) Volcano plot of mRNAs. Up-regulated and down-regulated genes are colored in red and blue, respectively.

Overview of relative differential expression of ncRNAs. (A) Histogram showing the number of dysregulated ncRNAs and mRNAs. (B–D) Venn diagram showing the overlap between the target mRNAs of dysregulated ncRNAs and dysregulated mRNAs.

Figure 3. Overview of relative differential expression of ncRNAs. (A) Histogram showing the number of dysregulated ncRNAs and mRNAs. (BD) Venn diagram showing the overlap between the target mRNAs of dysregulated ncRNAs and dysregulated mRNAs.

Validation of ncRNA and mRNA expression

To validate the reliability of the sequencing data, the changes in the expression of 12 DE ncRNAs and mRNAs, including three circRNAs (circ2464, circ7435, circ7010), three lncRNAs (Gm12840, Gm26809, H19), three miRNAs (miR-21a-5p-R+1, miR-92a-3p-R+1, miR-423-3p) and three mRNAs (Ftl1, Lyz2, Tmsb4x) in the lesion epicenter compared with the sham group were were randomly selected for qRT-PCR analysis (Figure 4B, 4D). All the validated qRT-PCR results of the DE ncRNAs and mRNAs were consistent with the corresponding sequencing data (Figure 4A, 4C).

Validation of differential ncRNA and mRNA expression. (A, C) Sequencing results of the ncRNAs and mRNAs. (B, D) Expression of corresponding ncRNAs and mRNAs validated by qRT-PCR.

Figure 4. Validation of differential ncRNA and mRNA expression. (A, C) Sequencing results of the ncRNAs and mRNAs. (B, D) Expression of corresponding ncRNAs and mRNAs validated by qRT-PCR.

Enrichment of biological functions and pathway networks

To detect the enrichment categories and to examine the underlying functions of ncRNAs DE after TSCI, DE ncRNAs and DE mRNAs were subjected to GO and KEGG pathway analyses. DE mRNAs and coexpressed or target mRNAs of DE ncRNAs were identified. The GO molecular function analysis showed that the dysregulated transcripts of ncRNAs were associated with cell division, focal adhesion, proteinaceous extracellular matrix, positive regulation of cell migration, extracellular matrix components, regulation of cell shape, integrin binding, defense response to bacterium and leukocyte cell-cell adhesion (Figure 5A). In addition, the significant GO items indicated that mRNAs DE after TSCI were significantly associated with cytoplasm, protein binding, extracellular exosome, extracellular space, cell surface, focal adhesion, innate immune response and mitotic nuclear division (Figure 5B). Correspondingly, the top 20 ncRNA-associated pathways were demonstrated by KEGG analysis, and the most significantly associated pathways were cytokine-cytokine receptor interaction, cell cycle, leukocyte transendothelial migration, phagosome, Leishmaniasis, Malaria and Systemic lupus erythematosus pathways, (Figure 5C). In addition, KEGG pathway analysis of DE mRNAs revealed significant associations with cytokine-cytokine receptor interaction, focal adhesion, phagosome, chemokine signaling, Regulation of actin cytoskeleton, lysosome, Cell cycle and Toll-like receptor signaling pathways, among others (Figure 5D).

Enriched GO terms and KEGG pathways of host genes of DE ncRNAs in SCI mice. (A) Top 20 significantly enriched GO terms of DE ncRNAs are shown in the scatterplot. (B) Top 20 significantly enriched GO terms of DE mRNAs. (C) The top 20 significantly enriched KEGG pathways of DE ncRNAs are shown in the scatterplot. (D) The top 20 significantly enriched KEGG pathways of DE mRNAs are listed.

Figure 5. Enriched GO terms and KEGG pathways of host genes of DE ncRNAs in SCI mice. (A) Top 20 significantly enriched GO terms of DE ncRNAs are shown in the scatterplot. (B) Top 20 significantly enriched GO terms of DE mRNAs. (C) The top 20 significantly enriched KEGG pathways of DE ncRNAs are shown in the scatterplot. (D) The top 20 significantly enriched KEGG pathways of DE mRNAs are listed.

Regulatory networks of ncRNAs and mRNAs

The network of interactions of the host genes of these DE ncRNAs was also examined to elucidate the molecular mechanisms underlying the pathogenesis of TSCI. Considering that an important biological function of competing endogenous RNAs (ceRNAs) is binding to miRNAs, the binding relationships between ceRNAs and miRNAs were preliminarily determined. The miRNA- binding sites of lncRNAs and circRNAs were identified to construct lncRNA/circRNA–miRNA–mRNA coexpression networks. miRNA was used as the center of each network, which clearly shows possible regulated target genes. The lncRNA/circRNA–miRNA–mRNA interaction networks were conveniently displayed using Cytoscape. Eight DE miRNAs, i.e., miR-23a-5p, miR-222-3p, miR-223-3p, miR-22-5p, miR-218-5p, miR-214-5p, miR-21a-3p, and miR-21a-5p, and their paired ceRNAs and mRNAs were selected as intuitive examples showing parts of the whole complicated network involved in the pathogenesis of TSCI (Figures 67). The results show the regulatory relationship between ncRNAs and mRNAs with regard to TSCI. Furthermore, two pairs of binding relationships between ceRNAs and miRNAs were verified with a dual-luciferase reporter system. We found that the overexpression of miR-21-5p significantly decreased the luciferase activity of reporter vectors containing the wild-type lncRNA Gm33755 and circRNA6370 3′-UTR (Figure 8). Collectively, these results establish lncRNA33755 and circRNA6370 as targets of miR-21-5p.

LncRNA–miRNA–mRNA regulatory interaction network analysis.

Figure 6. LncRNA–miRNA–mRNA regulatory interaction network analysis.

CircRNA–miRNA–mRNA regulatory interaction network analysis.

Figure 7. CircRNA–miRNA–mRNA regulatory interaction network analysis.

Confirmation of the relationships. (A) Relative luciferase expression of wild-type and mutant lncRNAGM33755 UTR-bearing luciferase vectors cotransfected with miR-135b expression vectors. (B) Relative luciferase expression of wild-type and mutant circRNA6370 UTR-bearing luciferase vectors cotransfected with miR-135b expression vectors. n=6, ***P

Figure 8. Confirmation of the relationships. (A) Relative luciferase expression of wild-type and mutant lncRNAGM33755 UTR-bearing luciferase vectors cotransfected with miR-135b expression vectors. (B) Relative luciferase expression of wild-type and mutant circRNA6370 UTR-bearing luciferase vectors cotransfected with miR-135b expression vectors. n=6, ***P<0.001.

Discussion

With the aging of the world population, increasing numbers of elderly persons are sustaining vertebral compression fractures (VCF) due to osteopenia, 25% of postmenopausal women are affected by a compression fracture during their lifetime, and spinal cord injury is one of the most serious complications of VCF in elderly patients leading to significant morbidity and mortality [1618].

Since there are no approved therapies for restoring sensation or mobility following TSCI, achieving functional rehabilitation has been among the primary research interests of experimental neuroscientists in recent decades [19]. TSCI is a two-step process that can cause a permanent loss or reduction in bodily function below the level of the lesion site. The primary damage is the mechanical injury itself, and the secondary damage results from biochemical processes following the primary damage [3]. Physical trauma causes rupture of the blood-spinal cord barrier in the lesion epicenter, leading to hemorrhage, ischemia and inflammation, followed by local neuronal and glial cell death [5]. The nonneural damage in the lesion core ultimately resolves into a cavity surrounded by astrocytic and fibrotic scar borders [20]. Axonal regeneration is a complex procedure that includes structural synapse remodeling, axonal sprouting and regrowth across the lesion. A reduced intrinsic growth capacity, the absence of external growth stimulation and the presence of external inhibitory factors could lead to the failure of axons to regrow spontaneously across severe tissue lesions [21, 22]. The lesion compartments consist of different cell types, and cell biology influences axonal growth and regrowth in different ways [23]. Alleviating these differences is fundamental for achieving or improving axonal regeneration and designing rationally targeted interventions. Although several biomolecules are being used as diagnostic or prognostic biomarkers and therapeutic targets, they do not have sufficient accuracy or sensitivity to recognize pathogenesis, guide therapy or evaluate prognosis [5]. Since TSCI is a multifaceted pathological process, it is unlikely that any one molecule or pathway can affect the large number of obstacles that occur following trauma. Indeed, this may be the reason why many disparate treatments generate similar levels of recovery in TSCI animal models; as such, constructing the regulatory network involved with TSCI is clinically significant. In this study, we demonstrated that the expression of related ncRNAs and mRNAs significantly changes in the spinal cord tissue after traumatic injury, and we predicted the structure and potential functions of the regulatory network associated with these DE ncRNAs and mRNAs.

Recent studies have revealed the involvement of some specific miRNAs in many types of neuronal function in diseases, such as axon regrowth and neurodegeneration. MiRNAs are considered to be one of the major factors in the pathogenesis of CNS injury because of their intrinsic properties in regulating several biological functions and their potentially large impact in RNA disorders [24]. In this study, a marked dysregulating occurs in the expression of axon regrowth- or regeneration-associated mRNAs after injury, such as STAT3, p53, c-Jun, FOXO, KLFs and Sox, as well as their target DE miRNAs, miR-125b, miR-9, miR-222, miR-21, miR-135b and miR-145, respectively. In addition, miRNAs are critical regulators of the main molecular cascades regulating axonal growth, i.e., miR-26a and miR-222 repress the PTEN pathway, miR-124 targets the GSK-3b pathway, miR-9 targets the MAP1B–Rac1 pathway, and miR-133 inhibits the rhoA–PI3K–AKT pathway. MiRNAs also participate in inflammation, apoptosis and myelination-related lesions in neurological damage disease, i.e., let-7 inhibits IL-6 during inflammation, miR-29b increases proapoptotic gene expression, and miR-138 regulates myelination-related lesions.

Most strikingly, ceRNAs, which include lncRNAs, circRNAs and pseudogenic RNAs, cross-regulate each other by competing for shared miRNAs on miRNA response elements (MREs) [25]. CeRNA crosstalk is a type of posttranscriptional regulation that is mediated by miRNAs and links the functions of coding and noncoding RNAs. LncRNA can directly regulate the structure of DNA and the transcription and translation of RNA; notably, lncRNA can act as an miRNA sponge to competitively bind miRNA [14]. CircRNAs were later identified and are more enriched in neuronal tissues than other tissues because the long introns of neuronal genes promote circRNA formation [26]. The ceRNA regulation network plays a critical role in central neuropathy, i.e., the lncRNA SNHG5–KLF4–eNOS axis enhances the viability of astrocytes and microglia [27], the circRNA2837–miR-34a axis protects neurons against injury by inducing autophagy [28], and the circRNA ZNF609–miR-615–METRN axis reverses retinal neurodegeneration [29]. Recent research has even generated a complicated ceRNA network demonstrating that lncCyrano–miR-7 prevents the cytoplasmic destruction of circCdr1while repressing miR-671–circCdr1splicing in the brain [12]. In addition, to explore the roles of ncRNAs through this potential mechanism, we performed GO and KEGG pathway analysis to annotate predicted target mRNAs and predominant pathways of the differentially expressed ncRNAs. In this study, we found that the target mRNAs are involved in multiple biological processes, cellular signaling pathways, protein activities and gene splicing after SCI. Strikingly, focal adhesion was noted to be one of the most significantly enriched and meaningful terms of biological processes in both ncRNAs and mRNAs after GO analysis, and phagosome pathways was found to be one of the most significantly enriched and meaningful pathway of ncRNAs and mRNAs groups after KEGG pathway analysis.

In previous work, we demonstrated that miR-21 regulates astrogliosis through the PI3K–Akt–mTOR pathway and regulates fibrosis through the TGF-β–Smad pathway after TSCI [33, 30]. The reactive astrocytes and fibrotic scar tissue formed by perivascular cells stabilize the outer borders of the initial lesion epicenter and act as a chronic, physical, and chemical-entrapping barrier that prevents axonal regeneration [31, 32]. Our previous results suggest that miR-21 knockdown significantly suppressed scar formation and improved motor functional recovery after TSCI [30, 33]. In this article, a ceRNA regulation network with miR-21 as the center was constructed and provides initial evidence of the binding relationships between lncRNA GM33755, circRNA 6730 and miR-21. However, further specific studies are needed to determine whether lncRNA GM33755 and circRNA 6730 play roles in pathological processes after neurological injury.

Recently, our understanding of the extrinsic and intrinsic factors that block axonal regeneration or neuroplasticity has immensely improved with further illumination of the molecular events that occur following TSCI in mouse models. However, it is important to keep in mind that variations in axonal regeneration responses exist between mice and humans. In this article, we choose 3 days post-SCI as our time point, acute phase is the original point of pathological process. In acute phase, cell death and eventually axonal die back. Inflammation and axon regrowth at later time point, in sub-acute or chronic phases, then astrogliosis and fibrotic scar are done, axon struggle to regeneration. TSCI has an acute phase and a chronic phase, what should we do next is to compare different injuries severity in acute time points and chronic time points. Future efforts should be made to mimic endogenous ncRNA deregulation and determine the functions of ncRNA interactions in the context of posttranscriptional regulation as a whole. Findings arising from such studies will expand our understanding of the potential of ncRNAs, offer novel insight into the identification of new biomarkers and help guide strategies toward the development of potential therapies to enhance axonal regeneration and functional recovery during acute and chronic stages following TSCI. The spinal cord rarely repairs itself after an injury, but methods for promoting axonal regeneration are on the horizon.

Materials and Methods

Mouse SCI model and experimental groups

All animal protocols were approved by the Research Ethics Committee of Shandong University (Jinan, China). All tissue samples of the SCI epicenter and parameters of Allen’s weight-drop apparatus were obtained as described in our previous study [33]. Briefly, 64 adult male C57BL/6 mice were randomly divided into the SCI (n=8)×4 and sham (n=8)×4 groups. Mice in the SCI group underwent T8–10 vertebral laminectomy to expose the spinal cord after anesthesia was induced with 3% pentobarbital. Then, moderate SCI was induced using a modified Allen’s weight-drop apparatus, subsequently, the muscles and skin were sutured layer by layer. Mice in the sham group underwent the same laminectomy but without SCI. The movement ability of randomly selected mice (SCI group n=3, sham group n=3) was evaluated using the Basso motor score (BMS) for 3 days. Three independent and well-trained investigators scored the animals according to the standard guidelines. Then, the final score was recorded as the average of the investigators’ scores. First batch of 16 C57BL/6 mice divided into the SCI group (n=4) and sham group (n=4), were humanely sacrificed on days 1 and 3 postsurgery, respectively, for the collection of T8–10 spinal cord tissues. The spinal cord lesions were analyzed after tissue samples were fixed, washed, dehydrated, cleared, embedded, frozen and stained with H&E as previously described [33]. The remaining 48 mice, SCI (n=8)×3 and sham (n=8)×3, were humanely sacrificed on day 3 postsurgery for the collection of T8–10 spinal cord samples, then used to extract total RNA for Sequencing and qRT-PCR validation.

RNA library construction and sequencing, transcript abundance estimation and differential expression analysis

Total RNA of spinal specimens was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The procedure mainly included homogenization, phase separation, RNA precipitation, washing, solubilization, and monitoring of RNA degradation. The RNA concentration and quality were measured by UV absorbance at 260/280 nm; then, a LabChip Kit (Agilent, CA, USA) was used to analyze the RNA integrity.

Small RNA single-end sequencing was performed on an Illumina HiSeq 2500 LC-BIO system (Hangzhou, China). CircRNA and lncRNA paired-end sequencing were performed on an Illumina HiSeq 4000 LC-BIO system (Hangzhou, China). The fragments per kilobase of exon per million fragments mapped (FPKM) was used to measure the relative abundance of the transcripts after aligned read files were processed by in-house scripts. The expression levels of lncRNAs and mRNAs were determined by StringTie (http://ccb.jhu.edu/software/stringtie). CIRCexplorer was used to measure the expression of circRNAs [34], unique circRNAs were generated from assemblies.

qRT-PCR validation

As previously described [35], total RNA was reverse-transcribed into cDNA; then qRT-PCR was performed using an Applied Biosystems (Wilmington, DE, USA) 7500 RT-PCR system. GAPDH was used as an internal control to normalize relative circRNA, lncRNA and mRNA expression levels. MiRNA expression levels were normalized using U6. The 2–ΔΔCT method was used for comparative quantitation. Three independent experiments were performed. The specific primers for each gene are listed in Table 5.

Table 5. Primers designed for qRT-PCR validation.

GenePrimer
H19F TTCACTTAGAAGAAGGTTCA
R TTCCATTCTCCAGTTATTGA
Gm12840F CCAAGGAGTTGACTGATTATCT
R ACACAAGCAAGACCAATACA
Gm26809F ATCTCTAAGCACACTCGTCCAC
R ACTAATCGCCGCCGTCAG
circRNA7010F CTGGAGACTGTGGAAAGC
R TGTAAGGACACTGGGGC
circRNA2464F CTGTCAAGTATGTGGAGTG
R CAACAGCACCATCACC
circRNA7435F ATGACATCCGCAGAAGG
R AGGCAAATACCGCACTC
mmu-miR-21a-5p_R+1F CGGGCGTAGCTTATCAGACTG
RT GTCGTATCCAGTGCAGGGTCC
GAGGTATTCGCACTGGATACGACGTCAAC
mmu-miR-423-3pF TTAGCTCGGTCTGAGGCCC
RT GTCGTATCCAGTGCAGGGTCC
GAGGTATTCGCACTGGATAC
GACACTGAG
mmu-miR-92a-3pF CCGTATTGCACTTGTCCCG
RT GTCGTATCCAGTGCAGGGTCC
GAGGTATTCGCACTGGATAC
GACACAGGC
Tmsb4xF AGAACTACTGAGCAGGAAGG
R GGACATCTTTGACCATCTTGAA
Lyz2F ATGAAGACTCTCCTGACTCTG
R ATAGTAGCCAGCCATTCCAT
Ftl1F TGGAGAAGAACCTGAATCA
R AGGAAGTCACAGAGATGAG
GAPDHF GGTGAAGGTCGGTGTGAACG
R CTCGCTCCTGGAAGATGGTG
U6F CTCGCTTCGGCAGCACATATACT
R ACGCTTCACGAATTTGCGTGTC

GO annotations and KEGG pathway analysis

GO annotations and KEGG pathway analysis were performed to investigate the potential roles of all DE ncRNAs. GO analysis includes three domains, cellular components, biological processes, and molecular functions, and provides a controlled vocabulary to describe DE mRNAs (P<0.05) in GO categories (http://www.geneontology.org) [36]. In addition, the KEGG database (http://www.genome.ad.jp/kegg/) was used to detect the potential functions of the target genes in the identified pathways [37], with significance indicated by P-values <0.05.

Analysis of ncRNA regulatory network

An ncRNA regulatory network was constructed to examine the interactions and functional links among dysregulated mRNAs and ncRNAs in the pathological process of SCI. The target mRNAs of miRNA were predicted by software programs as previously described [35]. We selected the dysregulated target RNAs correlating to DE ncRNAs. Cytoscape software (San Diego, CA, USA) was used to construct interaction networks for lncRNA–miRNA–mRNA and circRNA–miRNA–mRNA.

Luciferase assay

293T cells were cultured in 94-well plates and cotransfected with luciferase reporter constructs containing lncRNA GM33755 or circRNA 6370 (LC) and a Renilla luciferase construct (Invitrogen), miRNA-21 mimic or scrambled negative control (LC) were transfected using Lipofectamine 2000 (Invitrogen) for 6 h. After 48 h of culture at 37°C, the culture supernatant was mixed with LAR II and measured using an illuminometer. Then, a luciferase activity assay was performed using a dual luciferase reporter system (E1910, Promega, Madison, WI, USA). In addition, Stop&Glo Reagent used as an internal control. The results shown represent the means of three experiments and are presented as the mean ± standard deviation (SD).

Statistical analysis

SPSS 20.0 (IBM, Chicago, IL, USA) and GraphPad Prism software (La Jolla, CA, USA) were used to perform the statistical analysis. Data are presented as the mean ± SD. ANOVA and Student’s t-test were used for comparisons (P<0.05). The Chi-squared 2X2 test, Chi-squared nXn test and Fisher’s exact test were used to assess the differential expression of miRNA (P<0.05), differential lncRNAs expression was examined using the R package Ballgown (P<0.05) and CircRNA expression in the different samples and groups was calculated using scripts developed in-house (P<0.05).

Ethical disclosure

The authors state that the study was approved by the Ethics Committee of Jinan Central Hospital Affiliated to Shandong University.

Author Contributions

W.W. performed the experiments. Y.S. and T.S. were responsible for construction of the animal model, collection of the results and analysis of the ncRNA regulatory network. H.L. and W.X. were responsible for the PCR and luciferase assays. J.C., X.P. and L.S. participated in the data analysis. B.N. conceived of the study and participated in its design and coordination. All authors have read and approved the final submitted manuscript.

Conflicts of Interest

Grant support was provided by the National Natural Science Fund of China (Nos. 81401014, 81771346), the Chinese Postdoctoral Science Foundation (No. 2014M561935), the Chinese Postdoctoral Science Foundation (No. 2015T80725), and Technology Research and Development Program of Jinan City (No. 201704133).

References