Identifying the Biomarkers of Spinal Cord Injury and the effects of Neurotrophin-3 Based on MicroRNA and mRNA Signature

Background: To gain a better understanding of the molecular mechanisms of spinal cord injury and the effects of Neurotrophin-3, differentially expressed microRNAs (DEmiRNAs) and genes (DEGs) were analyzed. Methods: The miRNA transcription profile of GSE82195 and the mRNA transcription profile of GSE82196 including dorsal root ganglion (DRG) tissue samples of normal adult rat (none, n=6), 10 weeks post-pyramidotomy, intramuscular AAV-1 GFP (injury, n=6) and 10 weeks post-pyramidotomy, intramuscular AAV-1 prepro-neurotrophin-3 (NT-3, n=6) were downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/). Then, differentially expressed RNAs (DERs) including DEmiRNAs and DEGs were identified using limma. The noise-robust soft clustering of the intersection DERs was performed using Mfuzz package. Additionally, the integrated miRNAs–targets regulatory network was constructed using Cytoscape. Finally, the Comparative Toxicogenomics Database 2019 update (http://ctd.mdibl.org/) were used to search the central nervous system injury related pathway. Results: A total of 444 DERs including 382 DEGs and 62 DEmiRNAs were screened between group injury and group none whlie 576 DERs including 523 DEGs and 55 DEmiRNAs were screened between group NT-3 and group injury. Moreover, 80 intersections DERs were identified.Two clusters were obtained including cluster 1 (including rno-miR-3072, rno-miR-667-5p and so on) and cluster 2 (OPRL1, GHSR and so on). DREs in cluster 1 were firstly significantly downregulated in group injury and subsequently were significantly upregulated in group NT-3. DERs in cluster 2 were firstly upregulated in group injury and subsequently downregulated in group NT-3. OPRL1 and GHSR were enriched in the KEGG pathway of Neuroactive ligand-receptor interaction which is also found in the Comparative update. OPRL1 was involved in the chemical homeostasis and ion homeostasis while GHSR was related to the regulation of fatty acid metabolic process and regulation of cellular ketone metabolic process. Conclusion: DEmiRNAs rno-miR-3072 and rno-miR-667-5p and DEGs OPRL1 and GHSR might participate in the pathogenesis of neurological injury and the neurotrophin-3 treatment. expressed miRNA, expressed

interaction which is also found in the Comparative Toxicogenomics Database 2019 update. OPRL1 was involved in the chemical homeostasis and ion homeostasis while GHSR was related to the regulation of fatty acid metabolic process and regulation of cellular ketone metabolic process. Conclusion: DEmiRNAs rno-miR-3072 and rno-miR-667-5p and DEGs OPRL1 and GHSR might participate in the pathogenesis of neurological injury and the neurotrophin-3 treatment. Keywords: spinal cord injury, differentially expressed miRNA, differentially expressed gene, regulatory network, Neurotrophin-3 Background Spinal cord injury results in permanent disability due to the limited growth capacity of adult central nervous system axons [1]. Impairment of sensorimotor processing in the spinal cord, reduce mobility and upper limb function due to corticospinal tract lesions in spinal cord injury or stroke [2] is a major cause of disability in the world [3]. In the rat model, it has been found that complete transection of corticospinal pathways in the pyramids could lead to increased spasms, excessive mono-and polysynaptic spinal reflexes and impaired locomotion [4].
Neurotrophin-3 (NT-3) is a growth factor found in many body tissues including the heart, intestines, skin, nervous system and in skeletal muscles including muscle spindles [5]. It has been proven that NT-3 is required for the survival, correct connectivity and function of sensory afferents that innervate muscle spindles [6].
MicroRNAs (miRNAs) with the length from 18-25 nucleotides, is a class of nonprotein coding RNAs which have been shown to be involved in a wide variety of biological processes as regulatory molecules through suppress the mRNA expression or translation. In 2018, Claudia et al. [7] conducted RNA (miRNA and mRNA) sequence (now available at GSE82195 and GSE82196) describing transcriptional changes in cervical dorsal root ganglia after bilateral pyramidotomy and forelimb intramuscular gene therapy with an adeno-associated viral vector encoding human neurotrophin-3 and Many of the dysregulated genes are involved in axon guidance and plasticity. Intramuscular neurotrophin-3 treatment normalized many of those gene changes and may be one of the mechanisms how reflexes, functional recovery and molecular markers in the spinal cord are restored.
In our study, using the same microarray data by Claudia et al., we aimed to further screen the DEmiRNAs and DEGs with linear models for microarray data (limma) package based on a threshold FDR-value < 0.05 and |log2 FC| > 0.5) and clustering analysis using Mfuzz package, miRNAs-mRNAs regulatory network followed by Gene Ontology (GO) functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment as well as the central nervous system injury related KEGG pathway screening. It has been shown that analyses based on differential statistical tests may result in different outcomes [8]. Therefore, we believe some different results may be obtained.

