Aberrant LncRNA Expression Profile in a Contusion Spinal Cord Injury Mouse Model

Long noncoding RNAs (LncRNAs) play a crucial role in cell growth, development, and various diseases related to the central nervous system. However, LncRNA differential expression profiles in spinal cord injury are yet to be reported. In this study, we profiled the expression pattern of LncRNAs using a microarray method in a contusion spinal cord injury (SCI) mouse model. Compared with a spinal cord without injury, few changes in LncRNA expression levels were noted 1 day after injury. The differential changes in LncRNA expression peaked 1 week after SCI and subsequently declined until 3 weeks after injury. Quantitative real-time polymerase chain reaction (qRT-PCR) was used to validate the reliability of the microarray, demonstrating that the results were reliable. Gene ontology (GO) analysis indicated that differentially expressed mRNAs were involved in transport, cell adhesion, ion transport, and metabolic processes, among others. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that the neuroactive ligand-receptor interaction, the PI3K-Akt signaling pathway, and focal adhesions were potentially implicated in SCI pathology. We constructed a dynamic LncRNA-mRNA network containing 264 LncRNAs and 949 mRNAs to elucidate the interactions between the LncRNAs and mRNAs. Overall, the results from this study indicate for the first time that LncRNAs are differentially expressed in a contusion SCI mouse model.


Introduction
Spinal cord injury (SCI) has become a global burden, influencing the quality of life and triggering serious socioeconomic consequences. Having SCI means being confined to a wheelchair and a lifetime of medical disease [1]. Repairing SCI is challenging due to multiple factors, including extensive cell loss, axonal disruption, growth-inhibiting molecules in the scar, and a lack of growth-promoting molecules [2]. Several studies have indicated that changes in various cellular events are significantly associated with SCI. Among these changes is the dysregulated gene expression of specific molecules. Large-scale gene expression studies had revealed that numerous protein-coding genes are differentially expressed in SCI, a portion of which was demonstrated to play pivotal roles in SCI [3,4]. However, abnormal gene expression is a highly complex process.
Recently, increasing evidence has indicated that noncoding RNAs possess significant regulatory functions in SCI [5][6][7][8][9]. Long noncoding RNAs (LncRNAs) have various functions [9] and range from 200 bp up to several kilobases in length; this length was determined from a convenient practical cutoff in RNA purification protocols which excludes other RNAs [5]. Several studies have revealed that LncRNAs are especially enriched in the central nervous system and are implicated in several neurological diseases [10][11][12].
To date, no studies have focused on the differential expression of LncRNAs in an SCI model. In this study, we surveyed the temporal expression of LncRNAs in the spinal cord following contusive SCI in ICR mice by microarray analysis and validated the results with quantitative realtime polymerase chain reaction (qRT-PCR). Our research provides the first evidence of an aberrant LncRNA expression profile in SCI.

Ethics Statement.
All experimental procedures were conducted in conformity with institutional guidelines for the care and use of laboratory animals, and protocols were approved by the Institutional Animal Care and Use Committee in Soochow University, Jiangsu, China. This institution approved this study (Permit number: SYXK2012-0045). Although the license was Chinese version, we upload the figure of license in the attached file (S1 license pdf). Mice were provided with food, sterile water, and the appropriate ambient temperature in all experiment procedures and executed using the broken neck method under the deeply anesthetized condition. After death, the mice were incinerated and buried in the assigned place. A small funeral was held in the end.

Animals.
Male ICR mice (20-25 g) aged 6-8 weeks were used in these experiments (Carvens Laboratory Animal Technology Co., Ltd., Changzhou, China). All of the mice were allowed preoperative environmental adaptation for 1 week with normal circadian rhythms (one 12-hour light-dark cycle) and had free access to water and food. The mice were bred in cages (3 mice/cage).

