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

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.

AGING
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 proteincoding 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 closedloop 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.

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  1B-1D). To further determine whether ncRNAs were   chr4 14 down miRNAs and their target DE mRNAs are shown in a Venn diagram ( Figure 3B-3D).

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).

AGING
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 cytokinecytokine 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 cytokinecytokine receptor interaction, focal adhesion, phagosome, chemokine signaling, Regulation of actin cytoskeleton, lysosome, Cell cycle and Toll-like receptor signaling pathways, among others ( Figure 5D).

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 miRNAbinding 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- AGING 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 6-7). 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 wildtype lncRNA Gm33755 and circRNA6370 3′-UTR ( Figure  8). Collectively, these results establish lncRNA33755 and circRNA6370 as targets of miR-21-5p.

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 [16][17][18].
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 AGING 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 bloodspinal 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 AGING 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 AGING 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 myelinationrelated 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  AGING 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.

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 AGING Table 5. Primers designed for qRT-PCR validation.

Gene
Primer 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 reversetranscribed 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.

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, Chisquared 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