Differentilly expressed RNAs analysis
Limma package in R language [10] (Version 3.34.0, https://bioconductor.org/packages/release/bioc/html/limma.html) to screen the differentilly expressed RNAs (DREs) including differentially expressed genes (DEGs) and differentially expressed miRNA (DEmiRNAs) between group injury and group none, as well as the DEGs and DEmiRNAs between group NT-3 and group injury. We used the FDR-value < 0.05 and |log2 FC| > 0.5 as the cutoff criteria for DEGs and DEmiRNAs. Furthermore, the intersection DERs were also screened. The pheatmap package in R [11] (version 1.0.8, https://cran.r-project.org/package=pheatmap) was employed to conduct the bidirectional hierarchical clustering.

Clustering analysis
The noise-robust soft clustering of the intersection DERs was performed for the samples in three groups using Mfuzz package (version 2.42.0, http://www.bioconductor.org/packages/release/bioc/html/Mfuzz.html) [12]. According to the changes of gene expression at three groups, genes were clustered into different clusters. The default parameters were as follows: minimum standard deviation = 0.3, acore = 0.6.

miRNA-mRNA regulatory networks construction
By using DERs in the clusters, the the miRNA-mRNA regulatory network was constructed. Firstly, DEmiRNA targets was predicted using starBase database Version 2.0, http://starbase.sysu.edu.cn/ . Secondly, the miRNA-mRNA pair with the negative correlation were retained. lastly, the screened miRNA-mRNA pairs 6 were visualized using Cytoscape software with the DEmiRNAs in cluster 1 and their target DEGs in cluster 2, as well as using the DEmiRNAs in cluster 2 and their target DEGs in cluster 1 (Version 3.6.1, http://www.cyto scape.org) [13]. DAVID (version 6.8, https://david.ncifcrf.gov/) was used to perform the GO and KEGG pathway enrichment for the mRNA in the network [14].

Central nervous system injury related KEGG pathway screening
The Comparative Toxicogenomics Database 2019 update (http://ctd.mdibl.org/) were used to search the central nervous system injury related pathway with keyword of central nervous system and injury [15]. Then the searched pathways were further filterd based on the enriched pathways of the genes in the miRNA-mRNA regulatory network. Finally, the central nervous system injury related DERs were screened.

DERs analysis
A total of 444 DERs including 382 DEGs and 62 DEmiRNAs were screened between group injury and group none whlie 576 DERs including 523 DEGs and 55 DEmiRNAs were screened between group NT-3 and group injury. Moreover, 80 intersections DERs were identified. We used volcano plots for the visualization and assessment of the variation (or reproducibility) of RNAs expression (Fig. 1). The bidirectional hierarchical clustering revealed that DERs could distinguish these two groups (Fig.   2).

Clustering analysis
With the aforementioned cut-off criteria, two clusters were obtained including cluster 1 (44 DERs, containing 10 miRNAs and 4 mRNAs) and cluster 2 (36 DERs, containing 10 miRNAs and 26 mRNAs) ( Table 1). DREs in cluster 1 were firstly significantly downregulated in group injury and subsequently were significantly upregulated in group NT-3. DERs in cluster 2 were firstly upregulated in group injury and subsequently downregulated in group NT-3 (Fig. 3).

MiRNA-mRNA network construction
By using the DEmiRNAs in cluster 2 and their target DEGs in cluster 1, the miRNA-mRNA regulatory network was constructed which contained 45 miRNA-mRNA pairs, while by using the DEmiRNAs in cluster 1 and their target DEGs in cluster 2, the miRNA-mRNA regulatory network was constructed which contained 31 miRNA-mRNA pairs (Fig. 4). Functional enrichment showed that 7 GO terms and 3 KEGG patyways were enriched for the DEGs in cluster 1 while 9 GO terms and 1 KEGG pathways were enriched for the DEGs in cluster 2 (Fig. 5)

Central nervous system injury related KEGG pathway screening
In the CTD database, 13 KEGG pathways were searched, among which Neuroactive ligand-receptor interaction (including two genes: OPRL1 and GHSR in cluster 2) was intersected with enriched KEGG pathway for the DEGs in the miRNA-mRNA regulatory network. GHSR could be the target of rno-miR-1839-5p and rno-miR-344-3p while OPRL1could be the target of rno-miR-3072 and rno-miR-667-5p (Fig. 6).

Discussion
In order to further investigate the mechanism of neurological injury and the neurotrophin-3 treatment to improves mobility, we re-analyzed the mRNA and miRNA expression of DRG tissues from the normal adult rat, 10 weeks postpyramidotomy, intramuscular AAV-1 GFP adult rat and 10 weeks post-pyramidotomy, intramuscular AAV-1 prepro-neurotrophin-3 adult rat and screened the two genes OPRL1 and GHSR were firstly significantly upregulated in group injury and 8 subsequently downregulated the with the treatment of NT-3. Besides, rno-miR-1839-5p, rno-miR-3072 and rno-miR-667-5p were firstly downregulated in group injury and subsequently upregulated with the treatment of NT-3. GHSR (growth hormone secretagogue 1 receptor) is the receptor of Ghrelin which is a brain-gut peptide hormone secreted from the stomach to stimulate food intake.
Additionally, GHSR has been found to mediate the neuroprotection in rodents including the activation of UCP2 and decreasing in mitochondrial ROS production, suppression of the pro-inflammatory cytokines TNFa, IL-6 and IL-1b and augmentation of midbrain dopamine neuron electrical activity [16]. The functional enrichment showed that GHSR was related to the regulation of fatty acid metabolic process and regulation of cellular ketone metabolic process. Former study has demonstrated that the ketogenic diet improves forelimb motor function after spinal cord injury in rodents [17]. And Omega-3 polyunsaturated fatty acids could promote functional recovery in rats undergoing spinal cord injury [18]. So, GHSR could be related to neurological injury and the neurotrophin-3 treatment by regulating the cellular ketone metabolic process and fatty acid metabolic process.
OPRL1 encodes a G protein-coupled receptor for nociceptin, an endogenous opioidrelated neuropeptide which is one of four opioid receptors. Opioid receptor activation could attenuate the release of inhibitory neurotransmitters and changes in neuronal excitability [19] which indicated that OPRL1 could played a important role in the nervous system. The GO functional enrichment showed that OPRL1 was involved in the chemical homeostasis and ion homeostasis. It has been illustrated that through a return to homeostasis of chloride after spinal cord injury, exercise could help to contribute to reflex recovery [20]. Additionally, OPRL1 could be the target of rno-miR-3072 and rno-miR-667-5p. MiR-3072 has been found to be related to the decrease of damage and paralysis of lower limbs following spinal cord injury (SCI) [21]. rno-miR-667 may be associated with presence of nerve injury-induced hypersensitivity [22]. Therefore, we speculated that OPRL1 targeted by rno-miR-3072 and rno-miR-667-5p could play an important role in neurological injury and the neurotrophin-3 treatment.
Besides, OPRL1 and GHSR were enriched in the KEGG pathway of Neuroactive ligand-receptor interaction which is found to associated with central nervous system injury in CTD database.
In conclusion, the identified DEmiRNAs rno-miR-3072 and rno-miR-667-5p and DEGs OPRL1 and GHSR might participate in the pathogenesis of neurological injury and the neurotrophin-3 treatment. However, further research is required to validate the results.
Shanshan Yu and Shuang Qi conceived, designed, priformed the experiments and wrote the paper. Zinan Li analyzed and interpreted the data and contributed methods, materials, analysis tools or data. All authors read and approved the final manuscript.  Figure 1 A, the volcano plot of DERs between group injury and group none; B. the volcano plot of DER Figure 2 The two-way clustering of differentially expressed RNAs (DErs); A for the DERs between grou 15 Figure 3 The gene expression changes in two clusters. The color varying from orange to red represent The KEGG pathways and GO terms enrichment for the DEGs in cluster 1 and (A) cluster 2 (B); Figure 6 Pathway regulatory networks associated with central nervous injury.Triangle and circle repre