Abstract
-
Objective
This study aimed to identify differentially expressed genes (DEGs) by transcriptome analysis to elucidate a potential mechanism by which governor vessel electroacupuncture (GV-EA) promotes neuronal survival, axonal regeneration, and functional recovery after complete transection spinal cord injury (SCI).
-
Methods
Sham, control, or GV-EA group adult female Sprague Dawley rats underwent a complete transection SCI protocol. SCI area RNA-seq investigated the DEGs of coding and noncoding RNAs 7 days post-SCI. Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses were used to classify DEGs functions, to explain a possible molecular mechanism. Immunofluorescence and BBB (Basso, Beattie, and Bresnahan) score were used to verify a GV-EA treatment effect following SCI.
-
Results
GV-EA treatment could regulate the expression of 173 mRNA, 260 lncRNA, and 153 circRNA genes among these DEGs resulted by SCI. GO enrichment analysis showed that the DEGs were most enriched in membrane, actin binding, and regulation of Toll-like receptor signaling pathway. KEGG pathway analysis showed enriched pathways (e.g. , Toll-like receptors, MAPK, Hippo signaling). According to the ceRNA network, miR-144-3p played a regulatory role by interacting with lncRNA and circRNA. GV-EA also promoted the injured spinal cord neuron survival, axonal regeneration, and functional improvement of hind limb locomotion.
-
Conclusion
Results of our RNA-seq suggest that post-SCI GV-EA may regulate characteristic changes in transcriptome gene expression, potential critical genes, and signaling pathways, providing clear directions for further investigation into the mechanism of GV-EA in subacute SCI treatment. Moreover, we found that GV-EA promotes neuronal survival, nerve fiber extension, and motor function recovery in subacute SCI.
-
Keywords: Spinal cord injury, Governor vessel electroacupuncture, RNA sequencing, Bioinformatic analysis, CeRNA
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.
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.
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).
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.
SUPPLEMENTARY MATERIALS
Supplementary Fig. 2.
Correlation analysis of mRNA expression. (A) Principal components analysis (PCA) analysis with 3 principal components (PC1, PC2, and PC3) was performed to demonstrate the source of variance in our data. (B) Correlation coefficient cluster analysis, the data in the color block represents the specific value of the phase relationship between 2 samples, and the color change from light to dark indicates the correlation coefficient change from low to high. N1–N3: 3 samples of sham group; C1–C5: 5 samples of control group; EA1–EA5: 5 samples of governor vessel electroacupuncture group.
Supplementary Fig. 3.
Correlation analysis of lncRNA expression. (A) Principal components analysis (PCA) analysis with 3 principal components (PC1, PC2, and PC3) was performed to demonstrate the source of variance in our data. (B) Correlation coefficient cluster analysis, data in the color block represent the specific value of the phase relationship between 2 samples, and the color change from light to dark indicates the correlation coefficient change from low to high. N1–N3: 3 samples of sham group; C1–C5: 5 samples of control group; EA1–EA5: 5 samples of governor vessel electroacupuncture group.
Supplementary Fig. 4.
Correlation analysis of miRNA expression. (A) Principal components analysis (PCA) analysis with 3 principal components (PC1, PC2, and PC3) was performed to demonstrate the source of variance in our data. (B) Correlation coefficient cluster analysis, the data in the color block represent specific values of the phase relationship between 2 samples, and the color change from light to dark indicates the correlation coefficient change from low to high. N1–N3: 3 samples of sham group; C1–C5: 5 samples of control group; EA1–EA5: 5 samples of governor vessel electroacupuncture group.
Supplementary Fig. 5.
Correlation analysis of cicrRNA expression. (A) Principal components analysis (PCA) analysis with 3 principal components (PC1, PC2, and PC3) was performed to demonstrate the source of variance in our data. (B) Correlation coefficient cluster analysis, the data in the color block represent specific values of the phase relationship between 2 samples, and the color change from light to dark indicates the correlation coefficient change from low to high. N1–N3: 3 samples of sham group; C1–C5: 5 samples of control group; EA1–EA5: 5 samples of governor vessel electroacupuncture group.
NOTES
-
Conflict of Interest
The authors have nothing to disclose.
-
Funding/Support
This study was supported by grants from the National Natural Science Foundation of China (82074528; 81774397) and the National Key R&D Program of China (2017-YFA0104702).
-
Author Contribution
Conceptualization: XX, QD, XZ, BQL, YHM, GL, YSZ, YD; Data curation: XX, QD, XZ, BQL, YHM, GL, YSZ, YD; Formal analysis: XX, QD, XZ, BQL, YHM, GL, YSZ, YD; Funding acquisition: YD; Methodology: XX, QD, YD; Project administration: YD; Writing - original draft: XX, QD; Writing - review & editing: YD.
This study was supported by grants from the National Natural Science Foundation of China (82074528; 81774397) to Y. Ding. The funding bodies played no role in the study design, collection, analysis, and interpretation of data, in the writing of the report, or in the decision to submit the paper for publication.
Fig. 1.Differential mRNA expression, gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) pathway analysis. (A) Heatmap of differentially expressed genes (DEGs) between sham (n=3), control (n=5), and governor vessel electroacupuncture (GV-EA) (n=5) groups (p<0.05) with green and red indicating downregulated and upregulated expression, respectively. (B, C) Venn diagrams show DEGs overlap among experimental groups. Green and red represent downregulated and upregulated expression, respectively. (D) Top 30 GO terms in the enrichment analysis. The abscissa represents numbers of differentially expressed mRNAs and the ordinate represents the function in each GO term. Pink: BP; Green: cellular component; Blue: molecular function. (E) Bubble diagram of top 20 significantly enriched KEGG pathways. The abscissa shows the rich factor and ordinate shows the enrichment to KEGG pathways. The circle size marks the number of genes and the color change from red to blue indicates p-value change from low to high.
Fig. 2.Differential lncRNA expression in the injured spinal cord of sham (n=3), control (n=5), and governor vessel electroacupuncture (GV-EA) (n=5) groups. (A) Heatmap of differentially expressed genes (DEGs) between sham, control, and GV-EA samples (p<0.05) with green and red indicating downregulated and upregulated expression, respectively. (B, C) Venn diagrams show overlaps of DEGs between experimental groups. Green and red represent downregulated and upregulated expression, respectively.
Fig. 3.Differential circRNA expression in the injured spinal cord of sham (n=3), control (n=5), and governor vessel electroacupuncture (GV-EA) (n=5) groups. (A) Heatmap of differentially expressed genes (DEGs) between sham, control, and GV-EA samples (p<0.05) with green and red indicating downregulated and upregulated expression, respectively. (B, C) Venn diagram show overlaps of DEGs among experimental groups. Green and red represent downregulated and upregulated expression, respectively.
Fig. 4.CeRNA network construction. (A) LncRNA-miRNA-mRNA ceRNA network. Red ball, green circle and orange triangle indicate mRNA, miRNA and lncRNA, respectively. (B) CircRNA-miRNA-mRNA ceRNA network. Red ball, green circle, and blue diamond indicate mRNA, miRNA, and circRNA, respectively.
Fig. 5.Governor vessel electroacupuncture (GV-EA) can protect neurons survival, promote axon regeneration and improve function recovery. (A, B) The immunofluorescence staining of NeuN and NF in sham, control, and GV-EA group after spinal cord injury (SCI) 4 weeks, scale bars=1,000 µm; 40 µm in (a1–a3; b1–b3). (C) Bar chart shows the number of NeuN+ and Hoe+ in the area of spinal cord injury. (D) Bar chart shows the mean fluorescence intensity of NF positive expression in the area of spinal cord injury. (E) Basso, Beattie, and Bresnahan (BBB) behavioral scores in the different stages of sham, control, and GV-EA group after SCI 4 weeks. n=5/group; *p<0.05 indicates significant difference when compared with the sham group. #p<0.05 indicates significant difference when compared with the control group.
REFERENCES
- 1. McDonald JW, Sadowsky C. Spinal-cord injury. Lancet 2002;359:417-25.
- 2. Fiani B, Kondilis A, Soula M, et al. Novel methods of necroptosis inhibition for spinal cord injury using translational research to limit secondary injury and enhance endogenous repair and regeneration. Neurospine 2021;18:261-70.
- 3. Anjum A, Yazid MD, Fauzi Daud M, et al. Spinal cord injury: pathophysiology, multimolecular interactions, and underlying recovery mechanisms. Int J Mol Sci 2020;21:7533.
- 4. Stenudd M, Sabelström H, Frisén J. Role of endogenous neural stem cells in spinal cord injury and repair. JAMA Neurol 2015;72:235-7.
- 5. Ko CC, Tu TH, Wu JC, et al. Acidic fibroblast growth factor in spinal cord injury. Neurospine 2019;16:728-38.
- 6. Takami T, Shimokawa N, Parthiban J, et al. Pharmacologic and regenerative cell therapy for spinal cord injury: WFNS Spine Committee Recommendations. Neurospine 2020;17:785-96.
- 7. Xu H, Yang Y, Deng QW, et al. Governor vessel electro-acupuncture promotes the intrinsic growth ability of spinal neurons through activating calcitonin gene-related peptide/α-calcium/calmodulin-dependent protein kinase/neurotrophin-3 pathway after spinal cord injury. J Neurotrauma 2021;38:734-45.
- 8. Yang Y, Xu HY, Deng QW, et al. Electroacupuncture facilitates the integration of a grafted TrkC-modified mesenchymal stem cell-derived neural network into transected spinal cord in rats via increasing neurotrophin-3. CNS Neurosci Ther 2021;27:776-91.
- 9. Zeng YS, Ding Y, Xu HY, et al. Electro-acupuncture and its combination with adult stem cell transplantation for spinal cord injury treatment: a summary of current laboratory findings and a review of literature. CNS Neurosci Ther 2022;28:635-47.
- 10. Li WJ, Li SM, Ding Y, et al. Electro-acupuncture upregulates CGRP expression after rat spinal cord transection. Neurochem Int 2012;61:1397-403.
- 11. Li WJ, Pan SQ, Zeng YS, et al. Identification of acupuncturespecific proteins in the process of electro-acupuncture after spinal cord injury. Neurosci Res 2010;67:307-16.
- 12. Ding Y, Yan Q, Ruan JW, et al. Electro-acupuncture promotes survival, differentiation of the bone marrow mesenchymal stem cells as well as functional recovery in the spinal cordtransected rats. BMC Neurosci 2009;10:35.
- 13. Yan Q, Ruan JW, Ding Y, et al. Electro-acupuncture promotes differentiation of mesenchymal stem cells, regeneration of nerve fibers and partial functional recovery after spinal cord injury. Exp Toxicol Pathol 2011;63:151-6.
- 14. Wang JM, Zeng YS, Liu RY, et al. Recombinant adenovirus vector-mediated functional expression of neurotropin-3 receptor (TrkC) in neural stem cells. Exp Neurol 2007;203:123-7.
- 15. Stark R, Grzelak M, Hadfield J. RNA sequencing: the teenage years. Nat Rev Genet 2019;20:631-56.
- 16. Wu Y, Lu X, Shen B, et al. The therapeutic potential and role of miRNA, lncRNA, and circRNA in osteoarthritis. Curr Gene Ther 2019;19:255-63.
- 17. Cheng Y, Su Y, Wang S, et al. Identification of circRNA-lncRNA-miRNA-mRNA competitive endogenous RNA network as novel prognostic markers for acute myeloid leukemia. Genes (Basel) 2020;11:868.
- 18. Salmena L, Poliseno L, Tay Y, et al. A ceRNA hypothesis: the rosetta stone of a hidden RNA language? Cell 2011;146:353-8.
- 19. Shao R, Li C, Chen Y, et al. LncRNA-GAS5 promotes spinal cord repair and the inhibition of neuronal apoptosis via the transplantation of 3D printed scaffold loaded with induced pluripotent stem cell-derived neural stem cells. Ann Transl Med 2021;9:931.
- 20. Zhang Z, Li X, Chen F, et al. Downregulation of LncRNA Gas5 inhibits apoptosis and inflammation after spinal cord ischemia-reperfusion in rats. Brain Res Bull 2021;168:110-9.
- 21. Zhong Y, Du Y, Yang X, et al. Circular RNAs function as ceRNAs to regulate and control human cancer progression. Mol Cancer 2018;17:79.
- 22. Zhao L, Han T, Li Y, et al. The lncRNA SNHG5/miR-32 axis regulates gastric cancer cell proliferation and migration by targeting KLF4. FASEB J 2017;31:893-903.
- 23. Longhurst JC. Defining meridians: a modern basis of understanding. J Acupunct Meridian Stud 2010;3:67-74.
- 24. Zhang K, Liu Z, Li G, et al. Electro-acupuncture promotes the survival and differentiation of transplanted bone marrow mesenchymal stem cells pre-induced with neurotrophin-3 and retinoic acid in gelatin sponge scaffold after rat spinal cord transection. Stem Cell Rev Rep 2014;10:612-25.
- 25. Frendo ME, da Silva A, Phan KD, et al. The Cofilin/Limk1 pathway controls the growth rate of both developing and regenerating motor axons. J Neurosci 2019;39:9316-27.
- 26. Tedeschi A, Dupraz S, Laskowski CJ, et al. The calcium channel subunit Alpha2delta2 suppresses axon regeneration in the adult CNS. Neuron 2016;92:419-34.
- 27. Gay CM, Zygmunt T, Torres-Vázquez J. Diverse functions for the semaphorin receptor PlexinD1 in development and disease. Dev Biol 2011;349:1-19.
- 28. Clarke TB, Weiser JN. Intracellular sensors of extracellular bacteria. Immunol Rev 2011;243:9-25.
- 29. Zhou Y, Wen LL, Li YF, et al. Exosomes derived from bone marrow mesenchymal stem cells protect the injured spinal cord by inhibiting pericyte pyroptosis. Neural Regene Res 2022;17:194-202.
- 30. Huang W, Qu M, Li L, et al. SiRNA in MSC-derived exosomes silences CTGF gene for locomotor recovery in spinal cord injury rats. Stem Cell Res Ther 2021;12:334.
- 31. Li JL, Liang YL, Wang YJ. Knockout of ALOX12 protects against spinal cord injury-mediated nerve injury by inhibition of inflammation and apoptosis. Biochem Biophys Res Commun 2019;516:991-8.
- 32. Liu Z, Yao X, Jiang W, et al. Advanced oxidation protein products induce microglia-mediated neuroinflammation via MAPKs-NF-κB signaling pathway and pyroptosis after secondary spinal cord injury. J Neuroinflammation 2020;17:90.
- 33. Saikia JM, Chavez-Martinez CL, Kim ND, et al. A critical role for DLK and LZK in axonal repair in the mammalian spinal cord. J Neurosci 2022;42:3716-32.
- 34. Islamov RR, Izmailov AA, Sokolov ME, et al. Evaluation of direct and cell-mediated triple-gene therapy in spinal cord injury in rats. Brain Res Bull 2017;132:44-52.
- 35. Rong Z, Huang Y, Cai H, et al. Gut Microbiota disorders promote inflammation and aggravate spinal cord injury through the TLR4/MyD88 signaling pathway. Front Nutr 2021;8:702659.
- 36. Xie C, Shen X, Xu X, et al. Astrocytic YAP promotes the formation of glia scars and neural regeneration after spinal cord injury. J Neurosci 2020;40:2644-62.
- 37. Zhang A, Bai Z, Yi W, et al. Overexpression of miR-338-5p in exosomes derived from mesenchymal stromal cells provides neuroprotective effects by the Cnr1/Rap1/Akt pathway after spinal cord injury in rats. Neurosci Lett 2021;761:136124.
- 38. Dun XP, Parkinson DB. Role of Netrin-1 signaling in nerve regeneration. Int J Mol Sci 2017;18:491.
- 39. Peng Peng, Yu H, Yongjin L, et al. The emerging role of circular RNAs in spinal cord injury. J Orthop Translat 2021;30:1-5.
- 40. Hu YW, Hu YR, Zhao JY, et al. An agomir of miR-144-3p accelerates plaque formation through impairing reverse cholesterol transport and promoting pro-inflammatory cytokine production. PLoS One 2014;9:e94997.
- 41. Wu J, Niu P, Zhao Y, et al. Impact of miR-223-3p and miR-2909 on inflammatory factors IL-6, IL-1ß, and TNF-α, and the TLR4/TLR2/NF-κB/STAT3 signaling pathway induced by lipopolysaccharide in human adipose stem cells. PLoS One 2019;14:e0212063.
- 42. Chen F, Han J, Wang D. Identification of key microRNAs and the underlying molecular mechanism in spinal cord ischemia-reperfusion injury in rats. PeerJ 2021;9:e11454.
, Qingwen Deng1,*
