Systematic analysis of co-expression network and gene expression profile of differentially expressed lncRNA-mRNA one week after spinal cord injury

Abstract

Spinal cord injury (SCI) is a severe event that occurs in the central nervous system (CNS). Long noncoding RNAs (lncRNAs) have already been proven to play key roles in different diseases of the CNS. Nevertheless, the roles of messenger RNAs (mRNAs) and lncRNAs one week after SCI remain to be elucidated. We established a rat model of SCI and examined the differential levels of mRNAs and lncRNAs, and then, the differentially expressed mRNAs (DE mRNAs) and lncRNAs (DE lncRNAs) were investigated. Gene Ontology (GO) functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of DE lncRNAs and mRNAs were subsequently carried out to explore crucial biological processes and signaling pathways. Moreover, we constructed protein–protein interaction (PPI) networks, the STRING database, and identified the 10 high degree core nodes. We also established a co-expression network by calculating the correlation between lncRNA expression profile and mRNA expression profile. The identified DE lncRNAs and DE mRNAs may potentially be critical targets for the diagnosis and treatment of SCI in the future.

Abbreviations: SCI: spinal cord injury; CNS: central nervous system; YLDs: years of life lived with disability; lncRNAs: long noncoding RNAs; mRNAs: messenger RNAs; DE: differentially expressed; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; GSEA: gene set enrichment analysis; RIN: Renewable Identification Number; PCA: principal component analysis; PPI: protein–protein interaction; MF, molecular function; BP: biological process; CC: cellular component; MCODE: molecular complex detection

KEYWORDS:

• Spinal cord injury; microarray; co-expression of lncRNA and mRNA; neuronal development, inflammation

Introduction

Spinal cord injury (SCI) is one of the most severely devastating events that occur in the central nervous system (CNS) and may result in complete or incomplete loss of motor and sensory function (Kang et al. 2018; Yuan et al. 2018). SCI is associated with high incidence of disability, and mortality rates; it was reported that the age-standardized incidence rate of SCI was 13 (11–16) per 100,000, and the years of life lived with disability (YLDs) caused by SCI was 9.5 (6.7–12.4) million in 2016 (James et al. 2019). High-energy injuries, such as high falls and motor vehicle accidents, were the leading etiologies of new SCI cases (Jain et al. 2015). In 2016, SCI-related direct costs were reported to exceed 45 billion United States dollars (USDs) per year in the USA (Miller and Herbert 2016), which resulted in a great economic burden for patients, families, and society as a whole. However, there are no curative treatments available in the clinic (Ahuja et al. 2017; Fan et al. 2018).

The complicated mechanisms of SCI include primary and secondary injury. Mechanical traumatic injury causes primary injury, including the disruption of the blood-spinal cord barrier, ischemia, edema, compromised vascular supply, and hemorrhage, followed by secondary injury, including inflammatory cell infiltration, proinflammatory cytokine release, and necrotic cell death; these phenomena lead to the apoptosis and necrosis of neurons and glial cells, which may result in the loss of neural circuits and demyelination (Tran et al. 2018). Secondary injury can be divided into the acute, subacute, moderate, and chronic phases, and the subacute phase of SCI lasts from 2 days to 2 weeks (Ahuja et al. 2017). In this phase, edema, vasospasm, and vessel thrombosis cause further ischemia. The ongoing infiltration of inflammatory cells leads to the formation of cystic microcavities and more cell death owing to the damaged extracellular architecture of tissues and cells. Moreover, the proliferation of astrocytes is also important in the changes in the microenvironment of this stage post SCI (Silva et al. 2014). These pathological processes may play crucial roles in the secondary injury of SCI, and this stage may be a potential treatment window. Therefore, it is necessary to further understand the pathological mechanisms that occur after SCI to develop new therapeutic measures to better promote the neurological function of SCI patients.

Long noncoding RNAs (lncRNAs), noncoding RNAs with lengths of more than 200 nucleotides, participate in several regulatory processes related to gene expression, including gene transcription and translation modulation, chromatin modification and DNA methylation (Chu et al. 2011; Wang and Chang 2011; Rinn and Chang 2012). It has been reported that lncRNAs may play a significant role in several CNS degenerative diseases (Yu et al. 2015). However, there are few studies on the roles of these molecules after SCI.

