INTRODUCTION
Spinal cord injury (SCI) is a destructive trauma that causes a high rate of disability and serious adverse outcomes. Post-SCI pathological changes occur in 2 phases: primary injury and secondary injury, with the latter well established to cause more deleterious impacts [1]. Secondary injury occurs in acute, subacute, and chronic stages, referring to a dynamic, complex cascade of reactions [2]. The acute phase occurs from 2–48 hours post-SCI and leads to axonal injury and cell necrosis. The subacute phase begins on the second day post-SCI and lasts 2 weeks; during this period there is a phagocytic response to clear cellular debris, and early axonal growth occurs. The chronic phase is primarily characterized by glial scar formation and demyelination, beginning 2 weeks post-SCI and persisting long-term [3]. SCI treatment is often impossible during the acute phase, so treatment during the subacute stage is common. Advancing functional recovery requires that we find effective therapies that improve neuronal survival and promote axon regeneration [4-6].
Electroacupuncture (EA) at governor vessel (GV) acupoints (GV-EA), a traditional Chinese treatment in wide use for various SCI problems, has been reported to reduce secondary injury and promote axonal regeneration and functional recovery [7-9]. The GV route is similar to that of the spinal cord in anatomy and function, so that GV-EA produces a marked effect. In our previous research, we demonstrated that GV-EA may play an important role in neuroprotection by regulating annexin A5, collapsin response-mediated protein 2, and calcitonin gene-related peptide (CGRP) expression in SCI [10,11]. In particular, CGRP has been shown to promote regeneration of nerve fibers [12,13]. GV-EA also increases neurotophin-3 (NT-3) level in the injured spinal cord, ameliorating the microenvironment and protecting nerve cells by binding with corresponding receptor TrkC. NT-3 also has anti-inflammatory effects and reduces muscle atrophy [7,8,14]. However, the concrete mechanisms of GV-EA effects in SCI have yet to be clarified.
RNA sequencing (RNA-seq) is a powerful tool used in many fields to describe comprehensive genomic functions. Investigating differential gene expression is among its most commonly used functions [15]. Thus, this study aimed to identify differentially expressed genes (DEGs) by transcriptome analysis of RNA-seq to elucidate a potential mechanism of GV-EA neuroprotection for SCI, and evaluated the interaction networks among long noncoding RNAs (lncRNAs), circular RNAs (circRNAs) and mRNA with microRNAs (miRNAs) by using competitive endogenous RNA (ceRNA) fundamentals [16,17]. The results of bioinformatic analysis showed that GV-EA treatment could regulated 173 mRNAs, 260 lncRNAs, and 153 circRNAs among these differential expression genes after SCI. These DEGs are highly related to inflammation, oxidative stress, apoptosis signaling pathways. In addition, GV-EA also promoted the injured spinal cord neuron survival, axonal regeneration, and functional improvement of hind limb locomotion. To our knowledge, our study is the first to provide a comprehensive description of the mechanism underlying the reparative effects of GV-EA in rats with subacute SCI.
MATERIALS AND METHODS
1. Spinal Cord Injury
The complete transection SCI model used adult female Sprague Dawley rats (aged 2 months, 220–250 g), supplied by the Experimental Animal Center of Sun Yat-sen University, China. All animal experiments were approved by the Institutional Animal Care and Use Committee (IACUC), Sun Yat-Sen University (SYSU) (approved No. SYSU-IACUC-2019-B086). All procedures were compliant with the National Institutes of Health guidelines for the Care and Use of Laboratory Animals. Rats were divided into 3 groups: control (gelatin sponge [GS] graft only); GV-EA (GS graft plus GV-EA); and sham (laminectomy only, without SCI).
The animals were anesthetized with 1% pentobarbital sodium (40 mg/kg, intraperitoneal). For SCI, we exposed the T10 spinal cord segments via laminectomy by cutting the dura mater with sterile microforceps and microscissors. A pair of angled microscissors was used to completely transect the T10 spinal cord. Next, a thin slice of gelatin sponge (GS, 2 mm× 2 mm× 1 mm) was immediately transplanted into the injury gap. Postsurgical care included penicillin (50,000 U/kg/day, intramuscular) for 3 consecutive days, and manual emiction 3 times daily until the end of the experiment.
2. EA treatment
The EA stimulation method has been previously described [4,5,12]. Briefly, we used 2 pairs of GV acupoints: GV9 (Zhiyang)–GV6 (Jizhong) and GV2 (Yaoshu)–GV1 (Changqiang) (Supplementary Fig. 1). The acupuncture needles were manipulated for insertion into the specific acupoints at a depth of 5 mm, and the needles connected with a purpose-made device (model number G6805-2A, Shanghai Medical Electronic Apparatus Company, Shanghai, China), which provided an alternating pattern of sparse and dense wave pulsed stimulation. The stimulation frequencies and durations were applied in a repetition period of 2 Hz for 2.85 seconds and 60 Hz for 1.05 seconds. During EA treatments, the current intensity between acupoint pairs across the lesion site was ~5 μA. EA treatments lasted for 20 minutes every other day, beginning on the third postsurgical day.
3. Tissue Extraction and RNA Sequencing
Seven days postsurgery, rats in all groups were anesthetized with 1% pentobarbital sodium (40 mg/kg, intraperitoneal), then perfused with phosphate buffer saline (PBS). Spinal cords (1-cm segment containing the injury epicenter) were then removed from rats in the sham (n= 3), control (n= 5), and GV-EA (n= 5) groups.
Total RNA of each sample was extracted using TRIzol Reagent/RNeasy Mini Kit (Qiagen, Venlo, The Netherlands). Total RNA of each sample was quantified and qualified by Agilent 2100/2200 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), NanoDrop (Thermo Fisher Scientific Inc., Waltham, MA, USA), and 1% agarose gel. One microgram total RNA was used for following library preparation. Next-generation sequencing library preparations were constructed according to the manufacturer’s protocol. The ribosomal RNA (rRNA) was depleted from total RNA using Ribo-Zero Gold Kit. Then the remaining RNA was fragmented and reverse-transcribed. First strand cDNA was synthesized using ProtoScript II Reverse Transcriptase with random primers and Actinomycin D. The secondstrand cDNA was synthesized using Second Strand Synthesis Enzyme Mix (include dACG-TP/dUTP). The purified doublestranded cDNA by beads was then treated with End Prep Enzyme Mix to repair both ends and add a dA-tailing in one reaction, followed by a T-A ligation to add adaptors to both ends. Size selection of Adaptor-ligated DNA was then performed using beads, and fragments of ~400 bp (with the approximate insert size of 300 bp) were recovered. The dUTP-marked second strand was digested with Uracil-Specific Excision Reagent enzyme. Each sample was then amplified by polymerase chain reaction (PCR) using P5 and P7 primers and the PCR products were validated. Then libraries with different indexes were multiplexed and loaded on an Illumina NovaSeq 6,000 instrument for sequencing using a 2× 150 paired-end configuration according to manufacturer’s instructions. The sequences were processed and analyzed by GENEWIZ, Inc. (Suzhou, China).
4. Tissue Processing Immunofluorescence Staining
All rats were then deeply anesthetized with 1% pentobarbital sodium (50 mg/kg, intraperitoneal) and transcardially perfused with normal saline containing 0.002% NaNO2 and 0.002% heparin, followed by a fixative of 4% paraformaldehyde in 0.1 M PBS (pH, 7.4). The spinal cord was dissected and post-fixed for 24 hours in the same fixative and kept in 30% phosphate-buffered sucrose at 4°C for 48 hours, then frozen and embedded in optimal cutting temperature compound. Successive T8–12 spinal cord segments were cut into longitudinal 25-μm sections. These were then washed with 0.01 M PBS for 5 minutes and blocked with 10% goat serum for 30 minutes at 37°C. The tissue sections were then incubated with primary antibodies diluted with 0.01 M PBS containing 0.3% Triton X-100 at 4°C overnight. After this, they were washed with PBS 3 times for 5 minutes, then incubated with secondary antibodies for 1 hour at 37°C. Cell nuclei were marked with the Hoechst3342 (Hoe). After 3 rinses, the slides were observed with a confocal microscope (Carl Zeiss, Oberkochen, Germany). Primary antibodies included rabbit anti-neurofilament protein 200 (1:200; Sigma-Aldrich, St. Louis, MO, USA) and rabbit polyclonal anti-NeuN (1:500; Abcam, Cambridge, UK). Secondary antibodies included goat anti-rabbit IgG488 (1:1,000; Abcam) and goat anti-rabbit IgG555 (1:1,000; Abcam).
5. Behavioral Assessment
After surgery, hindlimb locomotor function was assessed weekly with the Basso, Beattie, and Bresnahan (BBB) open field locomotor test, on which a score of 0 points indicates severe neurological deficits and 21 points indicates normal performance. Behavioral assessments were performed by 2 experimenters who were blinded to the group information.
6. Computational Analysis for RNA-Seq Data
1) Quality control
To remove technical sequences, including adapters, PCR primers, or fragments thereof, and base quality < 20, pass filter data of fastq format were processed by Cutadapt (v1.9.1) to be high-quality clean data. Then Q20, Q30, and GC contents of cleaning data were calculated. For miRNA, we used the next-generation sequencing data quality statistical software Trimmatic (V0.32) to delete connectors and low-quality sequences in the raw (i.e., pass filter) data for subsequent information analysis. Next, expression data for lncRNA, circRNA, miRNA, and mRNA were analyzed by principal components analysis (PCA) and Pearson correlation coefficient clustering to judge sample repeatability and whether there were remaining abnormal candidates for elimination.
2) Expression analysis and differential expression analysis
We used FPKM method for standardization, which not only normalizes the sequencing depth, but also normalizes the gene length, so that the estimated values of gene expression levels of genes are comparable. The formula is:
Among them, total exon reads/mapped reads (Millions) can be regarded as the percentage of all reads is mapped to this gene, and then divided by the gene length, we can get the percentage of total mapped reads per unit length that is expressed. Then reference genome sequences (version: Rattus_norvegicus. Rnor_6.0.9) and gene model annotation files of relative species were downloaded from genome website, such as UCSC, NCBI, ENSEMBL. Hisat2 (v2.0.1) was used to index reference genome sequence and align clean data to reference genome. In the beginning transcripts in FASTA format are converted from known gff annotation file and indexed properly. Then, with the file as a reference gene file, RSEM (v1.2.15) estimated gene and isoform expression levels from the pair-end clean data. At last, differential expression analysis used the DESeq2 Bioconductor package, a model based on the negative binomial distribution. The estimates of dispersion and logarithmic fold changes incorporate data-driven prior distributions, Padj of genes were set < 0.05 and fold change were set ≥ 1.2 to detect differential expressed ones.
3) GO and KEGG enrichment analysis
Gene ontology (GO) enrichment analysis of DEGs was implemented by the GOseq (v1.34.1), in which gene length bias was corrected. GOSeq was used to identify GO terms on an annotated list of enriched genes with a significant Padj < 0.05. GO enrichment analysis identifies the biological process (BP), cellular component (CC), and molecular function (MF) of DEGs.
Kyoto encyclopedia of genes and genomes (KEGG) is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances (http://en.wikipedia.org/wiki/KEGG). We used in-house scripts to enrich significant DEGs in KEGG pathways.
4) lncRNA/circRNA-miRNA-mRNA interaction analysis
Based on the sample expression data of lncRNA (or circRNA) and mRNA, the Pearson correlation analysis of lncRNA (or circRNA) and mRNA was carried out by using the stats package cor ( ) method of R language, and the results with correlation coefficient greater than 0.95 and p-value of < 0.05 were screened to prepare for the study of regulatory relationship. First, the miRanda software (v3.3a) is used to predict the targeting relationship between miRNA and lncRNA (or circRNA), and then miRanda software (v3.3a) is used to predict the targeting relationship between miRNA and mRNA. Finally, miRNA is used as a connecting bridge to build an indirect regulatory relationship between mRNA and lncRNA (or circRNA) with a correlation coefficient greater than 0.95. Using Cytoscape software (V3.5.1) to draw the regulation network diagram.
7. Statistical Analyses
All data are presented as mean± standard deviation and were analyzed using 1-way analysis of variance (ANOVA) or repeated measures ANOVA with IBM SPSS Statistics ver. 20.0 (IBM Co., Armonk, NY, USA). When variance analysis was satisfied, the least significant difference method was applied for multiple comparisons of the average number of each group. When variance was not uniform, the Kruskal-Wallis test or Dunnett T3 was applied. BBB score data were analyzed by repeated measures ANOVA for homogeneity of variance, then ANOVA was used for between-groups comparisons. A p-value < 0.05 was considered statistically significant.
RESULTS
1. RNA-seq of Spinal Cord Transcriptome
To determine the possible mechanism of GV-EA in the treatment of subacute SCI, we performed whole-genome sequencing analysis to identify possible key regulatory genes following 7-day EA treatment of the injured spinal cord. First, we analyzed fragments per kilo base of exon model per million mapped fragments with DEG-seq software, and filtered out the low-quality spliced reads in the original sequence to eliminate the influence of gene length and sequencing quantity on gene expression, and obtain high-quality data. Q20 and Q30 values of clean data were obtained after processing. There were ~12 G/sample lncRNA for analysis, with Q20 > 97% and Q30 > 93% (Supplementary Table 1). There were ~350 M/sample miRNA data for analysis, with Q20> 99% and Q30> 98% (Supplementary Table 2). There were ~12 G/sample circRNA data for analysis, with Q20> 97% and Q30> 93% (Supplementary Table 3). Thus, the data quality and quantity were reliable.
Through the overall correlation analysis of sample expression data, we found no need to remove any sample. For mRNA, PC score plots showed that the contribution of PC1, PC2, and PC3 were 39.1%, 8.3%, and 6.9%, respectively (Supplementary Fig. 2). The EA5 sample was far from the treatment groups, but the correlation coefficient between the EA5 sample and the other 4 samples in this GV-EA group is higher than 0.95, so this EA5 sample were not excluded. For lncRNA, PC1, PC2, and PC3 were 48.1%, 28.1%, and 10.2%, respectively (Supplementary Fig. 3), and correlation analysis indicated that there were no abnormal samples. For miRNA, PC1, PC2, and PC3 were 15.9%, 10.7%, and 9.5%, respectively (Supplementary Fig. 4); Although the C4 sample in the control group was far from the others, the correlation between them was strong so that this C4 sample was not removed. The PC1, PC2, and PC3 of circRNA were 14.2%, 11.7%, and 10.7%, respectively (Supplementary Fig. 5), and no samples needed to be removed.
2. Identification of DEGs in Injured Spinal Cord following GV-EA Treatment
We analyzed DEGs of mRNA, lncRNA, miRNA, and cicrRNA between the control and GV-EA group according to foldchange ≥ 2 and Padj < 0.05 and found that there was almost no DEGs difference of mRNA, miRNA, and lncRNA. Therefore, we reanalyzed after adjusting the mRNA, miRNA, and lncRNA parameters for foldchange ≥ 1.2 and Padj < 0.05.
To determine the impact of GV-EA on regulation of mRNA expression in the SCI groups, we performed a cluster analysis of the genes that differed significantly. We used EB-seq algorithm to analyze and screen DEGs to identify those that are upregulated and downregulated. We found that 5266 genes were significantly altered among the experimental groups (Fig. 1A). Venn diagrams showed 4636 dramatically DEGs in the control group relative to sham, with 2431 upregulated and 2205 downregulated. As compared with the control group, 704 genes in the GV-EA group displayed differential expression, with 347 upregulated and 357 downregulated genes (Fig. 1B, C). Among the downregulated DEGs after SCI, 59 genes were upregulated after GV-EA treatment (Fig. 1B); 114 of the upregulated DEGs following SCI were downregulated by GV-EA treatment (Fig. 1C).
To provide more precise insight into the DEG functions and relative pathways, we carried out GO analysis and KEGG enrichment analysis of the 3 groups. We selected 30 GO terms with the most significant enrichments (Padj ≤0.05) for display in a GO enrichment histogram: membrane (GO:0016020), cytoplasm (GO:0005737), and integral component of plasma membrane (GO:0005887) were the most enriched in CC. Actin binding (GO:0003779), protein heterodimerization activity (GO:0046982), and protein kinase binding (GO:0019901) were the highest in MF. Toll-like receptor (TLR) signaling pathway (GO:0034123), cytokine secretion (GO:0050707), and cellular response to lipopolysaccharide (GO:0071222) were significantly enriched in BP (Fig. 1D). To further refine the signaling pathways involved in DEG systems, a KEGG bubble chart was used to show the top 20 pathways with the most significantly differentiated enrichment. Mitogen-activated protein kinase (MAPK) signaling pathway, TLR signaling pathway, and cytokine-cytokine receptor interaction were most strongly related to GV-EA therapy for SCI. Twelve DEGs were involved in MAPK signaling pathway, primarily enriching MAP3K13, NCAM1, RasGRP4, MEF2C, STMN1, CACNA2D2, IL1r, EGF, and others. TLR signaling pathway and others contained 68 DEGs, enriching TLR5, TLR2, TLR8, TLR9, and others (Fig. 1E).
LncRNAs, circRNAs, and miRNAs are all noncoding RNAs (ncRNAs), which lack the protein coding ability and participate in mediating various disease responses. LncRNAs are ncRNAs with > 200 nucleotides which participate in a variety of BPs [18]. In addition to analyzing the DEG transcriptome, we evaluated the differential expression of lncRNA. Using the same method described above, we found that compared with the sham group, 1,597 genes were upregulated, and 1,469 genes were downregulated in the control group. Further, 1,287 genes in the GV-EA group relative to the control group show differential expression, with 587 upregulated and 700 downregulated. Moreover, 260 genes coexisted between the 2 comparison groups, allowing further detailed assessment (Fig. 2A–FC).
CircRNA is a newly identified ncRNA, with a stable circular structure. Change in circRNA expression may affect the level of its origin mRNA [19]. DEG enrichment showed that 1,284 genes changed significantly. Compared with the sham group, 380 genes were upregulated and 536 genes were downregulated in the control group. Compared with the control group, 298 genes were upregulated and 323 genes were downregulated in the GV-EA group. Venn analysis showed 153 genes involved in GV-EA repair SCI, including 55 upregulated and 98 downregulated (Fig. 3A–C).
An miRNAs is a short stranded RNA of ~22nt length that can regulate target gene expression [20]. According to the ceRNA theory, competitive binding between sponge RNA and miRNA target genes can regulate miRNA target gene activity [18]. Thus, lncRNAs, or circRNAs can play a considerable role in regulating miRNA levels to change target mRNA expressions [21,22]. To explore the relations between mRNA expression and ncRNA regulation, miRanda, and RNA hybrid algorithms were used to predict the target miRNAs of circRNA, lncRNA, and mRNA, respectively. Then, the lncRNA, circRNA, and mRNA genes with the same expression trends were extracted by sequence clustering analysis, and combined to form an interactive network. The predicted results revealed 4 target genes related to miRNA–mRNA regulation: TLR2, TPM1, LIMK1, and SCARB1. In the lncRNA-miRNA-mRNA interaction network, 43 lncRNAs and 4 mRNAs had the same expression trend (belonging to the increase-decrease type/decrease-increase type). In the circRNA-miRNA-mRNA predicted interaction network, 89 circRNAs and 4 mRNAs had the same expression trend. Further, 9 miRNAs were predicted to act as junctions to interact with lncRNAs and circRNAs to regulate mRNA expressions: rno-miR-338-3p, rno-miR-125a-3p, rno-miR-293-5p, rno-miR-125a-5p, rno-miR-29c-3p, rno-miR-21-5p, rno-miR-144-3p, rno-mir-134-5p, and rno-mir-455-5p (Fig. 4A, B).
3. Effects of Long-term Treatment With GV-EA on Subacute SCI
Based on the RNA-seq results, we further investigated the effect of long-term GV-EA in subacute SCI. Injured rats were treated with GV-EA from post-SCI day 3 for 1 month. We use immunofluorescence to detect neurons (NeuN) and axons (NF200). Four weeks after GV-EA, compared with the control group, the GV-EA group had more surviving NeuN+ neurons (Fig. 5A, C). They also had more NF200+ fibers regenerated into the SCI area (Fig. 5B, D). The hindlimbs of rats were paralyzed after SCI. In the GV-EA group, BBB scores were significantly higher than those of the control group in the period from 2 weeks to 4 weeks post-SCI. At 4 weeks after SCI, the BBB score of GV-EA group reached about 4 points (Fig. 5E). These results indicated that long-term GV-EA may protect neuron survival and promote nerve fiber regeneration in the injured area, facilitating recovery of locomotor function.
DISCUSSION
EA is an ancient, widely used method for treating human disease. It has shown positive effects in SCI, playing a variety of beneficial roles by stimulating special acupoints and conducting them from the meridians and channels to the target site [23]. GV-EA can promote neural survival and differentiation of pre-induced mesenchymal stem cells within GS scaffolds at the SCI lesion site [24]. Further, GV-EA stimulates the internal growth of spinal cord neurons to improve neuronal survival rate, axonal regeneration, and synaptic function by upregulating the CGRP/α-CaMKII/NT-3 pathway [7]. Additionally, GV-EA combined with transplanted stem cells has advantageous anti-inflammatory and neuroprotective effects during early-stage SCI [9]. Morphologically, we found that GV-EA effectively protects NeuN+ neurons and enhances NF+ nerve fiber regeneration over a relatively long period. Those results verified that long-term GV-EA treatment of subacute SCI may play a significant role in promoting neuronal survival, nerve fiber extension, and motor function recovery [7-9]. However, the specific molecular mechanism by which GV-EA treatment impacts subacute SCI has remained unclear. Therefore, we used high-throughput sequencing technology to conduct RNA-seq analysis of injured spinal cord tissues after 7 days of GV-EA treatment for SCI, screening out the related DEGs for biological analysis.
To our knowledge, this is the first study to use RNA-seq to reveal the impacts of GV-EA on SCI; these results are therefore valuable advances and a basis for subsequent studies. We employed reverse transcription to establish cDNA libraries after removing rRNA from total RNA. Thus, further study of the interaction between coding and ncRNA will help better explain these effects. The Q20 and Q30 values, combined with PCA and correlation cluster analysis results demonstrate that our data are of sufficient quality for use in further analyses.
After 7 days of GV-EA treatment, 173 DEGs were regulated, among which 59 were upregulated and 114 were downregulated. The 5 most altered DEGs after SCI were LIMK1, MOBP, PLXND1, CACNA2D2, and Wnt7a. According to a recent study, LIMK1 deletion elevates cofilin activity to promote spinal motor axon growth, as confirmed using a LIMK1−/− mouse model; and following sciatic nerve injury, modulation of the cofilin/LIMK1 pathway accelerates sciatic nerve growth and improves recovery of some sensory and motor functions [25]. Tedeschi et al. [26] demonstrated that CACNA2D2 encodes the voltage-gated calcium channel alpha2delta2 subunit, which limits the development and regeneration of axon growth in central nervous system injury, whereas pharmacological blockade of this subunit by pregabalin administration can enhance axonal regeneration in adult mice following SCI. It has been shown that loss of SEMA3E-PLXND1 signaling changes the synaptic connection; specifically, SEMA3E-induced, PLXND1-mediated rejection may block synaptic formation related to sensorimotor connectivity [27].
To further explain the function and role of DEG mRNAs, we performed GO enrichment and KEGG pathway analysis. GO enrichment analysis suggested that positive regulation of cell death and regulation of cytokine secretion may be related to neuron survival and inflammation after GV-EA treatment in SCI. Consistent with this, the genes for positive cell death regulation include NOD1, CTGF, and ALOX12. NOD1 can activate caspase-1, leading to cell death [28]. Bone marrow mesenchymal-derived exosomes may protect neurons and extend nerve fibers, ultimately accelerating locomotor function recovery by decreasing NOD1 [29]. SiRNA-induced silencing of CTGF can attenuate astrogliosis and scar formation, facilitating recovery of locomotor function following SCI [30]. Transgenic mice with ALOX12 overexpression have promoted inflammation and apoptosis in SCI [31]. Furthermore, the GO results for regulation of cytokine secretion included TLR8, TLR2, TLR9, and TLR5, corresponding to the pathway for enrichment of oxidative stress. Pathway enrichment also showed that MAPK signaling and TLR signaling were notable. MAPK signaling pathway not only participated in inflammatory response, but also mediated spinal cord regeneration and inhibit the secondary injury of SCI [32]. MAP3K13 that involved in MAPK signaling pathway has been shown to regulate injury signals in damaged neurons to regulate regeneration, while in undamaged neurons to regulate sprouting after SCI [33]. NCAM has been demonstrated to promote axonal regeneration as a therapeutic gene [34]. TLR signaling pathway is related to oxidative stress. For example, one function of the TLR4/MyD88 signaling pathway is promoting nerve cell death by inducing inflammation. Gut microbiota dysregulation can also deteriorate SCI by activating the TLR/MyD88 signaling pathway [35]. We also identified some pathways highly related to neuroprotection and axon regeneration. For instance, yes-associated protein (YAP) was upregulated in male C57BL/6 astrocytes and activated in a Hippo pathway-dependent manner post-SCI. Conditional knockout of YAP in astrocytes can significantly inhibit post-SCI astrocyte proliferation in male C57BL/6, as well as reduce glial scar formation, inhibit axon regeneration, and damage behavior recovery [36]. Zhang et al. [37] found that overexpression of miR-338-5p in exosomes derived from mesenchymal stromal cells provides neuroprotective effects via the Cnr1/Rap1/Akt pathway after SCI in rats. It has also been reported that in vivo grafting of Netrin-1 overexpressing fibroblasts into the dorsal column lesion cavity demonstrates that Netrin-1 inhibits axonal regeneration. Further, Netrin-1 was the first axon guidance molecule to be discovered in vertebrates [38]. These cumulative results provide a more precise direction for exploring the mechanism of GV-EA in SCI recovery.
Based on the mRNA analysis results, we focused on ncRNA changes with GV-EA. Previous studies suggested that regulating lncRNA-GAS5 promotes SCI recovery and limits neurocyte apoptosis [19,20]. Our lncRNA results showed that 91 of the overlapping DEGs increased, and 169 decreased, with GV-EA. These differences were greatest in TCONS_00032537, TCONS_00180742, and TCONS_00138251 in upregulated genes, and in TCONS_00125346, TCONS_00236532 and TCONS_00128041 in downregulated genes. Additional recent studies have also shown that circRNAs are highly expressed in neural tissues, regulating neuronal and synaptic functions [39]. Our Venn diagram demonstrates that significantly different circRNAs cooccurred in the 3 study groups, indicating that there were 55 upregulated genes and 98 downregulated genes after GV-EA.
We also analyzed the DEGs of miRNA, predicting involvement of the lncRNA/circRNA-miRNA-mRNA interaction network. We found that 9 key DEG miRNAs (e.g., rno-miR-21-5p, rno-miR-144-3p) are involved in the lncRNA/circRNA-miRNA-mRNA network. These results show that TLR2, TPM1, SCARB1, and other mRNA are associated with multiple noncoding RNAs. From the interaction network, we found that miR-144-3p is associated with TLR2. It has been reported that miR-144-3p mimics (e.g., agomir) can also enhance expression of inflammatory factors, including interleukin (IL)-1β, IL-6, and TNFα [40,41]. Therefore, GV-EA may effectively regulate the interaction between TLR2 and miR-144-3p, further reducing the inflammatory response. After spinal cord ischemia-reperfusion injury, mir-144-3p is upregulated, with the main target gene involved in GMP-PKG and cAMP signaling pathways [42]. Our ceRNA regulatory network results indicate that TCONS_00125346-miR-125a-3p-SCARB1, TCONS_00128041-miR-338-3p-TPM1, and 1:163182618|163216364-miR-134-5p-LIMK1 are markedly altered in the post-SCI pathophysiological process, providing comprehensive insight toward identifying potential GV-EA therapeutic targets. In this case, additional, specific pathways and molecular mechanisms need to be further studied. However, our results suggests that lncRNA and circRNA play a nonnegligible role in the mechanism by which GV-EA reduces the post-SCI inflammatory response.
Additionally, the results of behavioral test and the immunostaining of NeuN and NF200 demonstrated that GV-EA treatment could promote the survival of neurons and axonal regeneration into the lesion site of the injured spinal cord, as well as partial locomotor functional recovery in the spinal cord-transected rats. GV-EA treatment promoted the hindlimbs locomotor functional recovery of the paralyzed rats, which is consistent with results of our previous studies [8,12]. BBB scores gradually increased over time in each group during 8 weeks. BBB scores were significantly higher in the treatment groups compared to the control group in the period from 2 weeks to 8 weeks postinjury, and the BBB score can reach about 9 points at the end of 8 weeks after bone marrow mesenchymal stem cells combined with GV-EA treatment of SCI in rats [8,12]. These results suggested that GV-EA treatment promotes the survival of neurons and stimulates axonal growth into the lesion site, and improves partial locomotor functional recovery in the transected spinal cord in rats.
CONCLUSION
We explored DEGs and the ceRNA network through RNA-seq to establish a clear direction for further study of the mechanism of GV-EA in subacute SCI treatment. The results of RNAseq analysis showed that GV-EA treatment could regulated 173 DEGs of mRNA, 260 DEGs of lncRNA, and 153 DEGs of circRNA among these differential expression genes after SCI. The DEGs changed by GV-EA are highly related to inflammation, oxidative stress, apoptosis, neuroprotection signaling pathways (e.g., TLR, MAPK, Hippo signaling pathway), which may warrant further investigation in the future. We also found that novel ceRNA network through RNA sequencing, which provided a clear direction for further study on the mechanism of GV-EA in the treatment of subacute SCI. Moreover, we found that long-term GV-EA beginning during the subacute SCI stage can promote neuronal survival, nerve fiber extension, and motor function recovery.
PDF Links
PubReader
ePub Link
XML Download
Full text via DOI
Full text via PMC
Download Citation
Supplement
Print