Spinal Cord Injury Surgery.
Before the surgery, the skin preparation by scissor and disinfection of the operative region by povidone iodine were performed. The mice were deeply anesthetized with intraperitoneal pentobarbital sodium (50 mg/kg body weight). According to the anatomy of the mice, T10 vertebrae were the highest of the back when the mouse was in prone position and then we labeled T8 and T12 vertebrae proximally and distally. Subsequently, an incision was made on the mid-line of the back over the spinous processes from T8 to T12 vertebrae. Then, we exposed and separated the paravertebral muscles from the vertebra. Following a laminectomy at the 10th thoracic spinal vertebrae (T10), the underlying spinal cord at T10 was exposed, and the spinal cord was contused with a Multicenter Animal Spinal Cord Injury Study (MASCIS) Impactor weight-drop device, which uses a 5 g weight impact rod dropped from a height of 25 mm to produce a reliable contused SCI model. Following the injury, the muscle and skin were closed with absorbable sutures. According to previous studies, the mice were divided into five groups (three mice/group, sham operation, 1 day after injury, 3 days after injury, 1 week after injury, and 3 weeks after injury) and provided with food, sterile water, and the appropriate ambient temperature after the surgery. Meanwhile, tramadol hydrochloride (50 mg/kg body weight) was intraperitoneal for the postoperative analgesia. The mice were given manual bladder evacuations twice per day. The spinal cord located 5 mm proximally and distally to the injury epicenter was removed for RNA extraction.

Total RNA Extraction and Quality Control Assay.
Total RNA was isolated using TRIzol reagent (Invitrogen, Canada) and purified using an RNeasy Mini Kit (Qiagen, Germany), which includes a DNase digestion treatment. RNA concentrations were determined based on the absorbance at 260 nm and quality control standards at 260/ 280 = 1.8-2.1 using a NanoDrop 2000 (Thermo, US). RNAs from mice in the same groups were mixed to obtain equal masses for the microarray preparation.

Microarray
Analysis. The Affymetrix GeneChip Mouse Exon 1.0 ST Array (Custom CDF), which contains 42,000 LncRNAs and 16,416 protein-coding transcripts, was employed in this study according to the manufacturer's protocols from the GeneChip platform by Affymetrix (Santa Clara, CA. US). The LncRNAs were collected from several well-known data sources, including RefSeq, Ensembl, UCSC, NOCODE, and the related data literature.
Affymetrix Expression Console (versions 1.3.1) was employed for data parsing, quality control, and LncRNA standardization according to the random variance model (RVM). The normalized signal intensities for the LncRNAs and mRNAs were filed as a log 2 ratio.

Quantitative Real-Time PCR Assay.
To further confirm the reliability of the microarray assay, qRT-PCR assays were performed. Briefly, single-stranded cDNA was synthesized using the RevertAid kit (Fermentas Life Science, Burlington, ON, Canada), according to the manufacturer's protocols, with random primers and 1 g RNA from the same samples used in the microarray. Real-time PCR was conducted using the SYBR Green q-PCR SuperMix (Bio-Rad, USA). The primers used are listed in Table 1. Each qRT-PCR reaction included 10 L SYBR Green q-PCR SuperMix, 0.5 L forward primer (10 M), 0.5 L reverse primer (10 M), and 1 L cDNA. The total volume was adjusted to 20 L with ddH 2 O. The following thermocycler parameters were used to generate the dissociation curve: (1) 95 ∘ C for 5 min; (2) 40 cycles of 95 ∘ C for 15 s, 60 ∘ C for 40 s, and 72 ∘ C for 20 s; and (3) 65 to 95 ∘ C. The 7500 System SDS software (ABI, USA) was used for acquiring the Ct values with manual thresholds. PCR amplifications were performed in triplicate for each sample. Gene expression levels were normalized relative to the expression of GAPDH using the ΔΔCT method. The gene expression levels between the groups were compared using 2 −ΔΔCT and Student's -test.
values < 0.05 were considered statistically significant.