Thus, for this study, microarray analyses were used to investigate the differential expression of mRNAs and lncRNAs one week post SCI in a rat contusion model. Subsequently, functional enrichment analysis of these differential mRNAs was performed. In addition, the interaction between differentially expressed lncRNA (DE lncRNA) and mRNA (DE mRNA) was analyzed by bioinformatics. This study may be helpful for exploring the pathophysiological mechanism underlying SCI and to develop and formulate new therapeutic measures.

Materials and methods

Animals

Eighteen female Sprague–Dawley (SD) rats (220 ±10 g) used in the experiments were obtained from the Laboratory Animal Center of the Academy of Military Medical Sciences (Beijing, China). The rat experiments were approved by the Ethics Committee of Tianjin Medical University. All of the animals were housed in a humidity- and temperature-controlled environment with a 12-hour light–dark cycle for 1 week before the operation and were allowed free access to food and water. The animal study was approved by the Experimental animal Welfare Ethics Committee, Tianjin Medical University General Hospital (Approval number: IRB2021-DW-74).

First, we used a random number table method to divide the experimental animals into two groups with equal numbers: the animals in the control group underwent laminectomy alone without SCI (n = 9), while the animals in the injury group underwent a T10 contusion of SCI (n = 9). The rats were weighed, deeply anesthetized with isoflurane (3% isoflurane for induction, 2% isoflurane for maintenance) at 25°C, and fixed in place. Laminectomy was carried out at the T10 vertebra level, and the contusive SCI model was established by an Impactor Model II (New York University, US) (a 10 g node fell freely from a height of 50 mm), as described previously (Ban et al. 2009). After hemostasis, the muscle and skin at the incision site was sutured. Successful SCl exhibits the following characteristics in the rats: formation of tail sway reflex, spinal cord ischemia, leg swings, and paralysis. If rats did not show the above symptoms after a single blow and the motor function of their hind limb remained normal on the second day after surgery, they were considered unsuccessful in the construction of a SCI model and would be excluded from the scope of this study. Cefuroxime was used for 3 days post injury to prevent infection. The bladders of the injured rats were manually emptied twice per day post surgery. One week after SCI, the animals were sacrificed after being deeply anesthetized with isoflurane, and a 1 cm section of damaged spinal cord was taken for microarray.

Extraction and purification of RNA

Segments of spinal cord tissues (1 cm long), including the injured epicenter, were collected for RNA extraction. Takara RNAiso Plus#9109 was used for the isolation of total RNA, and the Renewable Identification Number (RIN) was detected by the Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, US) to evaluate the integrity of the RNA. The purification of the qualified total RNAs was then performed using an RNase-Free DNase Set (Cat#79254, Qiagen, GmBH, Germany) and RNeasy mini kit (Cat#74106, Qiagen, GmBH, Germany).

Microarrays

