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.

However, there is no effective treatment for these two types of patients because of the difficulty of recovery of sensation and function. 1,[3][4][5] Acute physical trauma is usually brief, and repair begins at the moment of damage. Neuronal, glial and angiocellular forms of epicentre necrosis occur acutely, followed by neuronal apoptosis. After 24 hours, apoptosis begins in glial cells, followed by oligodendrocytes, resulting in a non-neuronal lesion core. [6][7][8] Next, haemorrhage, ischaemia, reperfusion injury and inflammation develop. Therefore, transcription and expression of ribonucleic acid (RNA) change dramatically in the acute period. 1 To date, scientists in the neural regeneration field have focused on dissecting out the cellular and molecular mechanisms of these processes in vivo.
Although 1%-2% of the human genome encodes proteins, most nucleotides are detectably transcribed detectably. 9 Raquel's investigation is the first to emphasize the regulatory function of non-coding transcripts in SCIs. 10 Long non-coding RNAs (lncRNAs), composed of >200 nucleotides, lack an open reading frame (ORF) and cannot be translated into protein; however, they have multiple functions, including cis or trans transcriptional regulation, nuclear domain organization and the regulation of proteins or RNA molecules. 11 An important and major function revealed in this study was its ability to serve as competing RNA for microRNA (miRNA) binding. lncRNA has a spatio-temporal expression pattern in embryonic stem cells and brain development; it also participates in physiological processes such as neurogenesis, cell differentiation and maturation, myelin sheath formation and synaptic plasticity. 12 To date, few potential regulatory competitive endogenous RNA (ceRNA) networks have been revealed after SCI. One regulatory construct of note, lncRNA6032/miR-330-3p/Col6a1, has a potential regulatory effect in chronic SCI. 10 The regulatory function of lncRNA has been widely acknowledged, but systemic genome-wide analysis of ATSCI-related lncRNA expression profiles and the precise latent regulatory networks remain unclear. Hence, this study aimed to use high-throughput sequencing combined with a dedicated bioinformatics platform to comprehensively identify lncRNAs, miRNAs and messenger RNAs (mRNAs) expressed after ATSCI and to profile precisely the corresponding latent regulatory ceRNA network.

| Mouse model, sample extraction and staining
Forty-eight clean-grade, healthy, male C57BL/6 mice from the same litter were used for experimentation at the age of 8 weeks, each weighing 18-22 g. They were purchased from the laboratory animal centre of Sichuan University and were randomly stratified into six groups of 8:3 groups constituted the ATSCI cohort and the other three groups constituted the Sham cohort. Eight tissue samples were harvested from each group, one of which was randomly selected for haematoxylin-eosin (HE) staining, and the remaining were combined for RNA extraction. A standard Allen's drop SCI injury was induced as described previously. 13 In brief, laminectomy was performed to expose the dorsal aspect of the spinal cord (T8-T10) in SCI and Sham mice groups. Allen's drop injury (weigh of 6 g and height of 60 mm) was induced in the SCI group. The animals were sacrificed by cervical dislocation euthanasia after inducing anaesthesia with 3% pentobarbital. The spinal cord tissues at the level of the contusion injury were harvested on postoperative days 1, 3. The HE and total RNA groups were assembled 3 days post-injury. In the HE group, as previously described previously, 13 after being ground on ice and treatment with Trizol (Takara), the spinal cord samples were combined with lysis buffer and chloroform (Sinopharm), the middle layer centrifuged, and isopropanol (Sinopharm) added. Next, the supernatant was decanted, centrifuged with ethanol, subjected to dry precipitation and combined with RNAase-free water (Takara). Another batch of 6 miR-21-Knockout (21KO) mice (Model organisms) and 6 C57BL/6 mice underwent spinal cord contusion, and 6 C57BL/6 mice underwent laminectomy without spinal cord contusion. Locomotor activity was evaluated in an open field for 10 days using the Basso Mouse Scale (BMS). 13 Mice were sacrificed and spinal cord tissues were harvested at the lesion epicentre of the contusion injury at 10 days post-surgery for immunohistochemical staining.