Series Test of Cluster Analysis.
To profile the gene expression time series and ascertain the most probable set of clusters generating the observed time series, the series test of cluster (STC) algorithm was employed. The primary advantages of this method include the ability to take the dynamic nature of gene expression time series into consideration during clustering with a reliable method for identifying the number of distinct clusters [13,14]. Briefly, the null hypothesis was defined such that the value of any past or future time point was independent. According to the null hypothesis, the number of genes in each trend conformed to a binomial distribution. Next, the number of predicted genes originally belonging to various trends was calculated by the means of replacement, with the various binomial distribution parameters under the null hypothesis trend being satisfied. Finally, the trend level significance was determined by calculating the  probability that a variable obeyed the binomial distribution and was greater than or equal to the actual number of genes trends. The figures for the significant and nonsignificant trends were generated based on these results.

GO and KEGG Enrichment
Analysis. GO analysis was utilized to identify the main processes and functional categories involved in the differentially expressed enrichment as previously described [15][16][17][18]. This method analyzed the number of genes present in a category on the microarray. Specifically, two-sided Fisher's exact test and 2 test were used to classify the GO categories and multiple comparison testing was performed by computing the false discovery rate (FDR) to correct value ( < 0.0001). KEGG enrichment analysis was performed as described previously [19,20] to investigate whether differentially expressed genes share similar biological functions. Two-sided Fisher's exact test and 2 test were selected to identify the significant pathways, and the threshold of significance was defined using the FDR and value ( < 0.0001).
In the network analysis, the degree was defined as the sum of the links that one node has to all of the other nodes. The K-core in the graph theory was defined as a method of simplifying the analysis of graph topologies. The clustering coefficient was defined as a measure of the relationship between a gene and its neighboring genes. We calculated Pearson's correlation coefficient for each pair of genes and selected significant correlation pairs to build the network.