Differential expression of mRNAs and lncRNAs was detected by SBC Rat (8*60 K) lncRNA microarray. The normalized mRNA and lncRNA expression profiles were listed in Supplementary Table 1 and 2. The amplification and labeling of the total RNA were performed by the Low Input Quick Amp WT Labeling Kit (Cat.# 5190-2943, Agilent Technologies, Santa Clara, CA, US). Purification of the labeled cRNAs was subsequently carried out using the RNeasy mini kit (Cat.# 74106, Qiagen, GmBH, Germany). Next, we used the Gene Expression Hybridization Kit (Cat.# 5188-5242, Agilent Technologies, Santa Clara, CA, US) to hybridize each slide with 600 ng of Cy3-labeled cRNAs. An Agilent Microarray Scanner (Cat# G2565CA, Agilent Technologies, Santa Clara, CA, US) was subsequently used to scan all the slides. R software 3.2.3 was then used for the identification of the differentially expressed lncRNAs (DE lncRNAs) and mRNAs. In our study, | log₂FC | > 1 and P value < 0.05 were considered to be mRNAs with significant differences. For lncRNA, we selected the criteria of | log₂FC | > 0.5 and P value < 0.05.

RT–PCR

RNA from spinal cord tissues was extracted using an RNA extraction kit (Takara RNAiso Plus#9109) and then converted into cDNA by a reverse transcription kit (Takara, Japan). Quantitative PCR was per­formed using the SYBR Green PCR Master Mix and QuantStudio 5 detection system (Thermofisher, USA). GAPDH, as an endogenous housekeeping gene, was used to normalize gene expression values. Each sample was measured three times for each gene. The primers used were shown in Supplementary Table 5.

GO, KEGG and GSEA analyses

Gene Ontology (GO) annotation terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were carried out to identify the enriched terms and pathways of the DE mRNAs. KEGG pathway analysis was used to investigate the significantly enriched pathways. R software with the cluster Profiler package was used to conduct GO and KEGG analyses. In addition, Fisher’s exact tests were carried out for the analysis of the distinct terms. A P value of < 0.05 was defined as statistically significant. In addition, we also performed GSEA analysis on all the mRNA expression matrices we detected, focusing on the biological process and KEGG pathway. We chose the | NSE | > 1, P value < 0.05 terms, with significant statistical significance.

PPI network and lncRNA and mRNA co-expression network analysis

Protein–protein interaction (PPI) networks of the DE mRNAs were constructed using a multiple protein web tool, the STRING database. Then, the PPI network results were analyzed using Cytoscape software (Smoot et al. 2011). Moreover, we used the Degree algorithm in Cytohubba plug-in of Cytoscape to calculate the 10 central genes. In addition, Mcode algorithm was used to get the first four clustering modules in the whole network. Pearson correlation coefficients for each pair of differential lncRNA and mRNA expression levels were calculated by R software. The relationship between Pearson correlation coefficient > 0.995 and P value < 0.05 suggested that lncRNA and mRNA are co expressed. The correlation heat map was constructed by R software. We visualized the co expression network of lncRNA and mRNA by Cluego, and analyzed the functional enrichment of mRNA.

Statistical analyses

SPSS 17.0 and R software were used to perform the statistical analysis. Statistical analyses of differences between groups were performed using a t-test, P value of < 0.05 was defined as statistically significant.

Results

Identified differential mRNA and lncRNA

Microarray technology was used to examine the differential levels of different mRNAs and lncRNAs one week after SCI. We first conducted principal component analysis (PCA) on all mRNA expression and lncRNA expression, and found that there was significant difference in overall data between the control group and the group one week after SCI, but no significant difference within the group (Figure 1A,B). For mRNA, we chose between the two groups of | log₂FC | > 1, in the case of lncRNA we chose between the two groups | log₂FC | > 0.5 the difference of mRNA and lncRNA P values < 0.05, the difference was statistically significant. We identified all 1312 mRNAs, including 296 differentially expressed mRNAs one week after SCI, of which 90 were up-regulated and 206 were down-regulated (Figure 1C). In addition, 110 lncRNAs were detected and their expression changes were studied one week after SCI, and 35 of them were found to have significant changes. Among them, 14 expressions were up-regulated and 21 were down-regulated (Figure 1D). We verified the top five lncRNAs with the most significant up-regulation and down-regulation trends by RT–PCR (Figure 1E,F). Finally, we showed the mRNA and lncRNA expressions in our samples by heat map (Figure 1G). In addition, the top 10 up-regulated and down-regulated lncRNAs are listed in Table 1.

Figure 1. Identified differential mRNA and lncRNA. A,B. Principal component analysis (PCA) of mRNA and lncRNA expression profiles. C,D. Volcano plot of DE mRNA and lncRNA. In the volcano plot, the X axis indicates the fold change, which is shown on a log₂ scale, and the Y axis indicates the P value, which is shown on a −log₁₀ scale. Red represents mRNA up-regulation and blue represents mRNA down-regulation. For mRNA, | log₂FC |>1 was selected, P < 0.05; for lncRNA, | log₂FC |>0.5 and P < 0.05 were selected, which were considered statistically different. E,F. RT-PCR to verify the top five lncRNAs with the most significant up-regulation and down-regulation trends. G. Heat map of DE mRNA and lncRNA. In the heat map, red indicates lncRNAs or mRNAs that were up-regulated, and blue indicates lncRNAs or mRNAs that were down-regulated.

Table 1. The top 10 up-regulated and down-regulated lncRNAs.

Up-regulated			Down-regulated
lncRNA	P value	log2(fold change)	lncRNA	P value	log2(fold change)
March1	1.35E-05	2.743295296	4930426D05Rik	0.000458608	−2.607306044
Dnm3os	0.002282278	2.180853923	Mirg	0.000433324	−2.412371428
Chmp4b	0.000222002	1.52105678	Ddc	1.22E-05	−2.208050112
Bend5	0.000137061	1.451610298	6720427I07Rik	0.001648962	−1.678405103
Wnk1	4.84E-05	1.403836495	Meg3	3.87E-05	−1.640526887
Gm26940	0.000302481	1.387629872	Shank1	2.63E-05	−1.513361431
Ski	0.000160067	0.95651036	Gm27028	0.002251338	−1.452638715
Tlk1	7.83E-05	0.694536169	Fam196a	0.00201839	−1.424264429
Gm26691	0.00095133	0.685618598	Mfsd2a	0.000432778	−1.220212534
Zfand6	0.000522904	0.617593732	Gm12843	9.21E-05	−1.205030073

Functional enrichment analysis of differential mRNA after SCI

GO annotation results showed that in DE mRNA, the function of up-regulated mRNA was mainly concentrated in the regulatory process of neuronal apoptosis and death, the response to external stimuli and the positive regulation of immune inflammation (Figure 2A). KEGG pathway analysis showed that these up-regulated mRNA may be related to PI3K-Akt signaling pathway and inflammatory signaling pathway (Figure 2B). The down-regulated mRNA was mainly enriched in the biological processes related to axon growth, indicating that the down-regulated genes of neurons one week after SCI may be related to the growth disorder after neuron injury (Figure 2C). However, in KEGG analysis, the down-regulated mRNA was related to calcium signal pathway, the pathway related to neuronal axon signal transmission and the inflammatory signal pathway of MAPK and Rap1 (Figure 2D).

Figure 2. Functional enrichment analysis of differential mRNA after SCI. A,B. GO enrichment analysis and KEGG pathway enrichment analysis were performed for up-regulated mRNA. C,D. GO enrichment analysis and KEGG pathway enrichment analysis were performed for down-regulated mRNA. E. GSEA analysis of all mRNA expression profiles after SCI revealed the biological process. F. GSEA analysis of all mRNA expression profiles for KEGG pathway after SCI.

Then we performed GSEA analysis on all the detected mRNAs, and focused on the biological process and KEGG pathway, so as to comprehensively understand the functional enrichment of all mRNAs one week after SCI. We found some possible main pathways, such as impaired growth and regeneration of neurons, activation of immune system, such as T cell activity and stimulation of corresponding cytokines at the same time, angiogenesis after SCI appeared (Figure 2E). In addition, through KEGG pathway analysis, we found inflammation related pathways such as TNF, activation of cytokine receptor related pathway, activation of IL-17 signaling pathway and activation of phagocytosis related pathway. The inhibition of calcium signaling pathway may lead to the decline of cell physiological function and cause cell injury (Figure 2F).

Construction of PPI network and GSEA analysis of hub gene and module gene function

According to PPI network results (Figure 3A), DE mRNA contained 228 nodes and 665 connections. Ten hub nodes, including Stat3, Snap25, Camk2a, Rbfox3, Dlg2, Fn1, Ngf, Vegfa, Isl1 and Erbb4, were calculated by degree algorithm of Cytohubba plug-in. In addition, we conducted cluster analysis on the whole network through Mcode plug-in, and obtained four core networks (Figure 3B).

Figure 3. Construction of PPI network and GSEA analysis of hub gene and module gene function. A. Protein-protein interaction networks of DE mRNAs. Nodes indicate mRNAs. Lines represent the interactions between DE mRNAs. Changes in DE mRNA expression level were indicated by color changes, with red indicating up-regulation and blue indicating down-regulation. The 15 hub genes obtained by Degree algorithm in Cytohubba plug-in are represented by yellow border. B. Cluster calculation through Mcode plug-in, and the first four core sub-networks with scores greater than 4 are obtained. Top left: Cluster1, Score 7.7; Top right: Cluster2, Score 4.25; Bottom left: Cluster3, Score 4; Bottom right: Cluster4, Score 4. C. Biological processes involving hub genes in GSEA analysis. D. The biological processes involved in GSEA analysis of mRNA in corresponding Cluster.

In GSEA analysis, we discussed the functions of hub gene and Mcode related core network related genes through statistical analysis. We found that the functional enrichment of hub gene was mainly related to the disorder of neural development and differentiation, and was involved in cell death and cytokine stimulation responses (Figure 3C). The biological process terms of GSEA enrichment analysis involving multiple hub genes are listed in Supplementary Table 3. Module 1 was associated with neural development and metabolic disorders, module 2 was involved in ion transport and cell signal transduction disorders, module 3 was involved in receptor-mediated endocytosis and regulation of cell death, and module 4 was involved in cell cycle and division disorders (Figure 3D). Table 2 lists the hub genes in detail.

Table 2. The hub genes and their corresponding degrees.

Gene	Degree	Gene	Degree	Gene	Degree	Gene	Degree
Stat3	27	Fn1	21	Shank1	16	Snap91	15
Snap25	26	Ngf	21	Camk2g	16	Ntrk3	14
Camk2a	25	Vegfa	20	Stxbp1	16	Nrxn1	14
Dlg2	23	Isl1	19	Bsn	16	Stx1b	13
Rbfox3	23	Erbb4	17	Cacng2	15	Syn2	13

Construction and functional analysis of lncRNA and mRNA co-expression network

In our data, by calculating the correlation between lncRNA and mRNA expression, we found the correlation between 14 lncRNA and 70 mRNA expression, including 87 relationship pairs (Pearson correlation coefficient > 0.995, P value < 0.05), which were shown by heat map (Figure 4A). The correlation between differentially expressed lncRNA and mRNA is listed in Supplementary Table 4. We take the function of mRNA co-expressed by these lncRNAs as the function of their regulation, and show the function enrichment through Cluego (Figure 4B). In addition, in order to understand the functions of these lncRNA mRNAs, similarly, we built the network through Rshiny (https://rshinyanalysisfigure.shinyapps.io/lncrna-mrna-net/). In the clustering of all enrichment pathways, the occurrence of cellular components accounted for the most, followed by learning and memory related functions, intercellular synapse related signal transmission, and angiogenesis (Figure 4C) Through the bubble chart, we listed the specific functions of lncRNA co-expressed mRNA, including the regulation of cell development, synaptic signal transmission and the regulation of cell differentiation (Figure 4D).

Figure 4. Construction and functional analysis of lncRNA and mRNA co-expression network. A. Pearson correlation coefficient and P value of lncRNA and mRNA were calculated to obtain the correlation between co-expressed lncRNA and mRNA. Displayed by correlation heat map. B. The lncRNA-mRNA co-expression network was established by Cluego, and the biological process and pathway of mRNA were analyzed. C. The proportion of genes enriched by each biological process pathway is displayed by pie chart. D. The first 15 GO enrichment items obtained from the biological process pathway analysis of mRNA in the co-expression network were displayed.

Discussion

SCI is one of the most serious traumatic diseases of the CNS that often leads to permanent disability and affects the functions of different systems, such as the nervous, respiratory, immune, and circulatory systems (Ahuja et al. 2017; Kang et al. 2018). Recently, with the continuous development of microarray technology, an increasing number of studies have shown that many lncRNAs participate in the process of SCI (Shi et al. 2018). Knockdown of lncRNA Neat1 can suppress the levels of inflammatory cytokines and the apoptosis of PC-12 cells in lipopolysaccharide-induced cell injury and promote the proliferation of PC-12 cells (Ban and Cui 2020). Furthermore, knockdown of lncRNA RNA Xist can reduce the apoptosis and inflammation of microglia after SCI through the miR-27a/Smurf1 axis (Zhao et al. 2020). However, to the best of our knowledge, the mechanisms by which lncRNAs affect SCI are still largely unknown. For this study, we identified the expression profile of mRNAs and lncRNAs one week post SCI by microarray technology. We found that 296 mRNAs and 35 lncRNAs were differentially expressed. Then, GO annotation enrichment analyses and KEGG pathway analyses were carried out to explore the different biological functions of the DE mRNAs, and the features of the DE lncRNAs were summarized. Moreover, core genes, which might be critical proteins one week after SCI, were investigated with the STRING database.

The GO analysis of differentially expressed mRNA and the comprehensive analysis of GSEA of all genes show that the obstacle of neuronal growth and differentiation, the regulation of apoptosis signal pathway, T cell activation, immune system activation and vascular development activation are the main terms. It is known that CNS axons constitute complex neural pathways within the spinal cord, where injured axons form bulbs at the terminals after SCI (Ruschel et al. 2015), but also where regrowing axons fail to cross due to local inflammatory responses as well as glial scar formation (Barnabé-Heider et al. 2010). Promoting endogenous axonal growth of neurons after SCI may be an effective means of restoring function after SCI by improving the microenvironment after SCI (Fan et al. 2018). Apoptosis is a caspase dependent programed cell death manner, and apoptosis related pathways can be detected to be activated after SCI (Springer et al. 1999). And in previous studies, apoptosis plays an important role in throwing in secondary injuries of SCI (Lu et al. 2000). While Gu et al found after knocking down lncRNA Xist could regulate the mTOR pathway to inhibit apoptosis after SCI (Gu et al. 2017), playing a neuroprotective role, suggesting that modulating lncRNAs can have an effect on apoptosis after SCI. The secretion of inflammation-related cytokines by activated microglia and the activation of inflammatory cells after SCI are major causes of secondary injury (Greenhalgh and David 2014). Tuning the chemokine CCL28 and its receptor ccr10 may have a potential role in reducing T cell recruitment after SCI to improve SCI recovery (Wang et al. 2019). Post SCI lesion sites trigger limited repair processes such as regeneration associated gene expression, endogenous vascularization, as well as axonal sprouting (Hagg and Oudega 2006). Whereas local angiogenesis due to local hypoxia as well as inflammatory stimuli is not sufficient to generate an effective supportive environment that promotes axonal regeneration (Loy et al. 2002; Ng et al. 2011). Interventions neuronal sprouting, inflammatory processes, as well as vascular remodeling after SCI are critical for functional recovery after injury.

In the combined GSEA and KEGG pathway enrichment results, we identified the inflammation-related pathways TNF, IL-17 signaling pathway cytokine receptor interaction activation, endocytosis activation, and calcium ion signaling pathway inhibition. As the spinal cord is subjected to mechanical violence that causes the structural destruction of local tissue cells, it triggers the release of inflammatory cytokines such as IFN-γ, TNF, interleukin (IL), etc. (Pineau and Lacroix 2007). These inflammatory factors can maintain an elevated state for about a week, stimulate the activation of microglia and induce macrophage infiltration in the injured area (David and Kroner 2011). In previous reports, the increased expression of TNF after SCI could not only induce the adhesion factor ICAM-1 and VCAM-1 to collect neutrophils into the lesion area, but also alter the permeability of vascular endothelial cells (Garcia et al. 2016). IL-17 can play a signal synergy with a variety of inflammatory cytokines after SCI (Ruddy et al. 2004). It can regulate the expression of VEGF and induce the expression of glial fibrillary acidic protein (GFAP) in astrocytes by activating JAK/STAT pathway in vivo and in vitro (You et al. 2017). Myelin disintegration and myelin fragments will be produced immediately after SCI. These fragments can act as inflammatory response stimuli (Sun et al. 2010) and inhibit axon regeneration and oligodendrocyte differentiation, resulting in secondary injury (Quarles 2009). The microglia that remove necrotic cells and fragments in the injured area can reach the peak about one week after injury, and these fragments can be degraded after they reach the lysosome through endocytosis (Wang et al. 2020). In our data, calcium pathway related differentially expressed mRNA was significantly down-regulated. Calcium in neurons regulates biological processes such as synaptic signal transduction, cell growth and differentiation (Nicholls 2004). The imbalance of calcium concentration in mitochondria or endoplasmic reticulum caused by cell injury will trigger downstream neurotoxicity and apoptosis (Bano et al. 2005; Calì et al. 2011).

After analyzing the PPI network, we identified Stat3, Snap25, Camk2a, Dlg2, Rbfox3, Fn1, Ngf, Vegfa, Erbb4 and Isl1 as the top 10 core genes. Stat3 plays an important role in nervous system development, neural regeneration, its expression is elevated after SCI, and it is involved in the activation of microglia, a neurotoxic cell after injury, as well as in the up-regulation of reactive astrocytes to promote inflammatory responses after SCI (Zhang et al. 2021). After SCI, it leads to sensory loss below the level of the injury and Wang et al showed that Snap25 is mainly distributed in the nervous system and that it is down-regulated in the injured area after spinal cord transection. This study also supported that Snap25 is highly related to neural signaling and that Snap25 can also promote neurite outgrowth and elongation (Osuka et al. 2007). Camk2a belongs to the serine/threonine protein kinase family and plays a major regulatory role in calcium ion transport and neurotransmitter transmission. After SCI, Camk2a expression is markedly down-regulated, with emergence of autonomic dysregulation dominated by parasympathetic activity (Shi et al. 2020). Dlg2 expression is down-regulated after SCI and influences neurotransmitter transport during the subacute phase of reinnervation injury (Lin et al. 2018). Rbfox3 acts as a neuronal splicing regulator and can act during neuronal synaptogenesis (Kim et al. 2013). Previous studies have shown that as a nerve growth factor, Ngf can promote the growth, development, differentiation, maturation of central and peripheral neurons, maintain the normal function of the nervous system, and accelerate repair after nervous system injury (Keefe et al. 2017; Hu et al. 2020). Vascular endothelial growth factor A (Vegfa) has recently been shown to be neuroprotective. Vegfa stimulates angiogenesis by signaling to matrix metalloproteinases (MMPs) to disrupt the blood-spinal cord barrier (BSCB) and matrix (Figley et al. 2014). And accumulating evidence suggests that Vegfa can also have neurotrophic, neuroprotective, and neuroproliferative effects (Liu et al. 2010). In the CNS, Ngr1/Erbb4 signaling is involved in a wide range of processes, including neuronal migration, axonal guidance, and synaptic plasticity and is essential for normal myelination and sensory function (Tao et al. 2013). SCI can affect Ngr1/Erbb4 signaling, which may contribute to the aggravation of pathophysiological changes in major brain regions during SCI (Xue et al. 2019).

By constructing lncRNA-mRNA co-expression network, we studied the function of mRNA to understand the function of differential lncRNAs. We found that these differentially expressed co-expressed lncRNAs are mainly involved in physiological processes such as cell development and differentiation, angiogenesis and synaptic signal transduction. lncRNA Mirg is derived from chromosome 12 and undergoes up-regulation in the embryo and liver, as well as in the limbs (Sekita et al. 2006). William et al reported that a decrease in the expression level of Mirg appeared in mouse embryonic stem cells (ESCs) with impaired developmental potential (Schacker et al. 2020). While Cao et al., in overexpressing lncRNA Mirg negatively regulated miRNA-1897 by sponge mechanism, promoting bone cell development (Ling et al. 2019). The lncRNA dynamin 3 alternative strand (Dnm3os) was independent transcription unit in the intronic region of the dynamin 3 gene (Mitra et al. 2017). Studies have found that lncRNA Dnm3os is up-regulated in neurodegenerative disease, Huntington's disease (HD) models, regulates the expression of target genes by sponging miRNA-196-5p, and plays a pathogenic role in HD models, whereas it can reduce cell apoptosis and exert protective effects on cells by down-regulating Dnm3os (Dong and Cong 2021). We established a regulatory relationship network based on the correlation relationship between lncRNAs and mRNAs. By modulating the expression levels of these lncRNAs in the co-expression network, they may then affect the expression of multiple protein coding mRNAs, potentially having an important impact on improving cell growth and differentiation effects after SCI.

Despite bioinformatics analysis of microarray data, we clarified lncRNA functions in the subacute phase of SCI. But limited by the complicated regulatory network between lncRNAs and mRNAs, we still cannot identify the core lncRNAs after SCI. And the regulation of lncRNAs by experiments to reverse the suppressed spinal cord function at subacute phase may still present difficulties for us. In addition, the sample size of the present study was relatively small. Therefore, future research focus will be to expand the sample size and screen the core DE mRNAs and lncRNAs, and further explore the pathological mechanism of SCI.

Conclusion

In this study, mRNA and lncRNAs that are differentially expressed one week after SCI were identified by microarray analysis, and bioinformatics analysis was carried out. Moreover, the key enriched GO annotation terms and KEGG pathways, which may be crucial components of the processes that occur one week post SCI, were investigated. Our study indicates that Stat3, Snap25, Camk2a, Rbfox3, Dlg2, Fn1, Ngf, Vegfa, Isl1 and Erbb4 may contribute to the pathogenesis one week post SCI. Moreover, the results in the current study need to be further verified in future work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data supporting the findings of this study are available in the figshare repository at https://doi.org/10.6084/m9.figshare.20065661.v3.

Authors’ contributions

SQF, HXZ and XYL conceived and designed the study. YK, YJS, ZJS and QL, performed the experiments. YJS also revised this manuscript. These four authors contributed equally. YGL and MFH performed the experiments. JX and CD analyzed the data and drafted the Abstract and Introduction; YFL and CZ interpreted the data and wrote the remaining parts of this manuscript; HJ performed the experiments. All the authors read and approved the manuscript.

Additional information

Funding

This study was funded by Science and Technology Development Fund of Tianjin Education Commission for Higher Education (2018KJ078) and Natural Science Foundation of Tianjin City (20JCYBJC00680).