| Construction of sequencing library and sequencing experiment
Total RNA purity was detected by a spectrophotometer (Invitrogen). Thereafter, miRNA, lncRNA, and mRNA sequencing and database construction were performed. A small RNA library was built with TruSeq Small RNA Sample Prep Kits (Illumina).

Single-end sequencing (36 or 50 bp) was performed on an Illumina
Hiseq 2500 high-throughput sequencing system (LC-BIO). For the lncRNA and mRNA library preparation, paired-end sequencing was performed using an Illumina Hiseq 4000 system (LC-Bio) in the >200 bp length range.  17 Coding-Non-Coding-Index (CNCI) 18 and Pfam, 19 the remaining transcripts with class codes i, j, o, u and x were considered to be lncRNAs. StringTie was used to determine expression levels for mRNAs and lncRNAs by calculating the fragments per kilobase of exon model per million reads mapped (FPKM).
Differentially expressed miRNAs, mRNAs and lncRNAs were selected with the inclusions log 2 (fold change) value >1, or log 2 (fold change) value <−1 and statistically significant (P < .05), using the R package Ballgown.

| Functional annotation enrichment analysis and prediction of networks
In Gene ontology (GO), 20

| Cell culture and mechanical treatment
Primary spinal fibroblasts (PriCells, Wuhan, China) were cultured as previously described. 13

| Luciferase reporter assay and fluorescence in situ hybridization
293T cells were inoculated into 96-well plates at a confluence

| Western blot analysis, immunofluorescence and immunohistochemistry
Total protein was harvested from cells and tissues as previously described, 13  Fibroblasts were stained for immunofluorescence as previously

| Statistical analysis
Statistical analysis was performed using GraphPad Prism software (GraphPad Software). Data are expressed as means ± standard deviation (SD). Student's t test was used to assess statistical significance, and a P-value <.05 was considered statistically significant. NAs are displayed as a heatmap ( Figure 2A) and a volcano plot ( Figure 2B). Expression profiling of mRNAs ( Figure 2C) and miR-NAs ( Figure 2D) are shown in a volcano plot based on P-values and fold changes.

| Chromosome distributions: comparison between long non-coding RNAs and messenger RNAs
Almost all known lncRNAs were distributed over all the chromo-  Figure 3B), while the ORFs of most lncRNAs were <100 nucleotides long ( Figure 3C). Moreover, within length >500 bp sets, lncRNAs were shorter than mRNAs, while within length <500 bp sets, lncRNAs were longer than mRNAs ( Figure 3D). Most lncRNAs had <6 exons, mRNAs had more exons and exon numbers distributed over a wider range. Moreover, numerous mRNAs had >9 exons ( Figure 3E). FPKM data indicated that many lncRNAs had transcription levels similar to mRNAs ( Figure 3F).

| Enriched ontology terms and KEGG pathways of differentially expressed lncRNAs
Gene ontology term enrichment can determine the main biological functions of genes via GO-function enrichment analysis. In this study, such analysis showed that dysregulated lncRNA were associated with the following: high-density lipoprotein particle receptor binding, negative regulation of very low-density lipoproteins, protein oxidation, lipase inhibitor activity, cellular responses to cobalt ions and cellular detoxification of nitrogen compounds ( Figure 4A).

| Prediction and construction of an lncRNA-microRNA interaction network
miR-21a-5p was the most highly significant expressed miRNA post-ATSCI, with log 2 (fold change) = 2.49 (Table 2). This is in agreement with preliminary research findings, where miR-21a-5p was identified as the key regulatory node of many pathophysiological processes post-ATSCI. 5 Table 2). After a comprehensive analysis of TargetScan and MiRanda scores, lncENS-MUST00000195880 was selected as a potential key lncRNA that could interact with mmu-miR-21a-5p and affect the pathophysiological processed following ATSCI (Table 2). Moreover, we reverse screened several differentially expressed miRNAs that could bind with lncENSMUST00000195880: miR-21a-5p, miR-28a-3p,