Microarray Data.
The genes with an FDR ≤ 0.05 and a fold change ≥ 2 were selected. The microarray data ( Table 2) indicated that the LncRNA and mRNA expression dynamically changed from the initial injury to 1 day, 3 days, 1 week, and 3 weeks after injury. We discovered that the number of differentially expressed LncRNAs and mRNAs at 1 day after injury was similar to day 0. However, the expression strongly changed at 3 days and 1 week and appeared to decline at 3 weeks after injury. We integrated the differential LncRNA and mRNA expression at five time points and created a union set (Tables S1 and S2

qRT-PCR Validation Array.
To further validate the reliability of the microarray data, qRT-PCR was performed. We used the NONCODE database (http://www.noncode.org/) to select eight LncRNAs that are highly expressed in the central nervous system (NONMMUT052085, NONMMUT-038925, NONMMUT067118, NONMMUT051225, NONM-MUT005924, NONMMUT070015, NONMMUT021928, and NONMMUT061607). As shown in Figure 2, the variation from the RT-PCR was similar to that of the microarray. Student's -test indicates an agreement between the microarray data and qRT-PCR results with the exception of LncRNA NONMMUT052085. This finding indicates that our microarray data were reliable and can be used for bioinformatics analysis in the subsequent steps.

Series Test of Cluster Analysis.
To further narrow the differential expression of the LncRNAs and mRNAs and seek the trend levels of significance among the differential expression profiles, STC was used to detect temporal expression patterns of significantly differentially expressed genes and to identify the cluster of diverse genes that have similar expression patterns after SCI (Tables S3 and S4). As illustrated in Figures 3 and 4, the differential expression profile consisted of eighty differentially expressed LncRNA clusters and eighty differentially expressed mRNA clusters. Each cluster represented a model trend profile and contained genes with similar expression trends. Among these clusters, eighteen mRNA clusters and thirteen LncRNA clusters exhibited significant expression trends. These significant trends exhibited a variety of model trends, including stabilization, gradual augment, and decline.

Functional
Analysis of the Differentially Expressed mRNAs. GO and KEGG analyses were used to further evaluate the function of the genes that have significant expression trends. As shown in Figure 5, GO enrichment analysis, which identifies biological processes that the enriched transcripts are involved in, indicated that transport, cell adhesion, ion transport, metabolic process, innate immune response, and other significant biological processes were involved in SCI (Table S5). Additionally, KEGG enrichment analysis, which was performed to verify the significant module functions, revealed that neuroactive ligand-receptor interactions, the PI3K-Akt signaling pathway, focal adhesions, metabolic pathways, osteoclast differentiation, lysosomes, and other significant signal pathways were related to SCI (Table S6).

Dynamic
LncRNA-mRNA Network. To discover the significant molecular mechanisms of the LncRNAs associated with the pathology of SCI, a dynamic LncRNA-mRNA network that contained 264 LncRNAs and 949 mRNAs was constructed ( Figure S1 and Table S7). Interaction between these genes is noted when the quantification of LncRNA and mRNA interaction coexpression is greater than or equal to 0.997. In these maps, we identified various key LncRNAs with high degrees and K-cores, such as NONM-MUT038518, ENSMUST00000145363, NONMMUT001318, NONMMUT035870, and NONMMUT054688, which might play important roles in SCI pathology. In this map, the size

Discussion
This study was the first to report the LncRNA expression profile in a mouse contusion SCI model. SCI was a complex biological process that involves many molecules and cell events. Although LncRNAs were once deemed the "noise" of the genome, LncRNAs have been shown to be enriched in the central nervous system and exerted a significant effect on its growth and development [24]. Hence, understanding the LncRNA expression profile was crucial for exploring potential LncRNA functions in SCI pathology.
As described in previous studies [25][26][27], the mechanical forces imparted on the spinal cord cause primary damage to the neural tissue, but a complex cascade of pathophysiologic processes imperiling adjacent, initially spared tissue rapidly caused secondary damage following the initial event. Such studies [28,29] have identified a number of interrelated The upper number in the coordinate axis first is the number of profile; genes assigned represents the number of actual differential expression genes in this trend; genes expected stands for the theoretical number of LncRNAs in this trend according to the random distribution; adj value represents the significant level of the ratio of the actual LncRNAs number to the theoretical number of genes in this trend. We defined the significance at adj value < 0.05.  3  4  0  1  2  3  4  0  1  2  3  4  0  1  2  3  4  0  1  2  3  4   0  1  2  3  4  0  1  2  3  4  0  1  2  3  4  0  1  2  3  4  0  1 Figure 4: STC analysis for differential expression mRNAs related to SCI. (a) as is similar to Figure 3; it represents the expression patterns of these differential expression mRNAs. The coordinate axis and upper number were the same as Figure 3 processes that are thought to contribute to the secondary damage after spinal cord primary injury, including alterations in microvascular perfusion, free radical generation and lipid peroxidation, necrotic and apoptotic cell death, and ionic homeostasis dysregulation. Among these changes, the most important alteration was necrotic and apoptotic cell death.
Several studies [29,30] have suggested that cell fate peaks at 1 day and 7 days after injury. Based on this observation, we added two time points at 3 days and 21 days after injury to identify the dynamical gene expression changes that occur.
In this current study, the differential peaked at 7 days after injury. This finding might offer a clue for future research.
Animal models continued to play critical roles in the development of experimental research for SCI [31]. Injury reproducibility was an important characteristic of experimental SCI models because it limited the variability in gene expression outcomes. In this study, the contusion spinal cord injury model was classic and could mimic a clinical situation by using a MASCIS Impactor weight-drop device that employs a compression that causes bony fragments or extruded disk materials [32]. We produced a steady and reliable animal model, pooled RNA from injured spinal cords from 3 animals in each group to obtain total cellular RNA, and mixed the total cellular RNA at equal masses for the GeneChip preparations. However, only one GeneChip was subsequently used for each group, and we expanded the sample size (3 samples in each group) to guarantee reliable and precise microarray results, which correlated with the qRT-PCR results. To increase the reliability of the microarray results, the GEO database was used to identify results similar to those of our study. GSE45006 supports the reliability of our microarray results by blasting the mRNAs expression levels.
We chose the significant profile 12 by series test of cluster analysis to blast to GSE45006 data and found they had 74.8% of the same change tendency of mRNAs. After determining that the microarray results were reliable, RNAs that had significant changes in expression level (fold change ≥ 2) were chosen for the bioinformatics analysis. Cluster analysis of gene expression dynamics could offer linearity change trends for increasing and decreasing gene expression for future studies. This method represents gene expression dynamics as autoregressive equations, uses an agglomerative procedure to search for the most probable set of clusters given the available data, and considers the dynamic nature of gene expression time series during clustering with a reliable method for identifying the number of distinct clusters. Through this STC approach, one can extract gene sets with increasing and decreasing expression for specific studies.
For significantly dysregulated mRNAs, we examined potential LncRNA functions using GO enrichment and pathway analysis. GO analysis revealed that many of the genes with changed expression profiles were involved in transport, cell adhesion, ion transport, metabolic process, and innate immune response. Similarly, KEGG enrichment analysis demonstrated that the changed genes were involved in neuroactive ligand-receptor interactions, the PI3K-Akt signaling pathway, focal adhesion, and metabolic pathways. Then, subsequently, the dynamic LncRNA-mRNA network was constructed. This network was helpful for understanding the possible interactions and relationships between the differentially expressed LncRNAs and mRNAs.
The PI3K-Akt signaling pathway was classical pathway which was involved in apoptosis, protein synthesis, metabolism, and cell cycle and evidences indicated that this activated pathway could help in nerve regeneration [33,34]. The mRNAs which associated with this pathway and have been contained in the dynamic LncRNA-mRNA network were extracted. For example, FGF1 was a neurotrophic factor which had high expression level in gliocyte cell and a powerful neuroprotective and neuroregenerative factor of the nervous system [34,35]. In our study, this mRNA expression level was declined and had a dynamic relation with two LncR-NAs (ENSMUST00000138093 and ENSMUST00000129688) in the network. Both had high expression level and declined in the SCI pathology. In humans, FGF1 was regulated by HOTAIR via upregulating miR-326 expression, forming the HOTAIR-miR-326-FGF1 axis [35]. This gave us inspiration for the future work.
However, there were many other signal pathways which could infect cell motility, cell proliferation, cell differentiation, regulation of gene expression, and cell survival. One of them was focal adhesion signal pathway. With the same method, the mRNAs were extracted, containing THBS1, SPP1, CAV1, MAPK10, and so on. THBS1 was involved in neuronal migration and adhesion, neurite outgrowth, and synaptogenesis, which draws our attention. THBS1 could decrease neuronal excitability via reducing AMPA receptors (AMPARs) and increasing glycine receptors (GlyRs) in synapses [36]. Some researches suggested that THBS1 helped the recovery of normal synaptic activity after injury responses [37] and even was a potential therapeutic target in the pathogenesis of Alzheimer's disease [38]. Many studies [39,40] indicated that LncRNAs played key role in the neuronal degeneration disease which was involved in the progressive loss of structure or function of neurons, including death of neurons. In our study, THBS1 gradually increased and had contrary variation trend with the LncRNA NONM-MUT027272 in the dynamic LncRNA-mRNA network. This could offer the idea in the future work.
However, some limitations in the current study should be noted. For example, the small numbers of samples in the microarray limit its reliability. We decreased the individual differences by mixing RNAs from different samples in the same group and increased the number of qRT-PCR validation genes. Additionally, given our research methods, we merely predicted the function of differentially expressed LncRNAs and were unable to determine exactly how these LncRNAs regulate target gene expression. LncRNA function can be validated by overexpression and RNA interference approaches. The molecular mechanisms can be investigated by RNA immunoprecipitation and chromatin immunoprecipitation, among others. In the future, we aim to further investigate and focus on LncRNA molecular mechanisms in SCI repair.

Conclusions
Overall, we demonstrated the differential expression profile of LncRNAs after SCI. To the best of our knowledge, our results will be helpful for understanding SCI pathology.