- Email: [email protected]
Integrated analysis of competing endogenous RNA (ceRNA) networks in subacute stage of spinal cord injury
Journal Pre-proofs Research paper Integrated analysis of competing endogenous RNA (ceRNA) networks in subacute stage of spinal cord injury Nanxiang Wang, Lei He, Yang Yang, Simin Li, Yuyong Chen, Zhenming Tian, Ye Ji, Yufu Wang, Mao Pang, Yang Wang, Bin Liu, Limin Rong PII: DOI: Reference:
S0378-1119(19)30830-3 https://doi.org/10.1016/j.gene.2019.144171 GENE 144171
To appear in:
Gene Gene
Received Date: Revised Date: Accepted Date:
8 June 2019 13 October 2019 15 October 2019
Please cite this article as: N. Wang, L. He, Y. Yang, S. Li, Y. Chen, Z. Tian, Y. Ji, Y. Wang, M. Pang, Y. Wang, B. Liu, L. Rong, Integrated analysis of competing endogenous RNA (ceRNA) networks in subacute stage of spinal cord injury, Gene Gene (2019), doi: https://doi.org/10.1016/j.gene.2019.144171
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 Elsevier B.V. All rights reserved.
Integrated analysis of competing endogenous RNA (ceRNA) networks in subacute stage of spinal cord injury Nanxiang Wanga, Lei Hea, Yang Yanga, Simin Lib, Yuyong Chena, Zhenming Tiana, Ye Jic, Yufu Wangc, Mao Panga, Yang Wanga, Bin Liua,1,*, Limin Ronga,1,* aDepartment
of Spine Surgery, The Third Affiliated Hospital of Sun Yat-Sen University, No.
600 Tianhe Road, Tianhe District, Guangzhou, Guangdong Province, People’s Republic of China bDepartment
of Cariology, Endodontology and Periodontology, University Leipzig, Liebigstr.
12, 04103, Leipzig, Germany cDepartment
of Orthopaedics, The Second Affiliated Hospital of Harbin Medical University,
No. 246 Xuefu Road, Nangang District, Harbin, Heilongjiang Province, People’s Republic of China 1These
authors contributed equally as the corresponding author.
*Corresponding
author.
(1) Limin Rong Email: [email protected] (2) Bin Liu E-mail: [email protected]
Highlights Integrated analysis of competing endogenous RNA networks involved in the subacute stage of spinal cord injury were lacking in the literature. Integrated bioinformatic analysis identified specific genes (Flna, ID3, HK2), transcription factors (Sp1, Trp53, Jun, Rela), lnRNAs (Neat1, Xist, Malat1), miRNAs (miR-16-5p, miR-1958, miR-185-5p) as involved in the pathogenesis of spinal cord injury at subacute stage. These genetic and epigenetic mechanisms could have potential for translation in diagnostic and therapeutic research concerning spinal cord injury.
Abstract This study aims to investigate the genetic and epigenetic mechanisms involved in the pathogenesis of subacute stage of spinal cord injury (SCI). Gene-expression datasets associated with SCI were downloaded from the Gene Expression Omnibus (GEO) database, and differential expression analyses were performed in order to identify differentially expressed genes (DEGs). Multiple network types were constructed and analyzed , including proteinprotein-interaction (PPI) network, miRNA-target network, lncRNA-associated competing endogenous RNA (ceRNA) network, and miRNA-transcription factor (TF)-target network. Cluster analyses were performed to identify significant modules. To verify the prediction accuracy of the in-silico identified molecules, qRT-PCR experiments were conducted. The results depicted the Ywhae gene as the hub gene with the highest degree in the PPI network. The ceRNA network identified specific genes (Flna, ID3, and HK2), miRNAs (miR-16-5p, miR-1958, and miR-185-5p), and lncRNAs (Neat1, Xist, and Malat1) as playing critical regulating roles in the pathological mechanisms of SCI. The miRNA-TF-gene interaction network identified four important TFs (Sp1, Trp53, Jun, and Rela). The miRNA-gene-TF interaction loops from the significant modules indicated that miR-325-3p can interact with the Asah1 gene and TF-Sp1 by forming a closed loop. The qRT-PCR experiments verified four selected genes (Flna, ID3, HK2, and Ywhae) and two selected TFs (Jun, and Sp1) as significantly up-regulated following SCI. The results indicated that four genes (Flna, ID3, HK2, and Ywhae), four transcription factors (Sp1, Trp53, Jun, and RelA), two miRNAs (miR-16-5p and miR-325-3p), and three lncRNAs (Neat1, Xist, and Malat1) are likely to be involved in the
molecular mechanisms underlying the subacute stage of SCI. These findings uncover putative pathogenic mechanisms involved in SCI and might bear translation significance for future research towards therapeutic development.
Abbreviations: ceRNA, competing endogenous RNA; SCI, spinal cord injury; GEO, Gene Expression Omnibus; DEGs, differentially expressed genes; PPI, protein-protein-interaction; TF, transcription factor; mRNAs, messenger RNAs; ncRNAs, non-coding RNAs; miRNAs, microRNAs; lncRNAs, long non-coding RNAs; GO, Gene Oncology; KEGG, Kyoto Encyclopedia of Genes and Genomes; WGCNA, weighted gene co-expression network analysis; ADEV, astrocyte derived extracellular vesicles; ESCs, embryonic stem cells; NSCs, neural stem cells; BMSCs, bone marrow stromal cells.
1. Introduction Spinal cord injury (SCI) is a devastating neurological disorder caused by accidents, violence, and falls, and can often lead to neurological deficits and disabilities [1]. The pathogenesis of SCI consists of three phases, including an acute phase, a subacute phase, and a chronic phase [2]. Cell-based therapy has been found to be most effective in improving functional outcomes when applied during the subacute stage as compared to other stages [3, 4]. Therefore, the subacute phase may comprise the best period for cell-based therapeutic interventions and result in improved rehabilitation [5]. As such, a molecular understanding of the mechanistic events that occur during this stage is of utmost importance. The pathophysiological events during SCI are understood to be tightly regulated by several specific genetic and epigenetic entities [6, 7], which include both protein-coding messenger RNAs (mRNAs) and non-coding RNAs
(ncRNAs). To unravel the biological processes underlying subacute SCI, the expression patterns of these molecules require clear elucidation.
Two main types of ncRNAs, MicroRNAs (miRNAs) and long ncRNAs (lncRNAs) can suppress each other and thus act as competing endogenous RNAs (ceRNAs). They also form a regulatory lncRNA-associated ceRNA network (lncRNAs-miRNAs-mRNAs) to regulate the target mRNAs of miRNAs [8]. Some microarray or sequencing studies [9-11] have identified certain miRNAs (miR-21, miR-486, miR-20) to be dysregulated and involved in SCI by regulating inflammation, cellular death, regeneration, or gliosis [12]. A recent review [13] that summarized the results of existing lncRNAs profiling studies in SCI [14-16], noted that several lncRNAs’ dysregulation plays a pivotal role in regulating inflammatory response (lncRNA LINC00305, LncRNA LINC00341) and angiogenesis (lncRNA HIF1A‐AS2, lncRNA MEG3). Although multiple gene expression studies have investigated ncRNA signatures following SCI, their results are variable, plausibly owing to heterogeneity in experimental animal models, severity, degree and location of lesions, experiment type and the platforms used. At the same time, an integrated view of the data generated in these studies would be valuable, which can be facilitated by systematic integration using bioinformatics approaches. A recent bioinformatic study [17] has performed ceRNA network analysis in SCI; however, this study did not address different phases of SCI. Considering its specific clinical significance, the present study investigated the genetic and epigenetic mechanisms implicated during the subacute stage of SCI, by constructing an integrated ceRNA network to identify critical genes, miRNAs, lncRNAs, transcription factors, and signaling pathways.
2. Materials and methods The flowchart illustrating the main steps of this analysis is shown in Fig S1. The flowchart consists of 5 steps – identification of differentially expressed genes and co-expressed genes, identification of pathways using functional enrichment analysis, screening for significant gene modules, construction of various networks (e.g., protein-protein interaction (PPI) network, miRNA-target network, lncRNA-miRNA-gene network, and miRNA-transcription factor (TF)target network), and experimental validation of significant molecules. The detailed method for
each step will be described as follows.
2.1 Acquisition of datasets Gene expression datasets from SCI in house mice (Mus musculus) were downloaded from the GEO (Gene Expression Omnibus) database (https://www.ncbi.nlm.nih.gov/geo/). Mice with induced spinal cord contusion were defined as the SCI group, while mice without SCI were defined as the sham group. For mice in both groups, spinal cord contusion was conducted at thoracic level 8 (T8) with a moderate impact applied with the IH Spinal Cord Impactor (Precision Systems and Instrumentation, USA).
2.2 Differential expression and functional enrichment analysis Differential expression analysis was performed to identify differentially expressed genes (DEGs) from each of the selected datasets by using the R package ‘limma’ and DEGs that were common to all 3 datasets were identified. Next, functional enrichment analysis was of the shared DEGs identified enriched GO (Gene Ontology) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways, performed using the online resource DAVID v6.7 (https://david.ncifcrf.gov).
2.3 Construction of miRNA-target network and ceRNA network The miRNA-target datasets for Mus musculus were downloaded from starBase v3.0 (http://starbase.sysu.edu.cn/). The interaction pairs between miRNAs and their target DEGs were obtained in order to construct a miRNA-target interaction network. Furthermore, target lncRNAs of these miRNAs were extracted from the starBase v3.0 database. The lncRNAsmiRNAs-DEGs interaction pairs were used to construct a lncRNA-associated ceRNA network. Cluster analysis of the ceRNA network was performed by using the ‘MCODE’ plugin in Cytoscape, following which, modules representing densely connected nodes were extracted.
2.4 Construction of the miRNA-TF-target network Transcription factor (TF)-gene interaction pairs were downloaded from the TRRUST (http://www.grnpedia.org/trrust/) database. TF genes were obtained from the FANTOM 5
(http://fantom.gsc.riken.jp/5/) database. TF genes were extracted to obtain miRNA-TF interaction pairs. Next, integrating the miRNA-TF interaction pairs with TF-gene pairs, a miRNA-TF-gene network was constructed.
2.5 Construction of the protein-protein interaction network Based on the BioGRID (https://thebiogrid.org/) and STRING (https://string-db.org/) databases, proteins interacting with the DEGs were extracted and used to construct a protein-protein interaction (PPI) network. Ten well-known signaling pathways were obtained from the KEGG database, including Hedgehog, Hippo, JAK-STAT, MAPK, mTOR, NF-kappa B, Notch, PI3KAkt, TGF-beta, and Wnt signaling pathways. These pathways were mapped onto the PPI network, and shown as edges connecting the DEGs. In addition, miRNAs that target the DEGs used to construct a miRNA-target-pathway network.
2.6 Screening for significant modules The original CEL files of included datasets were downloaded from the GEO database and integrated. A differential expression analysis identified 10,159 DEGs from this integrated dataset. A weighted gene co-expression network analysis (WGCNA) was performed to identify the genes that were co-expressed both pre-SCI and post-SCI. Next, the modules closely related to differential expression were screened. The correlation between the degree (k) of the gene in the module and -log10 (p) was calculated.
Additionally, significant modules in which hub DEGs existed were selected. Functional enrichment analysis of the genes in these selected modules was performed using the Cytoscape plug-in clueGO. The DEGs that were commonly differentially expressed in all datasets were obtained from these significant modules and miRNA-gene-TF networks for the significant modules were constructed.
2.7 In-vivo validation of key molecules’ expression patterns by qRT-PCR Six C57BL/6 mice were randomly assigned to two groups using a computer-generated randomization schedule: mice in the sham control group (n = 3) were treated with only
laminectomy, while rodents in the SCI group (n = 3) were subjected to laminectomy plus spinal cord contusion. The expression level of four genes (Flna, ID3, HK2, Ywhae) and two TFs (Jun and Sp1) identified from aforementioned bioinformatics analysis was validated by qRT-PCR, as described in detail by Qin et al [16]. In each sample, three replicates per gene were tested. Statistical analysis was performed using SPSS software (version 21.0, Chicago, IL, United States). Descriptive data were summarized as mean ± SD. After testing for normality, a student’s t-test was used to compare expression levels between the two groups. P < 0.05 was considered as statistically significant.
3. Results 3.1 Procurement of datasets and identification of DEGs Three datasets (GSE5296, GSE42828, and GSE47681) were included which had the same experimental type (mRNAs expression profiling) and were based on the same platform GPL1261 (Table S1). By performing differential expression analysis on individual datasets, we identified 2,061 DEGs in the GSE5296 dataset, including 1,413 up-regulated DEGs and 648 down-regulated DEGs; 1,474 DEGs the in GSE42828 dataset, including 1,083 up-regulated DEGs and 391 down-regulated DEGs, and 1,472 DEGs in the GSE67515 dataset, including 1,082 up-regulated DEGs and 390 down-regulated DEGs (Table S2). DEGs the overlapped among these three datasets included 836 DEGs with 756 up-regulated DEGs and 80 downregulated DEGs (Table S2, Fig S4). Functional enrichment analysis identified significantly enriched signaling pathways within these 836 DEGs, and was found to include PI3K-Akt, Tolllike receptor, chemokine, and NF-kappa B signaling pathways (Fig S4).
3.2 The construction of miRNA-target network A total of 636 miRNA-DEGs interaction pairs were used to construct the miRNA-target network involved in subacute SCI (Fig 1a). As evident in Fig 1a, each DEG was targeted by multiple miRNAs; and vice versa, each miRNA regulated the expression of many DEGs. By analyzing the topological characteristics of this network, the top 30 nodes with the highest degree were identified (Table S3). The top 5 genes found to be critical in the miRNA-target regulation network in SCI included Flna, Zfp36l1, Rbm47, Id3, and Hk2.
3.3 ceRNA network and miRNA-TF-gene interaction network The ceRNA network in SCI reflecting the interactions of lncRNAs-miRNAs-mRNAs was depicted in Fig 1b, and contained 655 nodes and 10564 edges. By analyzing its topological characteristics. the top 20 nodes with highest degree were identified (Table S4). Certain lncRNAs (1700020I14Rik, Neat1, Xist, and Malat1), DEGs (Rbm47, Zfp36l1, Hk2, and Id3), and miRNAs (miR-1958, miR-185-5p, miR-882, and miR-669f-3p) (Table S4) were identified to play critical roles in the ceRNA network as they regulated the greatest number of other factors. Furthermore, cluster analysis was performed and modules in the ceRNA network were extracted. A total of 10 modules that interacted with each other were identified (Fig S4). In addition, a miRNA-TF-gene interaction network (Fig 2a) was constructed, which was found to contain 1,589 nodes, and 2,927 edges. The transcription factors (Sp1, Trp53, Jun and Rela) regulated multiple genes and miRNAs in this network, showing their central roles.
3.4 PPI network and miRNA-target-pathway regulatory network The PPI network based on the DEGs can be seen in Fig 2b. By topological analysis, the top 10 nodes were determined (Table S5). Certain DEGs (Ywhae, Jun, Ywhaz, and Traf6) were of key importance in the PPI network as these were found to regulate the expression of highest numbers of genes. In addition, the PPI network indicated that the DEGs were involved only in 7 pathways, namely the Wnt, MAPK, PI3K-Akt, TGF-beta, mTOR, NF-kappa B, and JAKSTAT signaling pathways. The DEGs involved in these 7 pathways and miRNAs targeting these DEGs were extracted and an miRNA-target-pathway regulatory network was constructed (Fig 3). It was noted that the Flna gene was targeted by various miRNAs (i.e., miR-135b-5p, miR-135a-5p, miR-146a-5p, miR-146b-5p, and miR-128-3p) is involved in the MAPK signaling pathway, and many DEGs (i.e., Col4a1, Thbs1, Ifnar2, Spp1, Fn1, Eif4ebp1, and Igf2) were found as involved in PI3K-Akt signaling pathway.
3.5 Enriched pathways and miRNA-gene-TF network of significant modules WGCNA analysis revealed five co-expressed modules (Fig 4a). Table S6 lists the significance of enrichment for each gene module respectively. Observed from Table S6 and Fig 4b, it can
be found that the significance of enrichment for the turquoise module was the highest (log10(p)= 8.457161537) amongst all modules. This also indicates that hub DEGs existed with highest frequency in the turquoise module (Fig 4b). The functional enrichment analysis of DEGs in the turquoise module showed that these DEGs were involved in many signaling pathways, including NF-kappa B, Toll-like receptor, and B cell receptor signaling pathways (Fig 5).
A miRNA-gene-TF regulation network for the turquoise module was constructed (Fig 6a). Some DEGs also acting as transcription factors are shown to play critical roles in this network by interacting with many other transcription factors and genes. For example, Jun as a DEG that acts as a TF is shown to be targeted by miR-760-5p and seen to regulate the expression of many non-gene TFs (i.e., Mef2c, Mmp9, Rarb, Pelp1, and Hdac3). As another example, Stat1 as a DEG and that acts as a TF is shown to be targeted by miR-674-5p and also seen to regulate the expression of many non-gene TFs (i.e., Rara, Trim24, Eif2ak2, and Rbf2) and many non-TF genes (Cebpb, Birc5, and Irf9). In addition, a subnetwork of miRNA-gene-TF (Fig 6b) shows that miR-325-3p, Asah1 gene, and transcription factor Sp1 can interact with each other by forming a closed regulation loop.
3.6 Verification by qRT-PCR experiments qRT-PCR experiments verified that Flna, ID3, HK2, Ywhae, Jun and Sp1 in the SCI group were all significantly upregulated in-vivo, as compared with the sham control group, as predicted in-silico (Fig 6c). This finding supports a high prediction accuracy of the bioinformatics analysis. Table S7 lists the forward and reverse sequences, optimal temperature, CG%, and product size of primers used for these experiments.
4. Discussion This study based on bioinformatics analysis depicted genetic and epigenetic mechanisms in subacute SCI by identifying four genes (Flna, Id3, Hk2, and Ywhae), four transcription factors (Sp1, Trp53, Jun, and RelA), two miRNAs (miR-16-5p, and miR-325-3p), and three lncRNAs (Neat1, Xist, Malat1) as critical molecular players. Existing research reports that support our
predictions of several genes, transcription factors, miRNAs, and lncRNAs identified in this study as involved in the pathogenesis of spinal cord injury are summarized as follows, in context of our findings.
Three genes (Flna, Id3, and HK2) were identified as key players in the ceRNA network analysis, and the Ywhae gene emerged as the hub gene with the highest degree in PPI network analysis. Flna (filamin, alpha) has been shown to regulate the proliferation of a neural progenitor [18], necessary for neurogenesis and axonogenesis [19]. Flna can also induce the onset of neural precursor migration [20], essential for cellular and functional regeneration after SCI [21]. Id3 (inhibitor of DNA binding 3; also known as Idb3; Hlh462; bHLHb25, a member of the helixloop-helix (HLH) gene family, is critical to cellular responses following SCI by mediating the proliferation of neural cells [22]. Specifically, Id3 expression is upregulated in a time-dependent manner in astrocytes, oligodendrocytes, and neural progenitor subpopulations after SCI [22]. Id3 upregulation is notable even at eight days post-SCI, suggesting its relevance to the early healing period [22]. Hk2 (hexokinase 2), a member of the hexokinase family, may bear neuroprotective effects by acting as a potent promoter of neuronal cell survival [23]. Another member of the hexokinase family (Hk3) is upregulated in infiltrating leukocytes, activated microglia/macrophages, and astrocytes after SCI, and is suggested to be involved in its pathologic process by promoting glucose metabolism [24]. The expression pattern of Hk2 in SCI, however, has not been investigated so far. Ywhae (tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, epsilon polypeptide may play multiple roles in neurological diseases by regulating the neuronal morphogenesis [25] and neuronal migration [26], as well as mediating neurogenesis and differentiation of neuronal progenitor cells [27]. Again, its expression pattern and regulatory roles in SCI have not yet been investigated.
Four transcription factors (Sp1, Trp53, Jun, RelA) were identified as critical roles in the pathogenesis of subacute SCI. Sp1 (trans-acting transcription factor 1) can act as a scaffolding protein by recruiting three SCI-inducible transcription factors (ATF3, c-Jun, and STAT3), necessary in promoting nerve regeneration, and supporting synergistic transactivation activity [28]. Sp1-mediated mechanism may be critical to the intrinsic regeneration response after nerve
injury via its action on these TFs [28]. Trp53 (transformation related protein 53) encodes tumor protein p53, which mediates neuronal cell apoptosis [29]. p53 is dysregulated in neurons, oligodendrocytes, and astrocytes after SCI [30]. Jun (jun proto-oncogene) is considered a permanent marker for neurons and is potential mediator of neuronal cell death, survival and regeneration [31]. In neuronal injury, c-Jun can activate the repair program in Schwann cells and induce the generation of a repair cell [32]. Since clinical trials [33, 34] involving the transplantation of Schwann cells in subacute SCI patients have been performed with successful outcomes, transplantation of Schwann cells with activated c-Jun transcription factor could improve treatment outcomes. Another identified transcription factor RelA [v-rel reticuloendotheliosis viral oncogene homolog A (avian); also known as p65; p65 NFκB; p65 NF-kappa B], is an NF-κB (nuclear factor kappa B) family member and a subunit of the NFκB transcription factor complex [35]. SCI can induce the activation of RelA and NF-κB, thereby promoting the apoptosis of neuronal soma [36] and vice versa, neuron-specific inhibition of RelA can promote neuronal survival under ischemic conditions [37]. In theory, this finding suggests a potential application of RelA inhibition in controlling ischemic cell death after SCI.
Multiple miRNAs; miR-1958, miR-185-5p, miR-882, miR-669f-3p, and miR-669b-3p were identified as key players in the ceRNA network analysis. By searching the gene targets of these miRNAs using the TargetScan database [38], we found these miRNAs can target some genes related to neuropathic pain and regenerative response involved in spinal cord injury. As the first example, miR-1958 can target genes that encode the receptor of neuropeptides, including Npffr1 (neuropeptide FF receptor 1), Npy2r (neuropeptide Y receptor Y2), and Npbwr1 (neuropeptides B/W receptor 1), and Npy6r (neuropeptide Y receptor Y6). As neuropeptide has been shown to be related to SCI-induced neuropathic pain [39], targeting the miR1958/neuropeptide signalling axis could be used for prevention and treatment of chronic pain after spinal cord injury. In case of miR-185-5p, it was found to target gene Nptx1 (neuronal pentraxin 1), which is known to encode a protein with exclusive CNS (central nervous system) expression and has been related to neuropathic pain. The overexpression of Nptx1 has been shown in the brainstem rostral ventromedial medulla (RVM) after nerve injury; and inhibiting the expression of Nptx1 can alleviate the hypersensitivity of neuropathic pain [40]. Based on
this evidence, miR-185-5p/Nptx1 could comprise a significant signalling axis for treating persistent pain resulting from spinal cord injury. In case of miR-882, this miRNA is not confidently annotated in the database; however, it shares a seed with miR-185-5p and thereby has the same predicted gene targets with miR-185-5p, suggesting similar mechanisms at play. Considering, miR-669f-3p and miR-669b-3p, these two miRNAs can target neurotrophic factors (NFs)-related genes, such as Bdnf (brain derived neurotrophic factor), Ndnf (neuronderived neurotrophic factor), as well as Gfra1 and Gfra2 (glial cell line derived neurotrophic factor family receptor alpha 1 and 2). Neurotrophic factors have been shown to promote protection and regeneration of neuron during the pathogenesis of SCI [41], thus miR-882/NFs might be used an effective therapy for improving functional recovery after SCI. In addition to the miRNAs involved in ceRNA network, another miRNA (miR-325-3p) was found to be significant due to its involvement in a closed regulation loop formed by miR-325-3p, Asah1 gene, and transcription factor-Sp1. MiR-325-3p can target genes playing roles in the protection, growth, differentiation, and regeneration of neurons; these genes include Nrep (neuronal regeneration related protein), Negr1 (neuronal growth regulator 1), Neurod2 (neurogenic differentiation 2), and Adnp (activity-dependent neuroprotective protein). This indicates that the miR-325-3p/Asah1/Sp1 axis may be a novel therapeutic target in SCI. Although the miRNAs predicted in our study can target genes involved in SCI, the expression patterns and exact roles of these miRNAs in SCI haven’t been investigated yet and therefore need to be elucidated in the future research.
Six lncRNAs; 1700020I14Rik, Neat1, Xist, Malat1, Tmem250-ps, Gm26917, were identified to be involved in the ceRNA network analysis. Among these Neat1, Xist, Malat1 have been investigated in context of SCI pathogenesis. lncRNA-Neat1 (nuclear paraspeckle assembly transcript 1) is upregulated in chronic constriction type SCI [42]. The overexpression of Neat1 can induce neuropathic pain behaviors, including mechanical and thermal hyperalgesia, by regulating the miR-381/HMGB1 axis [42]. Additionally, the upregulation of Neat1 can promote the development of neuroinflammation by upregulating the expression of proinflammatory cytokines in spinal cord injury [42]. lncRNA-Xist (X-inactive specific transcript) was found as one of the most significantly upregulated lncRNAs in SCI [43]. The overexpression of Xist can
contribute to neuronal apoptosis by downregulating the PI3K/AKT signaling pathway, implying that gene-targeted therapy consisting of Xist knockdown may elicit a protective effect on neurofunctional repair [43]. Additionally, the up-regulation of lncRNA-Malat1 (metastasisassociated lung adenocarcinoma transcript 1) can contribute to the neuroinflammatory response in SCI by the activation of the miR-199b/IKKβ/NF-κB signaling pathway and consequent proinflammatory cytokine expression [44]. The expression pattern and functions lncRNAs 1700020I14Rik, Tmem250-ps, Gm26917 in SCI pathogenesis are, however, not yet investigated.
The strengths and limitations of the present study need to be acknowledged. The greatest strength is that multiple possible genetic and epigenetic mechanisms during the pathological process of subacute SCI were examined, including DEGs, miRNAs, lncRNAs, transcription factors, and signaling pathways. Significant entities and biological processes were identified through a series of comprehensive bioinformatic analyses, including differential expression analysis, PPI network, functional enrichment analysis, and multiple co-expression networks. Another strength is that the expression patterns of the in-silico determined molecules, including genes and transcription factors were validated by qRT-PCR in an in-vivo experiment. The experimental results were consistent with the in-silico findings, suggesting high prediction accuracy of the bioinformatics approaches adopted. A limitation of this study is the inclusion of a single time point, namely 7 days post-SCI, was investigated, considering the objective to focus on the subacute stage of SCI. While multiple time points have been included previous bioinformatic addressing SCI [45-48], there is need to perform integrated bioinformatics analysis of data from sequential stages of SCI, in order to clarify the temporal events that underpin SCI neuropathology.
The present study provides several experimental directions for future research. Firstly, the genetic and epigenetic mechanisms revealed as critical to SCI could serve as therapeutic targets and be incorporated in combinational therapeutic strategies that include cell-based therapies (e.g., embryonic stem cells, neural stem cells, bone marrow stromal cells transplantation) along with transgenic modulation. More specifically, transplanting stem cells alone was found as
unable to fully repair injured spinal cords [49]. In contrast, genetic modification of stem cells before transplantation can improve outcomes, as transfected genes may aid in inhibiting apoptosis, preventing neuronal damage, and promoting neural regeneration [50-55]. Nevertheless, the biological entities identified in this study are based on computational prediction, and their translational significance to stem cell therapy in SCI needs to be investigated. Secondly, despite the genetic similarity between humans and rats [56], these findings should be verified using human neural cells in models and in vivo data from human SCI. Overall, the current findings enable improved understanding of genetic and epigenetic mechanisms that underpin subacute SCI that can be of translational value.
5. Conclusion Bioinformatic analysis and experimental validation used unveiled multiple key genetic and epigenetic mechanisms involved in the subacute stage of SCI. Four genes (Flna, ID3, HK2, and Ywhae), four transcription factors (Sp1, Trp53, Jun, and RelA), two miRNAs (miR-16-5p, and miR-325-3p), and three lncRNAs (Neat1, Xist, Malat1) emerged as critical in the pathogenesis of subacute stage of SCI.
Author contributions NW, BL and LR conceptualized the presented idea; NW, SL, YY, and LH procured data from datasets; NW, YC, ZT, YJ, and YW (Yufu Wang) conducted the bioinformatic analysis; NW, SL, YY, LH, and YC contributed to the interpretation of the results; MP and YW (Yang Wang) performed PCR experiments for validation of selected biomarkers; NW wrote the original draft; BL and LR reviewed the draft and did editing and proofreading of this manuscript; BL and LR acquired the funding and supervised this project.
Acknowledgements This work was supported by the following grants: (1) The National Key Research and Development Program of China (2017YFA0105400); (2) The National Natural Science Foundation of China (31470949, 81772349，81772398，81472122);
(3) The Key Research
and Development Program of Guangdong Province (20180236); (4) The Clinical Innovation
Research Program of Guangzhou Regenerative Medicine and Health Guangdong Laboratory (2018GZR0201006); (5) The Guangzhou Science and Technology Project (201704020221, 201707010115); (6) The Guangdong Natural Science Foundation (2017A030313594); (7) The Medical Scientific Research Foundation of Guangdong Province (A2018547).
References [1] Silva NA, Sousa N, Reis RL, Salgado AJ. From basics to clinical: a comprehensive review on
spinal
cord
injury.
Prog
Neurobiol
2014;
114:
25-57.
DOI:
10.1016/j.pneurobio.2013.11.002. [2] Mazwi NL, Adeletti K and Hirschberg RE. Traumatic Spinal Cord Injury: Recovery, Rehabilitation, and Prognosis. Curr Trauma Rep 2015; 1(3): 182-192. DOI: 10.1007/s40719-015-0023-x. [3] Cheng I, Park DY, Mayle RE, Githens M, Smith RL, Park HY, Hu SS, Alamin TF, Wood KB, Kharazi AI. Does timing of transplantation of neural stem cells following spinal cord injury affect outcomes in an animal model? J Spine Surg 2017; 3(4): 567-571. DOI: 10.21037/jss.2017.10.06. [4] Nishimura S, Yasuda A, Iwai H, Takano M, Kobayashi Y, Nori S, Tsuji O, Fujiyoshi K, Ebise H, Toyama Y, Okano H, Nakamura M. Time-dependent changes in the microenvironment of injured spinal cord affects the therapeutic potential of neural stem cell transplantation for spinal cord injury. Mol Brain 2013; 6: 3. DOI: 10.1186/1756-6606-
6-3. [5] Rowland JW, Hawryluk GW, Kwon B, Fehlings MG. Current status of acute spinal cord injury pathophysiology and emerging therapies: promise on the horizon. Neurosurg Focus 2008; 25(5): E2. DOI: 10.3171/FOC.2008.25.11.E2. [6] Finelli MJ, Wong JK, Zou H. Epigenetic regulation of sensory axon regeneration after spinal cord injury. J Neurosci 2013; 33(50): 19664-76. DOI: 10.1523/JNEUROSCI.058913.2013. [7] Wintzer M, Mladinic M, Lazarevic D, Casseler C, Cattaneo A, Nicholls J. Strategies for identifying genes that play a role in spinal cord regeneration. J Anat 2004; 204(1): 3-11. DOI: 10.1111/j.1469-7580.2004.00258.x. [8] Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone
of
a
hidden
RNA
language?
Cell
2011;
146
(3):
353-358.
DOI:
10.1016/j.cell.2011.07.014. [9] Strickland ER, Hook MA, Balaraman S, Huie JR, Grau JW, Miranda RC. MicroRNA dysregulation following spinal cord contusion: implications for neural plasticity and repair. Neuroscience 2011; 186: 146-60. DOI: 10.1016/j.neuroscience.2011.03.063. [10]Nakanishi K, Nakasa T, Tanaka N, Ishikawa M, Yamada K, Yamasaki K, Kamei N, Izumi B, Adachi N, Miyaki S, Asahara H, Ochi M. Responses of microRNAs 124a and 223 following spinal cord injury in mice. Spinal Cord 2010; 48(3): 192-6. DOI: 10.1038/sc.2009.89. [11]Yunta M, Nieto-Díaz M, Esteban FJ, Caballero-López M, Navarro-Ruíz R, Reigada D, Pita-Thomas DW, del Águila A, Muñoz-Galdeano T, Maza RM. MicroRNA dysregulation in the spinal cord following traumatic injury. PLoS One 2012; 7(4): e34534. DOI: 10.1371/journal.pone.0034534. [12]Nieto-Diaz M, Esteban FJ, Reigada D, Muñoz-Galdeano T, Yunta M, Caballero-López M, Navarro-Ruiz R, Del Águila A, Maza RM. MicroRNA dysregulation in spinal cord injury: causes, consequences and therapeutics. 2014; 8: 53. Front Cell Neurosci 2014; 8: 53. DOI: 10.3389/fncel.2014.00053. [13]Shi Z, Pan B, Feng S. The emerging role of long non‐coding RNA in spinal cord injury. J Cell Mol Med 2018; 22(4): 2055-2061. DOI: 10.1111/jcmm.13515.
[14]Ding Y, Song Z and Liu J. Aberrant LncRNA expression profile in a contusion spinal cord injury mouse model. Biomed Res Int 2016; 2016: 9249401. DOI: 10.1155/2016/9249401. [15]Duran RC, Yan H, Zheng Y, Huang X, Grill R, Kim DH, Cao Q, Wu JQ. The systematic analysis of coding and long non-coding RNAs in the sub-chronic and chronic stages of spinal cord injury. Sci Rep 2017; 7: 41008. DOI: 10.1038/srep41008. [16]Wang J, Hu B, Cao F, Sun S, Zhang Y, Zhu Q. Down regulation of lncSCIR1 after spinal cord
contusion
injury
in
rat.
Brain
Res
2015;
1624:
314-320.
DOI:
10.1016/j.brainres.2015.07.052. [17]Wang L, Wang B, Liu J, Quan Z. Construction and analysis of a spinal cord injury competitive endogenous RNA network based on the expression data of long noncoding, micro and messenger RNAs. Mol Med Rep 2019; 19(4): 3021-3034. DOI: 10.3892/mmr.2019.9979. [18]Lian G, Lu J, Hu J, Zhang J, Cross SH, Ferland RJ, Sheen VL. Filamin a regulates neural progenitor proliferation and cortical size through Wee1-dependent Cdk1 phosphorylation. J Neurosci 2012; 32(22): 7672-84. DOI: 10.1523/JNEUROSCI.0894-12.2012. [19]Hui SP, Nag TC and Ghosh S. Characterization of proliferating neural progenitors after spinal cord injury in adult zebrafish. PLoS One 2015; 10(12): e0143595. DOI: 10.1371/journal.pone.0143595. [20]Qin C, Liu CB, Yang DG, Gao F, Zhang X, Zhang C, Du LJ, Yang ML, Li JJ. Circular RNA expression alteration and bioinformatics analysis in rats after traumatic spinal cord injury. Front Mol Neurosci 2019; 11: 497. DOI: 10.3389/fnmol.2018.00497. [21]Vandestadt C, Vanwalleghem GC, Castillo HA, Li M, Schulze K, Khabooshan M, Don E, Anko M-L, Scott EK, Kaslin J. Early migration of precursor neurons initiates cellular and functional regeneration after spinal cord injury in zebrafish. bioRxiv 2019; 539940. DOI: https://doi.org/10.1101/539940. [22]Tzeng SF, Bresnahan JC, Beattie MS, de Vellis J. Upregulation of the HLH Id gene family in neural progenitors and glial cells of the rat spinal cord following contusion injury. J Neurosci Res 2001; 66(6): 1161-72. DOI: 10.1002/jnr.10089. [23]Gimenez-Cassina A, Lim F, Cerrato T, Palomo GM, Diaz-Nido J. Mitochondrial hexokinase II promotes neuronal survival and acts downstream of glycogen synthase
kinase-3. J Biol Chem 2009; 284(5): 3001-11. DOI: 10.1074/jbc.M808698200. [24]Lin YH, Wu Y, Wang Y, Yao ZF, Tang J, Wang R, Shen L, Ding SQ, Hu JG, Lü HZ. Spatio-temporal expression of Hexokinase-3 in the injured female rat spinal cords. Neurochem Int 2018; 113: 23-33. DOI: 10.1016/j.neuint.2017.11.015. [25]Cornell B, Wachi T, Zhukarev V, Toyo-oka KJHmg. Regulation of neuronal morphogenesis by 14-3-3epsilon (Ywhae) via the microtubule binding protein, doublecortin. Hum Mol Genet 2016; 25(20): 4405-4418. DOI: 10.1093/hmg/ddw270. [26]Toyo-oka K, Shionoya A, Gambello MJ, Cardoso C, Leventer R, Ward HL, Ayala R, Tsai LH, Dobyns W, Ledbetter D, Hirotsune S, Wynshaw-Boris A. 14-3-3ε is important for neuronal migration by binding to NUDEL: a molecular explanation for Miller–Dieker syndrome. Nat Genet 2003; 34(3): 274-85. DOI: 10.1038/ng1169. [27]Toyo-oka K, Wachi T, Hunt RF, Baraban SC, Taya S, Ramshaw H, Kaibuchi K, Schwarz QP, Lopez AF, Wynshaw-Boris A. 14-3-3ε and ζ regulate neurogenesis and differentiation of neuronal progenitor cells in the developing brain. J Neurosci 2014; 34(36): 12168-81. DOI: 10.1523/JNEUROSCI.2513-13.2014. [28]Kiryu-Seo S, Kato R, Ogawa T, Nakagomi S, Nagata K, Kiyama H. Neuronal injuryinducible gene is synergistically regulated by ATF3, c-Jun, and STAT3 through the interaction with Sp1 in damaged neurons. J Biol Chem 2008; 283(11): 6988-96. DOI: 10.1074/jbc.M707514200. [29]Saito N, Yamamoto T, Watanabe T, Abe Y, Kumagai T. Implications of p53 protein expression in experimental spinal cord injury. J Neurotrauma 2000; 17(2): 173-82. DOI: 10.1089/neu.2000.17.173. [30]Kotipatruni RR, Dasari VR, Veeravalli KK, Dinh DH, Fassett D, Rao JS. p53-and Baxmediated apoptosis in injured rat spinal cord. Neurochem Res 2011; 36(11): 2063-74. DOI: 10.1007/s11064-011-0530-2. [31]Herdegen T, Skene P, Bähr M. The c-Jun transcription factor–bipotential mediator of neuronal death, survival and regeneration. Trends Neurosci 1997; 20(5): 227-31. [32]Arthur-Farraj PJ, Latouche M, Wilton DK, Quintes S, Chabrol E, Banerjee A, Woodhoo A, Jenkins B, Rahman M, Turmaine M, Wicher GK, Mitter R, Greensmith L, Behrens A, Raivich G, Mirsky R, Jessen KR. c-Jun reprograms Schwann cells of injured nerves to
generate a repair cell essential for regeneration. Neuron 2012; 75(4): 633-47. DOI: 10.1016/j.neuron.2012.06.021. [33]Anderson KD, Guest JD, Dietrich WD, Bartlett Bunge M, Curiel R, Dididze M, Green BA, Khan A, Pearse DD, Saraf-Lavi E, Widerström-Noga E, Wood P, Levi AD. Safety of autologous human Schwann cell transplantation in subacute thoracic spinal cord injury. J Neurotrauma 2017; 34(21): 2950-2963. DOI: 10.1089/neu.2016.4895. [34]Guest J, Santamaria AJ, Benavides FD. Clinical translation of autologous Schwann cell transplantation for the treatment of spinal cord injury. Curr Opin Organ Transplant 2013; 18(6): 682-9. DOI: 10.1097/MOT.0000000000000026. [35]Gilmore TD. The Rel/NF-κB signal transduction pathway: introduction. Oncogene 1999; 18(49): 6842-4. DOI: 10.1038/sj.onc.1203237. [36]Bethea JR, Castro M, Keane RW, Lee TT, Dietrich WD, Yezierski RP. Traumatic spinal cord injury induces nuclear factor-κB activation. J Neurosci 1998; 18(9): 3251-60. [37]Engelmann C, Weih F, Haenold R. Role of nuclear factor kappa B in central nervous system regeneration. Neural Regen Res 2014; 9(7): 707-11. DOI: 10.4103/16735374.131572. [38]Lewis BP, Shih IH, Jones-Rhoades MW, Bartel DP, Burge CB. Prediction of mammalian microRNA targets. Cell 2003; 115(7): 787-798. DOI: 10.1016/s0092-8674(03)01018-3. [39]Coronel MF, Villar MJ, Brumovsky PR, Gonzalez SL. Spinal neuropeptide expression and neuropathic behavior in the acute and chronic phases after spinal cord injury: Effects of progesterone administration. Peptides 2017; 88:189-195. DOI: 10.1016/j.peptides. 2017.01.001. [40]Zapata A, Pontis S, Schepers RJ, Wang R, Oh E, Stein A, Backman CM, Worley P, Enguita M, Abad MA, Trullas R, Shippenberg TS. Alleviation of neuropathic pain hypersensitivity by inhibiting neuronal pentraxin 1 in the rostral ventromedial medulla. J Neurosci 2012; 32(36):12431-12436. DOI: 10.1523/JNEUROSCI.2730-12.2012. [41]Hodgetts SI, Harvey AR. Neurotrophic Factors Used to Treat Spinal Cord Injury. Vitam Horm 2017; 104:405-457. DOI: 10.1016/bs.vh.2016.11.007. [42]Xia LX, Ke C, Lu JM. NEAT1 contributes to neuropathic pain development through targeting miR‐381/HMGB1 axis in CCI rat models. J Cell Physiol 2018; 233(9): 7103-
7111. DOI: 10.1002/jcp.26526. [43]Gu S, Xie R, Liu X, Shou J, Gu W, Che X. Long coding RNA XIST contributes to neuronal apoptosis through the downregulation of AKT phosphorylation and is negatively regulated by miR-494 in rat spinal cord injury. Int J Mol Sci 2017; 18(4). pii: E732. DOI: 10.3390/ijms18040732. [44]Zhou HJ, Wang LQ, Wang DB, Yu JB, Zhu Y, Xu QS, Zheng XJ, Zhan RY. Long noncoding RNA MALAT1 contributes to inflammatory response of microglia following spinal cord injury via the modulation of a miR-199b/IKKβ/NF-κB signaling pathway. Am J Physiol Cell Physiol 2018; 315(1): C52-C61. DOI: 10.1152/ajpcell.00278.2017. [45]Wang W, Liu R, Xu Z, Niu X, Mao Z, Meng Q, Cao X. Further insight into molecular mechanism underlying thoracic spinal cord injury using bioinformatics methods. Mol Med Rep 2015; 12(6): 7851-8. DOI: 10.3892/mmr.2015.4442. [46]Yang Z, Lv Q, Wang Z, Dong X, Yang R, Zhao W. Identification of crucial genes associated with rat traumatic spinal cord injury. Mol Med Rep 2017; 15(4): 1997-2006. DOI: 10.3892/mmr.2017.6267. [47]Zhu Z, Shen Q, Zhu L, Wei X. Identification of pivotal genes and pathways for spinal cord injury via bioinformatics analysis. Mol Med Rep 2017; 16(4): 3929-3937. DOI: 10.3892/mmr.2017.7060. [48]Liu Y, Han N, Li Q, Li Z. Bioinformatics analysis of microRNA time-course expression in brown rat (Rattus norvegicus): spinal cord injury self-repair. Spine (Phila Pa 1976) 2016; 41(2): 97-103. DOI: 10.1097/BRS.0000000000001323. [49]Hasnan N, Ektas N, Tanhoffer AI, Tanhoffer R, Fornusek C, Middleton JW, Husain R, Davis GM. Exercise responses during functional electrical stimulation cycling in individuals with spinal cord injury. Med Sci Sports Exerc 2013; 45(6): 1131-8. DOI: 10.1249/MSS.0b013e3182805d5a. [50]Tang X, Cai PQ, Lin YQ, Oudega M, Blits B, Xu L, Yang YK, Zhou TH. Genetic engineering neural stem cell modified by lentivirus for repair of spinal cord injury in rats. Chin Med Sci J 2006;21(2):120-4. [51]Wu MF, Zhang SQ, Gu R, Liu JB, Li Y, Zhu QS. Transplantation of erythropoietin genemodifed neural stem cells improves the repair of injured spinal cord. Neural Regen Res
2015;10(9):1483-90. DOI: 10.4103/1673-5374.165521. [52]Zhang L, Gu S, Zhao C, Wen T. Combined treatment of neurotrophin-3 gene and neural stem cells is propitious to functional recovery after spinal cord injury. Cell Transplant 2007;16(5):475-81. [53]Zhang C, Tu F, Zhang JY, Shen L. E-cadherin-transfected neural stem cells transplantation for spinal cord injury in rats. J Huazhong Univ Sci Technolog Med Sci 2014;34(4):554558. DOI: 10.1007/s11596-014-1314-0. [54]Cui X, Chen L, Ren Y, Ji Y, Liu W, Liu J, Yan Q, Cheng L, Sun YE. Genetic modification of mesenchymal stem cells in spinal cord injury repair strategies. Biosci Trends 2013;7(5):202-8. [55]Song JL, Zheng W, Chen W, Qian Y, Ouyang YM, Fan CY. Lentivirus-mediated microRNA-124 gene-modified bone marrow mesenchymal stem cell transplantation promotes the repair of spinal cord injury in rats. Exp Mol Med 2017;49(5):e332. DOI: 10.1038/emm.2017.48. [56]Emes RD, Goodstadt L, Winter EE, Ponting CP. Comparison of the genomes of human and mouse lays the foundation of genome zoology. Hum Mol Genet 2003;12(7):701-9. DOI: 10.1093/hmg/ddg078.
Fig. 1. (a) The miRNA-target network in SCI. The network contains 496 nodes and 636 edges. The purple triangles represent miRNAs, red circles represent up-regulated DEGs, and
green circles represent down-regulated DEGs. The nodes with bigger scale represent higher degree. (b) The ceRNA network in SCI. The purple triangles represent miRNAs; red circles: up-regulated DEGs; green circles: down-regulated DEGs; dark green squares: lncRNA. Fig. 2. (a) The miRNA-TF-gene interaction network. Triangles represent miRNAs, diamonds represent TF, red circles represent up-regulated genes, and green circles represent downregulated genes. (b) The PPI network of DEGs. The PPI network contains 2,614 nodes and 9,505 edges. The nodes represent DEG, and larger nodes represent DEG with a higher degree. The edges represent different pathways. Fig. 3. The miRNA-target-pathway network. Fig. 4. (a) Clusters of co-expressed modules determined by WGCNA analysis. (b) Significantly enriched DEGs in different modules. The horizontal axis represents the modules, and the vertical axis represents the significance of enrichment. Fig. 5. Functional enrichment analysis of genes in the turquoise module. Fig. 6. (a) miRNA-gene-TF network of the turquoise module; (b) The subnetwork of miRNAgene-TF in the turquoise module consisting of closed regulation loops; (c) Relative expression fold change of four selected genes (Flna, ID3, HK2, Ywhae) and two selected TFs (Jun and Sp1) following SCI. ∗P < 0.05.
S0378-1119(19)30830-3 https://doi.org/10.1016/j.gene.2019.144171 GENE 144171
To appear in:
Gene Gene
Received Date: Revised Date: Accepted Date:
8 June 2019 13 October 2019 15 October 2019
Please cite this article as: N. Wang, L. He, Y. Yang, S. Li, Y. Chen, Z. Tian, Y. Ji, Y. Wang, M. Pang, Y. Wang, B. Liu, L. Rong, Integrated analysis of competing endogenous RNA (ceRNA) networks in subacute stage of spinal cord injury, Gene Gene (2019), doi: https://doi.org/10.1016/j.gene.2019.144171
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 Elsevier B.V. All rights reserved.
Integrated analysis of competing endogenous RNA (ceRNA) networks in subacute stage of spinal cord injury Nanxiang Wanga, Lei Hea, Yang Yanga, Simin Lib, Yuyong Chena, Zhenming Tiana, Ye Jic, Yufu Wangc, Mao Panga, Yang Wanga, Bin Liua,1,*, Limin Ronga,1,* aDepartment
of Spine Surgery, The Third Affiliated Hospital of Sun Yat-Sen University, No.
600 Tianhe Road, Tianhe District, Guangzhou, Guangdong Province, People’s Republic of China bDepartment
of Cariology, Endodontology and Periodontology, University Leipzig, Liebigstr.
12, 04103, Leipzig, Germany cDepartment
of Orthopaedics, The Second Affiliated Hospital of Harbin Medical University,
No. 246 Xuefu Road, Nangang District, Harbin, Heilongjiang Province, People’s Republic of China 1These
authors contributed equally as the corresponding author.
*Corresponding
author.
(1) Limin Rong Email: [email protected] (2) Bin Liu E-mail: [email protected]
Highlights Integrated analysis of competing endogenous RNA networks involved in the subacute stage of spinal cord injury were lacking in the literature. Integrated bioinformatic analysis identified specific genes (Flna, ID3, HK2), transcription factors (Sp1, Trp53, Jun, Rela), lnRNAs (Neat1, Xist, Malat1), miRNAs (miR-16-5p, miR-1958, miR-185-5p) as involved in the pathogenesis of spinal cord injury at subacute stage. These genetic and epigenetic mechanisms could have potential for translation in diagnostic and therapeutic research concerning spinal cord injury.
Abstract This study aims to investigate the genetic and epigenetic mechanisms involved in the pathogenesis of subacute stage of spinal cord injury (SCI). Gene-expression datasets associated with SCI were downloaded from the Gene Expression Omnibus (GEO) database, and differential expression analyses were performed in order to identify differentially expressed genes (DEGs). Multiple network types were constructed and analyzed , including proteinprotein-interaction (PPI) network, miRNA-target network, lncRNA-associated competing endogenous RNA (ceRNA) network, and miRNA-transcription factor (TF)-target network. Cluster analyses were performed to identify significant modules. To verify the prediction accuracy of the in-silico identified molecules, qRT-PCR experiments were conducted. The results depicted the Ywhae gene as the hub gene with the highest degree in the PPI network. The ceRNA network identified specific genes (Flna, ID3, and HK2), miRNAs (miR-16-5p, miR-1958, and miR-185-5p), and lncRNAs (Neat1, Xist, and Malat1) as playing critical regulating roles in the pathological mechanisms of SCI. The miRNA-TF-gene interaction network identified four important TFs (Sp1, Trp53, Jun, and Rela). The miRNA-gene-TF interaction loops from the significant modules indicated that miR-325-3p can interact with the Asah1 gene and TF-Sp1 by forming a closed loop. The qRT-PCR experiments verified four selected genes (Flna, ID3, HK2, and Ywhae) and two selected TFs (Jun, and Sp1) as significantly up-regulated following SCI. The results indicated that four genes (Flna, ID3, HK2, and Ywhae), four transcription factors (Sp1, Trp53, Jun, and RelA), two miRNAs (miR-16-5p and miR-325-3p), and three lncRNAs (Neat1, Xist, and Malat1) are likely to be involved in the
molecular mechanisms underlying the subacute stage of SCI. These findings uncover putative pathogenic mechanisms involved in SCI and might bear translation significance for future research towards therapeutic development.
Abbreviations: ceRNA, competing endogenous RNA; SCI, spinal cord injury; GEO, Gene Expression Omnibus; DEGs, differentially expressed genes; PPI, protein-protein-interaction; TF, transcription factor; mRNAs, messenger RNAs; ncRNAs, non-coding RNAs; miRNAs, microRNAs; lncRNAs, long non-coding RNAs; GO, Gene Oncology; KEGG, Kyoto Encyclopedia of Genes and Genomes; WGCNA, weighted gene co-expression network analysis; ADEV, astrocyte derived extracellular vesicles; ESCs, embryonic stem cells; NSCs, neural stem cells; BMSCs, bone marrow stromal cells.
1. Introduction Spinal cord injury (SCI) is a devastating neurological disorder caused by accidents, violence, and falls, and can often lead to neurological deficits and disabilities [1]. The pathogenesis of SCI consists of three phases, including an acute phase, a subacute phase, and a chronic phase [2]. Cell-based therapy has been found to be most effective in improving functional outcomes when applied during the subacute stage as compared to other stages [3, 4]. Therefore, the subacute phase may comprise the best period for cell-based therapeutic interventions and result in improved rehabilitation [5]. As such, a molecular understanding of the mechanistic events that occur during this stage is of utmost importance. The pathophysiological events during SCI are understood to be tightly regulated by several specific genetic and epigenetic entities [6, 7], which include both protein-coding messenger RNAs (mRNAs) and non-coding RNAs
(ncRNAs). To unravel the biological processes underlying subacute SCI, the expression patterns of these molecules require clear elucidation.
Two main types of ncRNAs, MicroRNAs (miRNAs) and long ncRNAs (lncRNAs) can suppress each other and thus act as competing endogenous RNAs (ceRNAs). They also form a regulatory lncRNA-associated ceRNA network (lncRNAs-miRNAs-mRNAs) to regulate the target mRNAs of miRNAs [8]. Some microarray or sequencing studies [9-11] have identified certain miRNAs (miR-21, miR-486, miR-20) to be dysregulated and involved in SCI by regulating inflammation, cellular death, regeneration, or gliosis [12]. A recent review [13] that summarized the results of existing lncRNAs profiling studies in SCI [14-16], noted that several lncRNAs’ dysregulation plays a pivotal role in regulating inflammatory response (lncRNA LINC00305, LncRNA LINC00341) and angiogenesis (lncRNA HIF1A‐AS2, lncRNA MEG3). Although multiple gene expression studies have investigated ncRNA signatures following SCI, their results are variable, plausibly owing to heterogeneity in experimental animal models, severity, degree and location of lesions, experiment type and the platforms used. At the same time, an integrated view of the data generated in these studies would be valuable, which can be facilitated by systematic integration using bioinformatics approaches. A recent bioinformatic study [17] has performed ceRNA network analysis in SCI; however, this study did not address different phases of SCI. Considering its specific clinical significance, the present study investigated the genetic and epigenetic mechanisms implicated during the subacute stage of SCI, by constructing an integrated ceRNA network to identify critical genes, miRNAs, lncRNAs, transcription factors, and signaling pathways.
2. Materials and methods The flowchart illustrating the main steps of this analysis is shown in Fig S1. The flowchart consists of 5 steps – identification of differentially expressed genes and co-expressed genes, identification of pathways using functional enrichment analysis, screening for significant gene modules, construction of various networks (e.g., protein-protein interaction (PPI) network, miRNA-target network, lncRNA-miRNA-gene network, and miRNA-transcription factor (TF)target network), and experimental validation of significant molecules. The detailed method for
each step will be described as follows.
2.1 Acquisition of datasets Gene expression datasets from SCI in house mice (Mus musculus) were downloaded from the GEO (Gene Expression Omnibus) database (https://www.ncbi.nlm.nih.gov/geo/). Mice with induced spinal cord contusion were defined as the SCI group, while mice without SCI were defined as the sham group. For mice in both groups, spinal cord contusion was conducted at thoracic level 8 (T8) with a moderate impact applied with the IH Spinal Cord Impactor (Precision Systems and Instrumentation, USA).
2.2 Differential expression and functional enrichment analysis Differential expression analysis was performed to identify differentially expressed genes (DEGs) from each of the selected datasets by using the R package ‘limma’ and DEGs that were common to all 3 datasets were identified. Next, functional enrichment analysis was of the shared DEGs identified enriched GO (Gene Ontology) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways, performed using the online resource DAVID v6.7 (https://david.ncifcrf.gov).
2.3 Construction of miRNA-target network and ceRNA network The miRNA-target datasets for Mus musculus were downloaded from starBase v3.0 (http://starbase.sysu.edu.cn/). The interaction pairs between miRNAs and their target DEGs were obtained in order to construct a miRNA-target interaction network. Furthermore, target lncRNAs of these miRNAs were extracted from the starBase v3.0 database. The lncRNAsmiRNAs-DEGs interaction pairs were used to construct a lncRNA-associated ceRNA network. Cluster analysis of the ceRNA network was performed by using the ‘MCODE’ plugin in Cytoscape, following which, modules representing densely connected nodes were extracted.
2.4 Construction of the miRNA-TF-target network Transcription factor (TF)-gene interaction pairs were downloaded from the TRRUST (http://www.grnpedia.org/trrust/) database. TF genes were obtained from the FANTOM 5
(http://fantom.gsc.riken.jp/5/) database. TF genes were extracted to obtain miRNA-TF interaction pairs. Next, integrating the miRNA-TF interaction pairs with TF-gene pairs, a miRNA-TF-gene network was constructed.
2.5 Construction of the protein-protein interaction network Based on the BioGRID (https://thebiogrid.org/) and STRING (https://string-db.org/) databases, proteins interacting with the DEGs were extracted and used to construct a protein-protein interaction (PPI) network. Ten well-known signaling pathways were obtained from the KEGG database, including Hedgehog, Hippo, JAK-STAT, MAPK, mTOR, NF-kappa B, Notch, PI3KAkt, TGF-beta, and Wnt signaling pathways. These pathways were mapped onto the PPI network, and shown as edges connecting the DEGs. In addition, miRNAs that target the DEGs used to construct a miRNA-target-pathway network.
2.6 Screening for significant modules The original CEL files of included datasets were downloaded from the GEO database and integrated. A differential expression analysis identified 10,159 DEGs from this integrated dataset. A weighted gene co-expression network analysis (WGCNA) was performed to identify the genes that were co-expressed both pre-SCI and post-SCI. Next, the modules closely related to differential expression were screened. The correlation between the degree (k) of the gene in the module and -log10 (p) was calculated.
Additionally, significant modules in which hub DEGs existed were selected. Functional enrichment analysis of the genes in these selected modules was performed using the Cytoscape plug-in clueGO. The DEGs that were commonly differentially expressed in all datasets were obtained from these significant modules and miRNA-gene-TF networks for the significant modules were constructed.
2.7 In-vivo validation of key molecules’ expression patterns by qRT-PCR Six C57BL/6 mice were randomly assigned to two groups using a computer-generated randomization schedule: mice in the sham control group (n = 3) were treated with only
laminectomy, while rodents in the SCI group (n = 3) were subjected to laminectomy plus spinal cord contusion. The expression level of four genes (Flna, ID3, HK2, Ywhae) and two TFs (Jun and Sp1) identified from aforementioned bioinformatics analysis was validated by qRT-PCR, as described in detail by Qin et al [16]. In each sample, three replicates per gene were tested. Statistical analysis was performed using SPSS software (version 21.0, Chicago, IL, United States). Descriptive data were summarized as mean ± SD. After testing for normality, a student’s t-test was used to compare expression levels between the two groups. P < 0.05 was considered as statistically significant.
3. Results 3.1 Procurement of datasets and identification of DEGs Three datasets (GSE5296, GSE42828, and GSE47681) were included which had the same experimental type (mRNAs expression profiling) and were based on the same platform GPL1261 (Table S1). By performing differential expression analysis on individual datasets, we identified 2,061 DEGs in the GSE5296 dataset, including 1,413 up-regulated DEGs and 648 down-regulated DEGs; 1,474 DEGs the in GSE42828 dataset, including 1,083 up-regulated DEGs and 391 down-regulated DEGs, and 1,472 DEGs in the GSE67515 dataset, including 1,082 up-regulated DEGs and 390 down-regulated DEGs (Table S2). DEGs the overlapped among these three datasets included 836 DEGs with 756 up-regulated DEGs and 80 downregulated DEGs (Table S2, Fig S4). Functional enrichment analysis identified significantly enriched signaling pathways within these 836 DEGs, and was found to include PI3K-Akt, Tolllike receptor, chemokine, and NF-kappa B signaling pathways (Fig S4).
3.2 The construction of miRNA-target network A total of 636 miRNA-DEGs interaction pairs were used to construct the miRNA-target network involved in subacute SCI (Fig 1a). As evident in Fig 1a, each DEG was targeted by multiple miRNAs; and vice versa, each miRNA regulated the expression of many DEGs. By analyzing the topological characteristics of this network, the top 30 nodes with the highest degree were identified (Table S3). The top 5 genes found to be critical in the miRNA-target regulation network in SCI included Flna, Zfp36l1, Rbm47, Id3, and Hk2.
3.3 ceRNA network and miRNA-TF-gene interaction network The ceRNA network in SCI reflecting the interactions of lncRNAs-miRNAs-mRNAs was depicted in Fig 1b, and contained 655 nodes and 10564 edges. By analyzing its topological characteristics. the top 20 nodes with highest degree were identified (Table S4). Certain lncRNAs (1700020I14Rik, Neat1, Xist, and Malat1), DEGs (Rbm47, Zfp36l1, Hk2, and Id3), and miRNAs (miR-1958, miR-185-5p, miR-882, and miR-669f-3p) (Table S4) were identified to play critical roles in the ceRNA network as they regulated the greatest number of other factors. Furthermore, cluster analysis was performed and modules in the ceRNA network were extracted. A total of 10 modules that interacted with each other were identified (Fig S4). In addition, a miRNA-TF-gene interaction network (Fig 2a) was constructed, which was found to contain 1,589 nodes, and 2,927 edges. The transcription factors (Sp1, Trp53, Jun and Rela) regulated multiple genes and miRNAs in this network, showing their central roles.
3.4 PPI network and miRNA-target-pathway regulatory network The PPI network based on the DEGs can be seen in Fig 2b. By topological analysis, the top 10 nodes were determined (Table S5). Certain DEGs (Ywhae, Jun, Ywhaz, and Traf6) were of key importance in the PPI network as these were found to regulate the expression of highest numbers of genes. In addition, the PPI network indicated that the DEGs were involved only in 7 pathways, namely the Wnt, MAPK, PI3K-Akt, TGF-beta, mTOR, NF-kappa B, and JAKSTAT signaling pathways. The DEGs involved in these 7 pathways and miRNAs targeting these DEGs were extracted and an miRNA-target-pathway regulatory network was constructed (Fig 3). It was noted that the Flna gene was targeted by various miRNAs (i.e., miR-135b-5p, miR-135a-5p, miR-146a-5p, miR-146b-5p, and miR-128-3p) is involved in the MAPK signaling pathway, and many DEGs (i.e., Col4a1, Thbs1, Ifnar2, Spp1, Fn1, Eif4ebp1, and Igf2) were found as involved in PI3K-Akt signaling pathway.
3.5 Enriched pathways and miRNA-gene-TF network of significant modules WGCNA analysis revealed five co-expressed modules (Fig 4a). Table S6 lists the significance of enrichment for each gene module respectively. Observed from Table S6 and Fig 4b, it can
be found that the significance of enrichment for the turquoise module was the highest (log10(p)= 8.457161537) amongst all modules. This also indicates that hub DEGs existed with highest frequency in the turquoise module (Fig 4b). The functional enrichment analysis of DEGs in the turquoise module showed that these DEGs were involved in many signaling pathways, including NF-kappa B, Toll-like receptor, and B cell receptor signaling pathways (Fig 5).
A miRNA-gene-TF regulation network for the turquoise module was constructed (Fig 6a). Some DEGs also acting as transcription factors are shown to play critical roles in this network by interacting with many other transcription factors and genes. For example, Jun as a DEG that acts as a TF is shown to be targeted by miR-760-5p and seen to regulate the expression of many non-gene TFs (i.e., Mef2c, Mmp9, Rarb, Pelp1, and Hdac3). As another example, Stat1 as a DEG and that acts as a TF is shown to be targeted by miR-674-5p and also seen to regulate the expression of many non-gene TFs (i.e., Rara, Trim24, Eif2ak2, and Rbf2) and many non-TF genes (Cebpb, Birc5, and Irf9). In addition, a subnetwork of miRNA-gene-TF (Fig 6b) shows that miR-325-3p, Asah1 gene, and transcription factor Sp1 can interact with each other by forming a closed regulation loop.
3.6 Verification by qRT-PCR experiments qRT-PCR experiments verified that Flna, ID3, HK2, Ywhae, Jun and Sp1 in the SCI group were all significantly upregulated in-vivo, as compared with the sham control group, as predicted in-silico (Fig 6c). This finding supports a high prediction accuracy of the bioinformatics analysis. Table S7 lists the forward and reverse sequences, optimal temperature, CG%, and product size of primers used for these experiments.
4. Discussion This study based on bioinformatics analysis depicted genetic and epigenetic mechanisms in subacute SCI by identifying four genes (Flna, Id3, Hk2, and Ywhae), four transcription factors (Sp1, Trp53, Jun, and RelA), two miRNAs (miR-16-5p, and miR-325-3p), and three lncRNAs (Neat1, Xist, Malat1) as critical molecular players. Existing research reports that support our
predictions of several genes, transcription factors, miRNAs, and lncRNAs identified in this study as involved in the pathogenesis of spinal cord injury are summarized as follows, in context of our findings.
Three genes (Flna, Id3, and HK2) were identified as key players in the ceRNA network analysis, and the Ywhae gene emerged as the hub gene with the highest degree in PPI network analysis. Flna (filamin, alpha) has been shown to regulate the proliferation of a neural progenitor [18], necessary for neurogenesis and axonogenesis [19]. Flna can also induce the onset of neural precursor migration [20], essential for cellular and functional regeneration after SCI [21]. Id3 (inhibitor of DNA binding 3; also known as Idb3; Hlh462; bHLHb25, a member of the helixloop-helix (HLH) gene family, is critical to cellular responses following SCI by mediating the proliferation of neural cells [22]. Specifically, Id3 expression is upregulated in a time-dependent manner in astrocytes, oligodendrocytes, and neural progenitor subpopulations after SCI [22]. Id3 upregulation is notable even at eight days post-SCI, suggesting its relevance to the early healing period [22]. Hk2 (hexokinase 2), a member of the hexokinase family, may bear neuroprotective effects by acting as a potent promoter of neuronal cell survival [23]. Another member of the hexokinase family (Hk3) is upregulated in infiltrating leukocytes, activated microglia/macrophages, and astrocytes after SCI, and is suggested to be involved in its pathologic process by promoting glucose metabolism [24]. The expression pattern of Hk2 in SCI, however, has not been investigated so far. Ywhae (tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, epsilon polypeptide may play multiple roles in neurological diseases by regulating the neuronal morphogenesis [25] and neuronal migration [26], as well as mediating neurogenesis and differentiation of neuronal progenitor cells [27]. Again, its expression pattern and regulatory roles in SCI have not yet been investigated.
Four transcription factors (Sp1, Trp53, Jun, RelA) were identified as critical roles in the pathogenesis of subacute SCI. Sp1 (trans-acting transcription factor 1) can act as a scaffolding protein by recruiting three SCI-inducible transcription factors (ATF3, c-Jun, and STAT3), necessary in promoting nerve regeneration, and supporting synergistic transactivation activity [28]. Sp1-mediated mechanism may be critical to the intrinsic regeneration response after nerve
injury via its action on these TFs [28]. Trp53 (transformation related protein 53) encodes tumor protein p53, which mediates neuronal cell apoptosis [29]. p53 is dysregulated in neurons, oligodendrocytes, and astrocytes after SCI [30]. Jun (jun proto-oncogene) is considered a permanent marker for neurons and is potential mediator of neuronal cell death, survival and regeneration [31]. In neuronal injury, c-Jun can activate the repair program in Schwann cells and induce the generation of a repair cell [32]. Since clinical trials [33, 34] involving the transplantation of Schwann cells in subacute SCI patients have been performed with successful outcomes, transplantation of Schwann cells with activated c-Jun transcription factor could improve treatment outcomes. Another identified transcription factor RelA [v-rel reticuloendotheliosis viral oncogene homolog A (avian); also known as p65; p65 NFκB; p65 NF-kappa B], is an NF-κB (nuclear factor kappa B) family member and a subunit of the NFκB transcription factor complex [35]. SCI can induce the activation of RelA and NF-κB, thereby promoting the apoptosis of neuronal soma [36] and vice versa, neuron-specific inhibition of RelA can promote neuronal survival under ischemic conditions [37]. In theory, this finding suggests a potential application of RelA inhibition in controlling ischemic cell death after SCI.
Multiple miRNAs; miR-1958, miR-185-5p, miR-882, miR-669f-3p, and miR-669b-3p were identified as key players in the ceRNA network analysis. By searching the gene targets of these miRNAs using the TargetScan database [38], we found these miRNAs can target some genes related to neuropathic pain and regenerative response involved in spinal cord injury. As the first example, miR-1958 can target genes that encode the receptor of neuropeptides, including Npffr1 (neuropeptide FF receptor 1), Npy2r (neuropeptide Y receptor Y2), and Npbwr1 (neuropeptides B/W receptor 1), and Npy6r (neuropeptide Y receptor Y6). As neuropeptide has been shown to be related to SCI-induced neuropathic pain [39], targeting the miR1958/neuropeptide signalling axis could be used for prevention and treatment of chronic pain after spinal cord injury. In case of miR-185-5p, it was found to target gene Nptx1 (neuronal pentraxin 1), which is known to encode a protein with exclusive CNS (central nervous system) expression and has been related to neuropathic pain. The overexpression of Nptx1 has been shown in the brainstem rostral ventromedial medulla (RVM) after nerve injury; and inhibiting the expression of Nptx1 can alleviate the hypersensitivity of neuropathic pain [40]. Based on
this evidence, miR-185-5p/Nptx1 could comprise a significant signalling axis for treating persistent pain resulting from spinal cord injury. In case of miR-882, this miRNA is not confidently annotated in the database; however, it shares a seed with miR-185-5p and thereby has the same predicted gene targets with miR-185-5p, suggesting similar mechanisms at play. Considering, miR-669f-3p and miR-669b-3p, these two miRNAs can target neurotrophic factors (NFs)-related genes, such as Bdnf (brain derived neurotrophic factor), Ndnf (neuronderived neurotrophic factor), as well as Gfra1 and Gfra2 (glial cell line derived neurotrophic factor family receptor alpha 1 and 2). Neurotrophic factors have been shown to promote protection and regeneration of neuron during the pathogenesis of SCI [41], thus miR-882/NFs might be used an effective therapy for improving functional recovery after SCI. In addition to the miRNAs involved in ceRNA network, another miRNA (miR-325-3p) was found to be significant due to its involvement in a closed regulation loop formed by miR-325-3p, Asah1 gene, and transcription factor-Sp1. MiR-325-3p can target genes playing roles in the protection, growth, differentiation, and regeneration of neurons; these genes include Nrep (neuronal regeneration related protein), Negr1 (neuronal growth regulator 1), Neurod2 (neurogenic differentiation 2), and Adnp (activity-dependent neuroprotective protein). This indicates that the miR-325-3p/Asah1/Sp1 axis may be a novel therapeutic target in SCI. Although the miRNAs predicted in our study can target genes involved in SCI, the expression patterns and exact roles of these miRNAs in SCI haven’t been investigated yet and therefore need to be elucidated in the future research.
Six lncRNAs; 1700020I14Rik, Neat1, Xist, Malat1, Tmem250-ps, Gm26917, were identified to be involved in the ceRNA network analysis. Among these Neat1, Xist, Malat1 have been investigated in context of SCI pathogenesis. lncRNA-Neat1 (nuclear paraspeckle assembly transcript 1) is upregulated in chronic constriction type SCI [42]. The overexpression of Neat1 can induce neuropathic pain behaviors, including mechanical and thermal hyperalgesia, by regulating the miR-381/HMGB1 axis [42]. Additionally, the upregulation of Neat1 can promote the development of neuroinflammation by upregulating the expression of proinflammatory cytokines in spinal cord injury [42]. lncRNA-Xist (X-inactive specific transcript) was found as one of the most significantly upregulated lncRNAs in SCI [43]. The overexpression of Xist can
contribute to neuronal apoptosis by downregulating the PI3K/AKT signaling pathway, implying that gene-targeted therapy consisting of Xist knockdown may elicit a protective effect on neurofunctional repair [43]. Additionally, the up-regulation of lncRNA-Malat1 (metastasisassociated lung adenocarcinoma transcript 1) can contribute to the neuroinflammatory response in SCI by the activation of the miR-199b/IKKβ/NF-κB signaling pathway and consequent proinflammatory cytokine expression [44]. The expression pattern and functions lncRNAs 1700020I14Rik, Tmem250-ps, Gm26917 in SCI pathogenesis are, however, not yet investigated.
The strengths and limitations of the present study need to be acknowledged. The greatest strength is that multiple possible genetic and epigenetic mechanisms during the pathological process of subacute SCI were examined, including DEGs, miRNAs, lncRNAs, transcription factors, and signaling pathways. Significant entities and biological processes were identified through a series of comprehensive bioinformatic analyses, including differential expression analysis, PPI network, functional enrichment analysis, and multiple co-expression networks. Another strength is that the expression patterns of the in-silico determined molecules, including genes and transcription factors were validated by qRT-PCR in an in-vivo experiment. The experimental results were consistent with the in-silico findings, suggesting high prediction accuracy of the bioinformatics approaches adopted. A limitation of this study is the inclusion of a single time point, namely 7 days post-SCI, was investigated, considering the objective to focus on the subacute stage of SCI. While multiple time points have been included previous bioinformatic addressing SCI [45-48], there is need to perform integrated bioinformatics analysis of data from sequential stages of SCI, in order to clarify the temporal events that underpin SCI neuropathology.
The present study provides several experimental directions for future research. Firstly, the genetic and epigenetic mechanisms revealed as critical to SCI could serve as therapeutic targets and be incorporated in combinational therapeutic strategies that include cell-based therapies (e.g., embryonic stem cells, neural stem cells, bone marrow stromal cells transplantation) along with transgenic modulation. More specifically, transplanting stem cells alone was found as
unable to fully repair injured spinal cords [49]. In contrast, genetic modification of stem cells before transplantation can improve outcomes, as transfected genes may aid in inhibiting apoptosis, preventing neuronal damage, and promoting neural regeneration [50-55]. Nevertheless, the biological entities identified in this study are based on computational prediction, and their translational significance to stem cell therapy in SCI needs to be investigated. Secondly, despite the genetic similarity between humans and rats [56], these findings should be verified using human neural cells in models and in vivo data from human SCI. Overall, the current findings enable improved understanding of genetic and epigenetic mechanisms that underpin subacute SCI that can be of translational value.
5. Conclusion Bioinformatic analysis and experimental validation used unveiled multiple key genetic and epigenetic mechanisms involved in the subacute stage of SCI. Four genes (Flna, ID3, HK2, and Ywhae), four transcription factors (Sp1, Trp53, Jun, and RelA), two miRNAs (miR-16-5p, and miR-325-3p), and three lncRNAs (Neat1, Xist, Malat1) emerged as critical in the pathogenesis of subacute stage of SCI.
Author contributions NW, BL and LR conceptualized the presented idea; NW, SL, YY, and LH procured data from datasets; NW, YC, ZT, YJ, and YW (Yufu Wang) conducted the bioinformatic analysis; NW, SL, YY, LH, and YC contributed to the interpretation of the results; MP and YW (Yang Wang) performed PCR experiments for validation of selected biomarkers; NW wrote the original draft; BL and LR reviewed the draft and did editing and proofreading of this manuscript; BL and LR acquired the funding and supervised this project.
Acknowledgements This work was supported by the following grants: (1) The National Key Research and Development Program of China (2017YFA0105400); (2) The National Natural Science Foundation of China (31470949, 81772349，81772398，81472122);
(3) The Key Research
and Development Program of Guangdong Province (20180236); (4) The Clinical Innovation
Research Program of Guangzhou Regenerative Medicine and Health Guangdong Laboratory (2018GZR0201006); (5) The Guangzhou Science and Technology Project (201704020221, 201707010115); (6) The Guangdong Natural Science Foundation (2017A030313594); (7) The Medical Scientific Research Foundation of Guangdong Province (A2018547).
References [1] Silva NA, Sousa N, Reis RL, Salgado AJ. From basics to clinical: a comprehensive review on
spinal
cord
injury.
Prog
Neurobiol
2014;
114:
25-57.
DOI:
10.1016/j.pneurobio.2013.11.002. [2] Mazwi NL, Adeletti K and Hirschberg RE. Traumatic Spinal Cord Injury: Recovery, Rehabilitation, and Prognosis. Curr Trauma Rep 2015; 1(3): 182-192. DOI: 10.1007/s40719-015-0023-x. [3] Cheng I, Park DY, Mayle RE, Githens M, Smith RL, Park HY, Hu SS, Alamin TF, Wood KB, Kharazi AI. Does timing of transplantation of neural stem cells following spinal cord injury affect outcomes in an animal model? J Spine Surg 2017; 3(4): 567-571. DOI: 10.21037/jss.2017.10.06. [4] Nishimura S, Yasuda A, Iwai H, Takano M, Kobayashi Y, Nori S, Tsuji O, Fujiyoshi K, Ebise H, Toyama Y, Okano H, Nakamura M. Time-dependent changes in the microenvironment of injured spinal cord affects the therapeutic potential of neural stem cell transplantation for spinal cord injury. Mol Brain 2013; 6: 3. DOI: 10.1186/1756-6606-
6-3. [5] Rowland JW, Hawryluk GW, Kwon B, Fehlings MG. Current status of acute spinal cord injury pathophysiology and emerging therapies: promise on the horizon. Neurosurg Focus 2008; 25(5): E2. DOI: 10.3171/FOC.2008.25.11.E2. [6] Finelli MJ, Wong JK, Zou H. Epigenetic regulation of sensory axon regeneration after spinal cord injury. J Neurosci 2013; 33(50): 19664-76. DOI: 10.1523/JNEUROSCI.058913.2013. [7] Wintzer M, Mladinic M, Lazarevic D, Casseler C, Cattaneo A, Nicholls J. Strategies for identifying genes that play a role in spinal cord regeneration. J Anat 2004; 204(1): 3-11. DOI: 10.1111/j.1469-7580.2004.00258.x. [8] Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone
of
a
hidden
RNA
language?
Cell
2011;
146
(3):
353-358.
DOI:
10.1016/j.cell.2011.07.014. [9] Strickland ER, Hook MA, Balaraman S, Huie JR, Grau JW, Miranda RC. MicroRNA dysregulation following spinal cord contusion: implications for neural plasticity and repair. Neuroscience 2011; 186: 146-60. DOI: 10.1016/j.neuroscience.2011.03.063. [10]Nakanishi K, Nakasa T, Tanaka N, Ishikawa M, Yamada K, Yamasaki K, Kamei N, Izumi B, Adachi N, Miyaki S, Asahara H, Ochi M. Responses of microRNAs 124a and 223 following spinal cord injury in mice. Spinal Cord 2010; 48(3): 192-6. DOI: 10.1038/sc.2009.89. [11]Yunta M, Nieto-Díaz M, Esteban FJ, Caballero-López M, Navarro-Ruíz R, Reigada D, Pita-Thomas DW, del Águila A, Muñoz-Galdeano T, Maza RM. MicroRNA dysregulation in the spinal cord following traumatic injury. PLoS One 2012; 7(4): e34534. DOI: 10.1371/journal.pone.0034534. [12]Nieto-Diaz M, Esteban FJ, Reigada D, Muñoz-Galdeano T, Yunta M, Caballero-López M, Navarro-Ruiz R, Del Águila A, Maza RM. MicroRNA dysregulation in spinal cord injury: causes, consequences and therapeutics. 2014; 8: 53. Front Cell Neurosci 2014; 8: 53. DOI: 10.3389/fncel.2014.00053. [13]Shi Z, Pan B, Feng S. The emerging role of long non‐coding RNA in spinal cord injury. J Cell Mol Med 2018; 22(4): 2055-2061. DOI: 10.1111/jcmm.13515.
[14]Ding Y, Song Z and Liu J. Aberrant LncRNA expression profile in a contusion spinal cord injury mouse model. Biomed Res Int 2016; 2016: 9249401. DOI: 10.1155/2016/9249401. [15]Duran RC, Yan H, Zheng Y, Huang X, Grill R, Kim DH, Cao Q, Wu JQ. The systematic analysis of coding and long non-coding RNAs in the sub-chronic and chronic stages of spinal cord injury. Sci Rep 2017; 7: 41008. DOI: 10.1038/srep41008. [16]Wang J, Hu B, Cao F, Sun S, Zhang Y, Zhu Q. Down regulation of lncSCIR1 after spinal cord
contusion
injury
in
rat.
Brain
Res
2015;
1624:
314-320.
DOI:
10.1016/j.brainres.2015.07.052. [17]Wang L, Wang B, Liu J, Quan Z. Construction and analysis of a spinal cord injury competitive endogenous RNA network based on the expression data of long noncoding, micro and messenger RNAs. Mol Med Rep 2019; 19(4): 3021-3034. DOI: 10.3892/mmr.2019.9979. [18]Lian G, Lu J, Hu J, Zhang J, Cross SH, Ferland RJ, Sheen VL. Filamin a regulates neural progenitor proliferation and cortical size through Wee1-dependent Cdk1 phosphorylation. J Neurosci 2012; 32(22): 7672-84. DOI: 10.1523/JNEUROSCI.0894-12.2012. [19]Hui SP, Nag TC and Ghosh S. Characterization of proliferating neural progenitors after spinal cord injury in adult zebrafish. PLoS One 2015; 10(12): e0143595. DOI: 10.1371/journal.pone.0143595. [20]Qin C, Liu CB, Yang DG, Gao F, Zhang X, Zhang C, Du LJ, Yang ML, Li JJ. Circular RNA expression alteration and bioinformatics analysis in rats after traumatic spinal cord injury. Front Mol Neurosci 2019; 11: 497. DOI: 10.3389/fnmol.2018.00497. [21]Vandestadt C, Vanwalleghem GC, Castillo HA, Li M, Schulze K, Khabooshan M, Don E, Anko M-L, Scott EK, Kaslin J. Early migration of precursor neurons initiates cellular and functional regeneration after spinal cord injury in zebrafish. bioRxiv 2019; 539940. DOI: https://doi.org/10.1101/539940. [22]Tzeng SF, Bresnahan JC, Beattie MS, de Vellis J. Upregulation of the HLH Id gene family in neural progenitors and glial cells of the rat spinal cord following contusion injury. J Neurosci Res 2001; 66(6): 1161-72. DOI: 10.1002/jnr.10089. [23]Gimenez-Cassina A, Lim F, Cerrato T, Palomo GM, Diaz-Nido J. Mitochondrial hexokinase II promotes neuronal survival and acts downstream of glycogen synthase
kinase-3. J Biol Chem 2009; 284(5): 3001-11. DOI: 10.1074/jbc.M808698200. [24]Lin YH, Wu Y, Wang Y, Yao ZF, Tang J, Wang R, Shen L, Ding SQ, Hu JG, Lü HZ. Spatio-temporal expression of Hexokinase-3 in the injured female rat spinal cords. Neurochem Int 2018; 113: 23-33. DOI: 10.1016/j.neuint.2017.11.015. [25]Cornell B, Wachi T, Zhukarev V, Toyo-oka KJHmg. Regulation of neuronal morphogenesis by 14-3-3epsilon (Ywhae) via the microtubule binding protein, doublecortin. Hum Mol Genet 2016; 25(20): 4405-4418. DOI: 10.1093/hmg/ddw270. [26]Toyo-oka K, Shionoya A, Gambello MJ, Cardoso C, Leventer R, Ward HL, Ayala R, Tsai LH, Dobyns W, Ledbetter D, Hirotsune S, Wynshaw-Boris A. 14-3-3ε is important for neuronal migration by binding to NUDEL: a molecular explanation for Miller–Dieker syndrome. Nat Genet 2003; 34(3): 274-85. DOI: 10.1038/ng1169. [27]Toyo-oka K, Wachi T, Hunt RF, Baraban SC, Taya S, Ramshaw H, Kaibuchi K, Schwarz QP, Lopez AF, Wynshaw-Boris A. 14-3-3ε and ζ regulate neurogenesis and differentiation of neuronal progenitor cells in the developing brain. J Neurosci 2014; 34(36): 12168-81. DOI: 10.1523/JNEUROSCI.2513-13.2014. [28]Kiryu-Seo S, Kato R, Ogawa T, Nakagomi S, Nagata K, Kiyama H. Neuronal injuryinducible gene is synergistically regulated by ATF3, c-Jun, and STAT3 through the interaction with Sp1 in damaged neurons. J Biol Chem 2008; 283(11): 6988-96. DOI: 10.1074/jbc.M707514200. [29]Saito N, Yamamoto T, Watanabe T, Abe Y, Kumagai T. Implications of p53 protein expression in experimental spinal cord injury. J Neurotrauma 2000; 17(2): 173-82. DOI: 10.1089/neu.2000.17.173. [30]Kotipatruni RR, Dasari VR, Veeravalli KK, Dinh DH, Fassett D, Rao JS. p53-and Baxmediated apoptosis in injured rat spinal cord. Neurochem Res 2011; 36(11): 2063-74. DOI: 10.1007/s11064-011-0530-2. [31]Herdegen T, Skene P, Bähr M. The c-Jun transcription factor–bipotential mediator of neuronal death, survival and regeneration. Trends Neurosci 1997; 20(5): 227-31. [32]Arthur-Farraj PJ, Latouche M, Wilton DK, Quintes S, Chabrol E, Banerjee A, Woodhoo A, Jenkins B, Rahman M, Turmaine M, Wicher GK, Mitter R, Greensmith L, Behrens A, Raivich G, Mirsky R, Jessen KR. c-Jun reprograms Schwann cells of injured nerves to
generate a repair cell essential for regeneration. Neuron 2012; 75(4): 633-47. DOI: 10.1016/j.neuron.2012.06.021. [33]Anderson KD, Guest JD, Dietrich WD, Bartlett Bunge M, Curiel R, Dididze M, Green BA, Khan A, Pearse DD, Saraf-Lavi E, Widerström-Noga E, Wood P, Levi AD. Safety of autologous human Schwann cell transplantation in subacute thoracic spinal cord injury. J Neurotrauma 2017; 34(21): 2950-2963. DOI: 10.1089/neu.2016.4895. [34]Guest J, Santamaria AJ, Benavides FD. Clinical translation of autologous Schwann cell transplantation for the treatment of spinal cord injury. Curr Opin Organ Transplant 2013; 18(6): 682-9. DOI: 10.1097/MOT.0000000000000026. [35]Gilmore TD. The Rel/NF-κB signal transduction pathway: introduction. Oncogene 1999; 18(49): 6842-4. DOI: 10.1038/sj.onc.1203237. [36]Bethea JR, Castro M, Keane RW, Lee TT, Dietrich WD, Yezierski RP. Traumatic spinal cord injury induces nuclear factor-κB activation. J Neurosci 1998; 18(9): 3251-60. [37]Engelmann C, Weih F, Haenold R. Role of nuclear factor kappa B in central nervous system regeneration. Neural Regen Res 2014; 9(7): 707-11. DOI: 10.4103/16735374.131572. [38]Lewis BP, Shih IH, Jones-Rhoades MW, Bartel DP, Burge CB. Prediction of mammalian microRNA targets. Cell 2003; 115(7): 787-798. DOI: 10.1016/s0092-8674(03)01018-3. [39]Coronel MF, Villar MJ, Brumovsky PR, Gonzalez SL. Spinal neuropeptide expression and neuropathic behavior in the acute and chronic phases after spinal cord injury: Effects of progesterone administration. Peptides 2017; 88:189-195. DOI: 10.1016/j.peptides. 2017.01.001. [40]Zapata A, Pontis S, Schepers RJ, Wang R, Oh E, Stein A, Backman CM, Worley P, Enguita M, Abad MA, Trullas R, Shippenberg TS. Alleviation of neuropathic pain hypersensitivity by inhibiting neuronal pentraxin 1 in the rostral ventromedial medulla. J Neurosci 2012; 32(36):12431-12436. DOI: 10.1523/JNEUROSCI.2730-12.2012. [41]Hodgetts SI, Harvey AR. Neurotrophic Factors Used to Treat Spinal Cord Injury. Vitam Horm 2017; 104:405-457. DOI: 10.1016/bs.vh.2016.11.007. [42]Xia LX, Ke C, Lu JM. NEAT1 contributes to neuropathic pain development through targeting miR‐381/HMGB1 axis in CCI rat models. J Cell Physiol 2018; 233(9): 7103-
7111. DOI: 10.1002/jcp.26526. [43]Gu S, Xie R, Liu X, Shou J, Gu W, Che X. Long coding RNA XIST contributes to neuronal apoptosis through the downregulation of AKT phosphorylation and is negatively regulated by miR-494 in rat spinal cord injury. Int J Mol Sci 2017; 18(4). pii: E732. DOI: 10.3390/ijms18040732. [44]Zhou HJ, Wang LQ, Wang DB, Yu JB, Zhu Y, Xu QS, Zheng XJ, Zhan RY. Long noncoding RNA MALAT1 contributes to inflammatory response of microglia following spinal cord injury via the modulation of a miR-199b/IKKβ/NF-κB signaling pathway. Am J Physiol Cell Physiol 2018; 315(1): C52-C61. DOI: 10.1152/ajpcell.00278.2017. [45]Wang W, Liu R, Xu Z, Niu X, Mao Z, Meng Q, Cao X. Further insight into molecular mechanism underlying thoracic spinal cord injury using bioinformatics methods. Mol Med Rep 2015; 12(6): 7851-8. DOI: 10.3892/mmr.2015.4442. [46]Yang Z, Lv Q, Wang Z, Dong X, Yang R, Zhao W. Identification of crucial genes associated with rat traumatic spinal cord injury. Mol Med Rep 2017; 15(4): 1997-2006. DOI: 10.3892/mmr.2017.6267. [47]Zhu Z, Shen Q, Zhu L, Wei X. Identification of pivotal genes and pathways for spinal cord injury via bioinformatics analysis. Mol Med Rep 2017; 16(4): 3929-3937. DOI: 10.3892/mmr.2017.7060. [48]Liu Y, Han N, Li Q, Li Z. Bioinformatics analysis of microRNA time-course expression in brown rat (Rattus norvegicus): spinal cord injury self-repair. Spine (Phila Pa 1976) 2016; 41(2): 97-103. DOI: 10.1097/BRS.0000000000001323. [49]Hasnan N, Ektas N, Tanhoffer AI, Tanhoffer R, Fornusek C, Middleton JW, Husain R, Davis GM. Exercise responses during functional electrical stimulation cycling in individuals with spinal cord injury. Med Sci Sports Exerc 2013; 45(6): 1131-8. DOI: 10.1249/MSS.0b013e3182805d5a. [50]Tang X, Cai PQ, Lin YQ, Oudega M, Blits B, Xu L, Yang YK, Zhou TH. Genetic engineering neural stem cell modified by lentivirus for repair of spinal cord injury in rats. Chin Med Sci J 2006;21(2):120-4. [51]Wu MF, Zhang SQ, Gu R, Liu JB, Li Y, Zhu QS. Transplantation of erythropoietin genemodifed neural stem cells improves the repair of injured spinal cord. Neural Regen Res
2015;10(9):1483-90. DOI: 10.4103/1673-5374.165521. [52]Zhang L, Gu S, Zhao C, Wen T. Combined treatment of neurotrophin-3 gene and neural stem cells is propitious to functional recovery after spinal cord injury. Cell Transplant 2007;16(5):475-81. [53]Zhang C, Tu F, Zhang JY, Shen L. E-cadherin-transfected neural stem cells transplantation for spinal cord injury in rats. J Huazhong Univ Sci Technolog Med Sci 2014;34(4):554558. DOI: 10.1007/s11596-014-1314-0. [54]Cui X, Chen L, Ren Y, Ji Y, Liu W, Liu J, Yan Q, Cheng L, Sun YE. Genetic modification of mesenchymal stem cells in spinal cord injury repair strategies. Biosci Trends 2013;7(5):202-8. [55]Song JL, Zheng W, Chen W, Qian Y, Ouyang YM, Fan CY. Lentivirus-mediated microRNA-124 gene-modified bone marrow mesenchymal stem cell transplantation promotes the repair of spinal cord injury in rats. Exp Mol Med 2017;49(5):e332. DOI: 10.1038/emm.2017.48. [56]Emes RD, Goodstadt L, Winter EE, Ponting CP. Comparison of the genomes of human and mouse lays the foundation of genome zoology. Hum Mol Genet 2003;12(7):701-9. DOI: 10.1093/hmg/ddg078.
Fig. 1. (a) The miRNA-target network in SCI. The network contains 496 nodes and 636 edges. The purple triangles represent miRNAs, red circles represent up-regulated DEGs, and
green circles represent down-regulated DEGs. The nodes with bigger scale represent higher degree. (b) The ceRNA network in SCI. The purple triangles represent miRNAs; red circles: up-regulated DEGs; green circles: down-regulated DEGs; dark green squares: lncRNA. Fig. 2. (a) The miRNA-TF-gene interaction network. Triangles represent miRNAs, diamonds represent TF, red circles represent up-regulated genes, and green circles represent downregulated genes. (b) The PPI network of DEGs. The PPI network contains 2,614 nodes and 9,505 edges. The nodes represent DEG, and larger nodes represent DEG with a higher degree. The edges represent different pathways. Fig. 3. The miRNA-target-pathway network. Fig. 4. (a) Clusters of co-expressed modules determined by WGCNA analysis. (b) Significantly enriched DEGs in different modules. The horizontal axis represents the modules, and the vertical axis represents the significance of enrichment. Fig. 5. Functional enrichment analysis of genes in the turquoise module. Fig. 6. (a) miRNA-gene-TF network of the turquoise module; (b) The subnetwork of miRNAgene-TF in the turquoise module consisting of closed regulation loops; (c) Relative expression fold change of four selected genes (Flna, ID3, HK2, Ywhae) and two selected TFs (Jun and Sp1) following SCI. ∗P < 0.05.