| Validation of differentially expressed ncRNA and confirmation of interaction relationships
To verify the accuracy of high-throughput sequencing, we used qRT-PCR analysis on the significant differentially expressed RNAs miR-21a-5p and ENSMUST00000195880. High-throughput sequencing showed that miR-21a-5p and ENSMUST00000195880 were significantly upregulated after ATSCI (Table 2); qRT-PCR results were consistent with sequencing results (Figure 5A,B). Binding was predicted between ENSMUST00000195880 and miR-21a-5p ( Figure 5C,D). We constructed ENSMUST00000195880 (mutant; MT) and (wild-type; WT) luciferase plasmids ( Figure 5E) to verify the targeted binding relationship. We found a significant signal difference upon binding between the WT lncRNA and miR-21a-5p; however, the difference was abolished by the mutation, confirming that ENSMUST00000195880 could bind to miR-21a-5p ( Figure 5F). Smad7 expression was significantly decreased after SCI ( Figure 5G), and Smad7 sequence contains the binding targets of miR-21a-5p ( Figure 5H). Prediction of targeted interaction of lncENSMUST00000195880-miR-21a-5p-smad7 in the lesion areas after SCI ( Figure 5I).

| Ethological and histological analysis of miR-21a-5p knockout mice
Hindlimb motor function was assessed using the BMS from day an improvement on day 2 and gradually returned to normal from day 3 to day 4 post-surgery. Five days after surgery, there was no significant difference in recovery between SCI and 21KO + SCI groups, but the 21KO + SCI group started to show a better improvement in motor function 6 days after SCI, which continued to increase over the next 3 days, gradually augmenting the distance with the SCI group ( Figure 6A). Immunohistochemistry on day 10 post-surgery showed that fibronectin was downregulated by miR-21a-5p inhibition, but Smad7 was simultaneously upregulated ( Figure 6B-H).

| The validation of regulatory network in vitro
Identification of primary spinal cord fibroblasts ( Figure 7A-B), and predominant localization of lncENSMUST00000195880 in primary spinal cord fibroblasts ( Figure 7C). In vitro SCI models of TGF-β1 stimulation were constructed: the expression of miR-21a-5p decreased after lncENSMUST00000195880 overexpression ( Figure 7D). Simultaneously, the expression of fibrosis-related proteins fibronectin and collagen I decreased after lncRNA overexpression ( Figure 7E). Overexpression of lncENSMUST00000195880 promoted Smad7 expression and inhibited Smad's pathway activation ( Figure 7F). Working model of the target interaction axis after SCI ( Figure 8).

| D ISCUSS I ON
This study describes lncRNA, miRNA and mRNA expression profiles, obtained via high-throughput RNA-seq analysis, of ATSCI spinal lesion epicentre samples from 3 ATSCI and 3 Sham groups.
Despite strong scientific interest in lncRNAs, only a few studies have focused on ATSCI and lncRNAs networks. In this study, differentially expressed lncRNAs, miRNAs and mRNAs in ATSCI spinal lesion epicentres 3 days post-SCI were identified using RNA-seq.
Sequencing data showed that ATSCI can dysregulate lncRNA, miRNA and mRNA expression. Chromosome distributions, ORFs, transcript abundance, exon numbers and lengths were compared between lncRNAs and mRNAs to study their origins. Functional GO and KEGG analyses indicated that some differentially expressed lncRNAs might play crucial regulatory roles in several mechanisms: D-glutamine metabolism, D-glutamate metabolism and high-density lipoprotein particle receptor binding. A lncRNA/miRNA interaction network was predicted, constructed and preliminary edified.
Thus, these results suggest that regulatory effects of altered lncR-NAs and their networks might contribute to the pathophysiology of ATSCI.
ATSCIs are mainly caused by traffic accidents and falling from heights but they can also be caused by vascular lesions, tumours and iatrogenic injuries and often result in paralysis, involuntary movements, incontinence and depression. 5 Recently, SCI patients have shown an ageing trend with the average age of injury rising from 28.7 to 37.6 years and the proportion of elderly patients rising from 4.7% to 10%. 23 Trauma leads to the death of nerve cells and the destruction of nerve connections, thereby disrupting the excitability of upper and lower conduction nerves. 24 The most common form of ATSCI is an immediate traumatic injury due to direct force on the spinal cord, disrupting the blood-spinal cord barrier and leading to vasogenic spinal cord oedema, haemorrhagic F I G U R E 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 transformation and disruption of axons and cell membranes. 25,26 ATSCIs include acute, subacute and chronic phases; the pathophysiology post-SCI is biphasic and can be divided into primary and secondary phases. 27 After the acute phase (within 3 days post-SCI), secondary injury processes become dominant; the acute phase is likely to be the most amenable to neuroprotective interventions as it is typically the earliest point at which patients arrive at an appropriate centre to receive treatment. 28 The delayed onset of the secondary injury phase is related to expression changes in many genes involved in vascular dysfunction, oedema, ischaemia, excitotoxicity, electrolyte shifts, free radical production, inflammation and delayed apoptotic cell death. 1 There is increasing evidence that ncRNA plays an important role in injury progression. 25,26 Over 70% of human genes are transcribed, but <2% are translated into proteins; most of the remainder are transcribed into ncRNA. 9 NcRNAs can be divided into two categories according to function: housekeeping ncRNAs, including small nuclear, ribosomal and transport RNA; and regulatory ncRNAs, including miRNA, circular RNA and lncRNA. [29][30][31] Increasing emerging evidence indicates that an lncRNA/miRNA network has critical roles in biology and aetiology. 32 For instance, lncRNA GAS8-AS1 suppresses papillary thyroid carcinoma growth through the miR-135b/CCND2 axis. 33 Particularly in the nervous system, miRNA and lncRNA can act as regulatory or hormone-like factors to affect communication between target cells through autocrine or paracrine pathways, thus exerting considerable influence on neurophysiology and axon regeneration. 26,34 The lncRNA SNHG5 enhances astrocyte and microglia viability by upregulating KLF4 in SCIs. 35 A network containing XR_350851 that regulates autophagy after SCI has also been discovered. 36 In our previous studies, we constructed in vitro TGF-β stimulation models to explore the pathophysiological changes of spinal fibroblasts after SCI, finding the optimal duration and concentration of TGF-β1-induced fibroblast activation to be 10 ng/mL for 48 hours. 13 We further demonstrated that miR-21a-5p acts as a positive factor for SCI recovery in the acute phase and regulates astrocyte activation of the PI3K/Akt/mTOR signalling pathway. [37][38][39] Additionally, we discovered that miR-21a-5p could bind with Smad7, and verified that miR-21a-5p promotes spinal fibrosis post-SCI via the TGF-β/Smad signalling pathway. 22,40 Most strikingly, the astrocytes and fibroblasts in the spinal cord have been found to be key players in axon regeneration after injury. 4,41 In the present study, the lncRNAs and miRNAs differentially research. 13 Further, miR-21a-5p KO mice were used to detect the effect of miR-21a-5p on post-ASCI motor function recovery. The observation period of 10 days post-ATSCI was selected because about the subinitial stage of fibrotic scar formation takes places 10 days post-injury. 1,24 The results showed that the motor function recovery of miR-21a-5p KO group was better than that of the con- lncENSMUST00000195880 overexpression inhibited miR-21a-5p but promoted Smad7 expression, thereby inhibiting activation of the TGF-β/Smad signalling pathway and fibrosis. The limitation of this study is that it only preliminarily proves the possibility of the existence of this regulatory network, but the precise regulatory relationship needs to be properly examined. In future studies, we will explore the specific biological roles, mechanisms, and signalling pathways of the lncRNA/miRNA/mRNA interaction network predicted by this study, with special focus on exploring and verifying inflammation-related pathophysiological functions in the acute phase of ATSCI.

ACK N OWLED G EM ENTS
This study received funding from: the National Natural Science Fund

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interests. The sponsors had no role in the design, execution, interpretation or writing of the study, or in the decision to publish the results.

E TH I C A L A PPROVA L
All experiments were reviewed and approved by the Ethics Committee of Sichuan University.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.