Electroacupuncture Alleviates Mechanical Allodynia of a Rat Model of CRPS-I and Modulates Gene Expression Profiles in Dorsal Root Ganglia

Complex regional pain syndrome type-I (CRPS-I) is chronic neurological disorder accompanied with devastating pain. Most conventional medical treatments lack effectiveness, making CRPS-I a challenging clinical condition. Electroacupuncture (EA) showed effectiveness in alleviating the pain symptoms of CRPS-I patients. However, the molecular mechanisms underlying EA's therapeutic effect are still not well-understood. Here, we established the rat chronic post-ischemic pain (CPIP) model to mimic CRPS-I and performed repetitive EA on bilateral hind limbs of the CPIP model rats. We then performed RNA-sequencing (RNA-Seq) to study the differences in gene expression, gene networks, and molecular pathways in ipsilateral DRGs innervating the hind limb of the CPIP model rats with and without repetitive EA treatment. Our results found that repetitive EA treatment significantly alleviated mechanical allodynia in bilateral hind limbs of CPIP model rats. RNA-Seq analysis indicated that EA modulated the expression of multiple genes and gene networks in the DRGs of CPIP model rats. Further bioinformatics analysis identified the up-regulation of an array of genes involved in biological process such as neutrophil chemotaxis and immune response in the DRGs of CPIP model rats after EA treatment. Thus, these results suggest that EA may alleviate pain response in CPIP model rats via regulating multiple genes. Our work may help to further advance the understandings of the molecular mechanisms underlying EA's therapeutic effects on CRPS-I and help to identify novel targets for CRPS-I treatment.

Introduction

Complex regional pain syndrome (CRPS) is a chronic neurological disorder with clinical characters of spontaneous pain, hyperalgesia, limb edema, vasomotor instability, and impairment of motor function (1–3). It is subdivided into two categories: CRPS-I (without identifiable nerve injury) and CRPS-II (with identifiable nerve injury) (4, 5). CRPS-I is more prevalent than CRPS-II. Bone fractures, limb trauma, surgery, or ischemia are common triggers for CRPS-I (6, 7). Of note, pain is among the clinical features of CRPS-I that mostly affect the patients, which severely impairs life qualities. It is estimated that the prevalence of CRPS-I is around 30–40% among patients with stroke or bone fracture (8, 9). Unfortunately, there is still no effective medical treatment available for CRPS-I at present (10). Clinical trials have also failed to support the efficacy of many commonly used interventions, making it a challenging neurological disorder in clinic (7, 10). Therefore, alternative therapeutic options for CRPS-I pain are urgently needed.

In order to study the pain mechanisms of CRPS-I, the rat chronic post-ischemic pain (CPIP) model was developed by applying prolonged ischemia and reperfusion to the hind limb (11). The CPIP rat model developed several important characters that recapitulate the clinic features of CRPS-I, including early swelling and hyperemia in the hind limb, followed with chronic mechanical, chemical, and thermal pain hypersensitivities (11, 12). The CPIP rat model has been widely utilized for mechanistic studies of CRPS-I. Recently, we and other peers have identified several important mechanisms that are involved in CPIP-related pain, including peripheral TRPV1/TRPA1 channels, NMDA receptor, reactive oxygen species, spinal CSF1, CXCL12, and NLRP3 inflammasome, etc. (13–18).

One potential therapeutic option for CRPS pain is electroacupuncture (EA), which incorporates traditional manual acupuncture with modern electrotherapy. EA has shown efficacy for treating a variety of pain conditions (19). Recent study demonstrates that EA can significantly reduce pain symptoms and improve life quality of CRPS-I patient, suggesting EA may serve as an alternative method for relieving CRPS-I-related pain (20). Furthermore, our recent work screened the optimal EA frequency and identified 100 Hz as an effective and reliable method for alleviating the pain response of CPIP model rats (15). Nevertheless, the molecular mechanisms underlying EA's therapeutic effects on CRPS pain is still largely unknown.

Dorsal root ganglia (DRGs) contain sensory neurons that receive and convey pain signals from the peripheral to the spinal cord dorsal horn (21). Thus, DRGs have important roles in primary pain signal generation, integration, transduction, and peripheral sensitization (21, 22). In this study, we performed RNA-Seq to study differences in gene expression, gene networks, and molecular pathways in ipsilateral L4-L6 DRGs innervating the hind limb of CPIP model rats with and without 100 Hz EA treatment. We found that EA can modulate the expression of multiple genes in DRGs of CPIP model rats. Interestingly, our data identified the up-regulation of an array of genes involved in neutrophil chemotaxis and immune response in CPIP model rats after EA treatment. Thus, our work may help to further advance the understandings of the molecular mechanisms underlying EA's therapeutic effect on CRPS-I.

Materials and Methods

Animals

Male Sprague-Dawley rats (8–10 weeks, 300–320 g) were purchased from Shanghai Laboratory Animal Center, Chinese Academy of Sciences. Rats were randomly allocated, and five rats were housed per cage. The rats were housed in the Laboratory Animal Center of Zhejiang Chinese Medical University accredited by the Association for Assessment and Accreditation of Laboratory Animal Care (AAALAC) under standard environmental conditions (12 h light–dark cycles and 24 ± 2°C). Food and water were provided ad libitum. The rats were given minimum of 1 week to adapt to the new environment before experiment. All animal care and experimental studies were approved by the Laboratory Animal Management and Welfare Ethical Review Committee of Zhejiang Chinese Medical University (Permission Number: ZSLL-2017-183).

CPIP Rat Model Establishment

CPIP was established through prolonged hind paw ischemia and reperfusion as described previously (11, 12). Anesthesia was induced in all rats with an intraperitoneal injection of 50 mg/kg of sodium phenobarbital and was maintained with an infusion of sodium phenobarbital at 20 mg/kg/h. A Nitrile 70 Durometer O-ring with a 7/32″ (5.5 mm) internal diameter was slide into the larger side of a 1.5 mL Eppendorf tube (with the snap-cap cut off before use). The O-ring was cut off 3 h later for reperfusion. Sham rats received the same anesthetic procedure but the ankle was surrounded with a cut O-ring which did not block blood flow.

EA Treatment

The procedure of EA was carried out according to our previous study (23). Briefly, rats were loosely immobilized, and acupuncture needles of 0.25 mm in diameter were inserted at a depth of 5 mm into bilateral Zusanli (ST36, 5 mm lateral to the anterior tubercule of the tibia) and Kunlun (BL60, at the ankle joint level and between the tip of the external malleolus and calcaneus) acupoints. The needles were connected with HANS acupuncture point nerve stimulator (HANS-200A Huawei Co., Ltd., Beijing, China). The parameters were set as follows: 100 Hz, square wave current output (pulse width: 0.2 ms). Intensities ranging from 0.5 to 1.5 mA (increased by 0.5 mA every 10 min, for a total of 30 min) were delivered for a period of 30 min. The treatment was conducted once daily for 10 consecutive days.

Determination of Mechanical Allodynia

Before baseline test, rats were habituated to the test environment daily for a consecutive 3 days. Rats were individually placed in transparent Plexiglas chambers on an elevated mesh floor and were habituated for 30 min before the test. The mechanical allodynia was determined using a series of von Frey filaments (UGO Basile, Italy) applied perpendicularly to the mid-plantar surface of the hind paws, with sufficient force to bend the filament slightly for 3–5 s according to methods previously used (14, 24). An abrupt withdrawal of the paw and licking and vigorously shaking in response to stimulation were considered pain-like responses. The threshold was determined using the up-down testing paradigm, and the 50% paw withdrawal threshold (PWT) was calculated by the non-parametric Dixon test (25). The behavior tests are conducted by an experimenter blinded to experimental conditions.

Determination of Hind Paw Swelling

Swelling was observed as an increase in hind paw diameter, as measured by a digital caliper, and was calculated as the difference between the basal value and the test value as in our previous study (13). Each rat was measured 3 times, and the mean value was calculated.

Tissue Collection and RNA Extraction

At day 10, rats were deeply anesthetized with sodium pentobarbital (40 mg/kg) and were perfused transcardially with 200 mL 0.9% saline (4°C). After the perfusion, the right (ipsilateral) side L4-6 DRG segments were collected and were immediately preserved in the RNA later solution (Invitrogen, Carlsbad, USA). Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, USA) following the manufacturer's instructions and treated with DNase I to degrade contaminating DNA. The purity and concentration of the samples was assessed by the Nanodrop Spectrophotometer (NanoDrop Products, CA, USA), and the RNA integrity was assessed by the Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA).

RNA-Seq Library Establishment and RNA-Seq

Total mRNAs from three rats of each group were isolated and used to construct sequencing libraries. DNase I was used to degrade double-stranded and single-stranded DNA contaminant in RNA samples. mRNA molecules were purified from total RNA using oligo(dT)-attached magnetic beads. mRNA were fragmented into small pieces using fragmentation reagent. First-strand cDNA was generated using random hexamer-primed reverse transcription, followed by a second-strand cDNA synthesis. The synthesized cDNA was subjected to end-repair and then was 3′ adenylated. Adaptors were ligated to the ends of these 3′ adenylated cDNA fragments. This process is to amplify the cDNA fragments with adaptors from previous step. PCR products are purified with the SPRI beads, and dissolved in EB solution. The double stranded PCR products were heat denatured and circularized by the splint oligo sequence. The single-strand circle DNA (ssCir DNA) was formatted as the final library. Library was validating on the Agilent Technologies 2100 bioanalyzer. The library was amplified with phi29 to make DNA nanoball (DNB) which have more than 300 copies of one molecular. The DNBs were loaded into the patterned nanoarray and single end 50 base reads were generated in the way of sequenced by synthesis. Finally, the fragments were enriched by PCR amplification to construct a library ready for sequencing using BGISEQ-500 by BGI Group (Shenzhen, China).

Bioinformatics Analysis

Primary sequencing data produced by RNA-Seq (raw reads) were subjected to quality control (QC). The information of total reads and mapping ratio reads were shown in Table 1. Raw reads were filtered into clean reads by internal software SOAPnuke (version 1.5.2), followed by: Remove reads in which unknown bases (N) are more than 10%; Remove reads with adaptors; Remove low quality reads (we define the low quality read as the percentage of base which quality is <15 and >50% in a read). QC of alignment was performed to determine if re-sequencing was needed. If the alignment result passed QC, downstream analysis including gene expression, differentially expressed genes, cluster analysis, Gene Ontology (GO) enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, etc. was proceeded with as we described previously (26).

TABLE 1

Table 1. The information of total reads and mapping ratio for Sham, CPIP, and CPIP+EA groups in RNA-Seq.

Cluster Analysis and Screening of Differentially Expressed Genes

Distances of expressed genes were calculated using the Euclidean method (27). The sum of the squared deviations algorithm was used to calculate distance. The cluster analysis and heat map visualization of gene expression patterns was performed using the “pheatmap” package in the R software of Bioconductor. Differentially expressed mRNAs with statistical significance were identified through Volcano Plot filtering. The threshold required for the results to be considered significant was as follows: Q < 0.01 and fold change ≥1.5 or ≤0.667.

Functional Enrichment Analysis of DEGs

Functional enrichment analysis was performed by functional annotation package “clusterProfiler” in R studio software (RStudio, Boston, MA, USA). GO and KEGG enrichment analysis was also conducted. For each enriched function term, P-value of enriched functions and the P-value by multiple testing corrections were calculated by “clusterProfiler” package in R studio software. The GO functional and KEGG pathway enrichment analysis were performed for DEGs using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tools (http://www.geneontology.org and http://www.genome.jp/kegg).

Real-Time Quantitative PCR (qPCR)

The extracted total RNA from the DRGs was reverse transcribed into cDNA using random hexamers primers (TaKaRa Bio Inc., Shiga, Japan) according to the manufacturer's instruction. The sequences of all primers used were shown in Table 2. qPCR was performed using the Fast Start Universal SYBR Green Master kit (TaKaRa Bio Inc., China) 20 μL reaction system according to the manufacturer's protocol by LightCycler480 System (Roche, America). Each reaction was performed in triplicates and normalized to β-actin gene expression. The CT value of each well was determined using the LightCycler480 System software and the average of the triplicates was calculated. The relative quantification was determined by the ΔΔCT method (28, 29).

TABLE 2

Table 2. Sequences of the primers used for qPCR validation of RNA-Seq data.

Protein-Protein Interaction (PPI) Network Analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) is used to provide information regarding predicted and experimental interactions of proteins and the prediction method of this database is from neighborhood, gene fusion, co-occurrence, co-expression experiments, databases, and text mining. By setting the Combination score >0.7 as the reliability threshold value, the web based STRING database (http://string-db.org/) was used to produce PPI predictions after uploading the union gene list to the search bar (30). Based on the interplayed relationships, a PPI network was established and then visualized using the Cytoscape software (31). The connectivity degree of each protein, namely the number of proteins it connected, was calculated to evaluate its importance in this network.

Statistical Analysis

One- or two-way ANOVA followed by Tukey's post-hoc test was used for comparison among groups ≥3. Student's t-test was used for comparisons between two groups. Data in graphs are expressed as means ± SEM. Comparison is considered significantly different if the p < 0.05.

Results

The Establishment of the CPIP Rat Model and Nocifensive Behavior Evaluation

At beginning, we set to establish the rat CPIP model, which mimics human CRPS-I, according to methods described before (11, 12). The hind paw of the CPIP model rats exhibited cyanosis, a sign of tissue hypoxia, and edema when the O-ring was placed (Figure 1). After reperfusion, the hind paw was filled with blood. The edema persisted 3 days and gradually returned back to normal. Furthermore, the CPIP model rats developed obvious mechanical allodynia in both ipsilateral and contralateral hind paws, which persisted until 10 days of our observation time frame. The above observations are consistent with previous results and indicated successful establishment of the CPIP rat model (12).

FIGURE 1

Figure 1. The rat chronic post-ischemic pain model (CPIP) exhibited long-lasting bilateral mechanical allodynia. (A) Representative photographs of rat hind paw taken during CPIP model establishment, 10 min and 7 day after model establishment. The red arrow indicates the paw with the O-ring (ipsilateral side). (B) 50% paw withdraw threshold (PWT) of ipsilateral hind paw. (C) 50% PWT of contralateral hind paw. (D) Summary of the normalized area under the curve (AUC) as in (B). (E) Summary of the normalized AUC as in (C). (F) Thickness evaluation of the ipsilateral hind paw. (G) Thickness evaluation of the contralateral hind paw. n = 6 rats/group. **p < 0.01 vs. Sham group. Student's t-test or two-way ANOVA followed by Tukey post-hoc test was used for statistical analysis. *p < 0.05.

EA Significantly Ameliorated Mechanical Allodynia of CPIP Model Rats

We proceeded to examine the effects of EA on the mechanical allodynia of CPIP model rats. In our recent study, we screened the optimal EA frequency for CPIP treatment and found that daily 100 Hz EA treatment showed effective and reliable analgesic effect on mechanical allodynia of CPIP model rats (15). Therefore, in the present study, we selected 100 Hz EA for the treatment. The procedure of the EA treatment was indicated in Figures 2A,B. EA significantly ameliorated mechanical allodynia in both ipsilateral and contralateral hind paws of CPIP model rats on each time point for 10 days of our observation time frame (Figures 2C,D). The analysis of the area under the curve (AUC) further indicated an overall anti-allodynic effect on bilateral hind paws (Figures 2E,F). These data confirmed that 100 Hz EA is able to attenuate the mechanical allodynia of CPIP model rats.

FIGURE 2

Figure 2. Evaluation of the anti-allodynic effects of electroacupuncture (EA) on CPIP model rats. (A) Experimental protocol for EA treatment. (B) Schematic picture of the locations of acupoints ST36 and BL60 on rat. (C) The effect of bilateral 100 Hz EA on 50% PWT of ipsilateral hind paws of CPIP rats. (D) The effect of bilateral 100 Hz EA on 50% PWT of contralateral hind paws of CPIP rats. (E) Summary of the normalized AUC as in (C). (F) Summary of the normalized AUC as in (D). n = 6 rats/group. **p < 0.01 vs. Sham group. ^#p < 0.05 and ^##p < 0.01 vs. CPIP group. One-way ANOVA or two-way ANOVA followed by Tukey post-hoc test was used for statistical analysis.

Exploring the Differentially Expressed Genes in Ipsilateral Dorsal Root Ganglions (DRGs) After EA Treatment

To determine the effects of EA on gene expression in DRGs of the ipsilateral side of CPIP model rats, we performed RNA-Seq to explore gene expression profiles of ipsilateral L4–6 DRGs 10 days after CPIP model establishment. RNA-Seq generated approximately 21.9 million raw reads in each sample. The ratio of clean reads reached above 99.0% (Table 1). Over 90% of bases had a quality score ≥Q30 and >95% of the clean reads were mapped to rat genome. In all, 20,413 genes were successfully mapped and identified from RNA-Seq. The differentially expressed genes (DEGs) of CPIP group vs. sham group and CPIP+EA group vs. CPIP group are displayed in the volcano plot as in Figures 3A,B. In all, 768 DEGs, including 506 up- and 262 down-regulated, were identified in CPIP vs. sham group, whereas 552 DEGs, including 294 up- and 258 down-regulated, were identified in CPIP+EA vs. CPIP group (Figures 3A,B). The top 20 up- and down-regulated DEGs identified in CPIP+EA vs. CPIP group were further illustrated in Tables 3, 4.

FIGURE 3

Figure 3. RNA-sequencing (RNA-Seq) reveals transcriptome profiles of gene expression changes in the dorsal root ganglion (DRGs) induced by CPIP and EA treatment. (A) Volcano plot showing gene expression profiles in DRGs of CPIP group vs. Sham group. (B) Volcano plot showing gene expression profiles in DRGs of CPIP+EA group vs. CPIP group. Red and blue spots indicate up and down regulated DEGs, respectively, whereas gray spots indicate non-DEGs. (C) Heat map illustration of the hierarchical clustering of DEGs from Sham, CPIP, and CPIP+EA groups.

TABLE 3

Table 3. The detailed information of the top 20 upregulated DEGs.

TABLE 4

Table 4. The detailed information of the top 20 downregulated DEGs.

These DEGs were further summarized in the heat map and hierarchical clustering analysis was performed. We found that the samples from sham group, CPIP model group, and CPIP+EA group were clustered into separate groups, indicating a clear segregation between but not within these 3 groups (Figure 3C). This data indicated that there are significant gene expression changes in ipsilateral DRGs of CPIP model rats vs. sham rats, whereas EA treatment can modify gene expressions in CPIP model rats.

The Analysis of DEGs in DRGs of CPIP Model Rats Receiving EA Treatment

Venn analysis was performed to study the overlapping of DEGs among groups. Venn diagram showed that among the 262 DEGs that were down-regulated in CPIP group, 85 were reversed, whereas 6 were further down-regulated by EA treatment (Figures 4A,B, Supplementary Table 1). Among the 506 DEGs that were up-regulated in CPIP group, 87 were reversed, whereas 15 were further up-regulated by EA treatment (Figures 4C,D, Supplementary Table 2). Figure 4E further summarized the overall Venn diagram results showing the overlapping DEGs among groups.

FIGURE 4

Figure 4. Overlapping of DEGs with genes involved in Sham, CPIP and CPIP+EA group. (A–E) Venn diagrams show overlaps of DEGs between different experimental groups. Eighty-five genes were decreased in CPIP group but increased in CPIP+EA group. Six genes were decreased in both CPIP and CPIP+EA group. Eighty-seven genes were increased in CPIP group but decreased in CPIP+EA group. Fifteen genes were increased in both CPIP and CPIP+EA group.

In order to further explore the mechanisms underlying EA's therapeutic effects, we carried out GO analysis of the DEGs that are up- or down-regulated in CPIP+EA vs. CPIP group of rats. GO analysis showed that the most significantly enriched biological process of upregulated DEGs was neutrophil chemotaxis, followed by defense response to bacterium and fungus and immune response, etc. (Figure 5A, Supplementary Table 3). The most significantly enriched molecular function of upregulated DEGs included oxygen carrier activity, oxygen binding, and heme binding, etc. (Figure 5B, Supplementary Table 3). The most significantly enriched cellular function of upregulated DEGs included hemoglobin complex, extracellular space and extracellular region, etc. (Figure 5C, Supplementary Table 3). In contrast, the most significantly enriched biological process of down-regulated DEGs included muscle contraction, regulation of muscle contraction and retinoic acid biosynthetic process, etc. (Figure 5D, Supplementary Table 3). The most significantly enriched molecular function of downregulated DEGs included retinal dehydrogenase activity, 3-chloroallyl aldehyde dehydrogenase activity, and N-acetyltransferase activity, etc. (Figure 5E, Supplementary Table 3). The most significantly enriched cellular function of downregulated DEGs included extracellular matrix, troponin complex, and collagen trimer, etc. (Figure 5F, Supplementary Table 3).

FIGURE 5

Figure 5. Gene ontology (Go) pathway analysis of DEGs. (A–C) Bubble plots showing the top 10 significant biological processes, molecular functions and cellular components of upregulated DEGs. (D–F) Bubble plots showing the top 10 significant biological processes, molecular functions, and cellular components of downregulated DEGs. Larger bubbles indicate higher number of genes. The color of each bubble reflects the significance (p-value).

Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed to further analyze the DEGs. KEGG analysis indicated that the up-regulated DEGs in CPIP+EA vs. CPIP group were exclusively involved in phagosome, antigen processing and presentation and cytokine-cytokine receptor interaction, etc. (Figure 6A). In addition, the down-regulated DEGs were exclusively involved in phagosome, antigen processing and presentation, cell adhesion molecules (CAMs), etc. (Figure 6B).

FIGURE 6

Figure 6. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs. (A) Bubble plots showing the top 10 significant pathways of upregulated DEGs. (B) Bubble plots showing the top 10 significant pathways of downregulated DEGs. Larger bubbles indicate higher number of genes. The color of each bubble reflects the significance (p-value).

RNA-Seq Data Validation via qPCR

Next, we set to examine the reliability of our RNA-Seq data using qPCR. We evaluated the up-regulated DEGs which were implicated in immune response identified by GO analysis. We found that the expression of 6 genes, including Ccl9 [chemokine (C-C motif) ligand 9], Pf4 (platelet factor 4), Il1rn (interleukin 1 receptor antagonist), Ccl3 (C-C motif chemokine ligand 3), RT1-A1 (RT1 class Ia, locus A1), and Pglyrp1 (peptidoglycan recognition protein 1), were significantly up-regulated by EA treatment (Figure 7A). The trend was consistent with the results derived from RNA-Seq as shown in Figure 7B. Thus, the results of qPCR provide evidence showing that the RNA-Seq data for gene expression profiling was reliable.

FIGURE 7

Figure 7. The validation of RNA-Seq results using qPCR. (A) The expression of representative genes in immune response detected by qPCR. (B) The expression of representative genes by RNA-Seq. n = 7 rat/group. *p < 0.05 and **p < 0.01 vs. with CPIP group. Student's t-test was used for statistical analysis.

PPI Network Analysis of the DEGs Modified by EA Treatment

We further performed protein-protein interaction (PPI) network analysis of the DEGs modified by EA treatment (including 294 up-regulated and 258 down-regulated DEGs as shown in Figures 4A,B). The major hub genes derived from the PPI network analysis included: RatNP-3b (defensin RatNP-3 precursor), Defa11 (defensin, alpha, 11), Defa5 (defensin alpha 5), Elane (elastase, neutrophil expressed), Ppbp (pro-platelet basic protein), Np4 (defensin NP-4 precursor), Cxcr2 (C-X-C motif chemokine receptor 2), and Mmp8 (matrix metallopeptidase 8), etc. (Figure 8).

FIGURE 8

Figure 8. Protein-protein interaction (PPI) network analysis of the DEGs identified from CPIP+EA group vs. CPIP group. PPI network analysis of the DEGs increased or decreased in CPIP+EA group vs. CPIP group. Larger and red color circle reflects more interactions, whereas smaller and purple color circle means fewer interactions.

Discussion

In this work, we studied the anti-allodynic effects of EA on pain-related behaviors of the rat CPIP model. We performed RNA-Seq to study the gene expression profiles of ipsilateral L4-L6 DRGs innervating the hind limbs of CPIP model rats with and without EA treatment. We successfully identified a number of DEGs that are up- or down-regulated by EA treatment. We then studied the cellular and molecular functions of these identified DEGs via GO and KEGG enrichment analysis. Our results indicated that these DEGs were mostly involved in neutrophil chemotaxis, defense response to bacterium/fungus, immune response, and phagosome. We further validated the expression of some representative genes involved in neutrophil chemotaxis and immune response via q-PCR. Finally, PPI network analysis identified a bunch of hub genes that were modulated by EA treatment.

Mechanical allodynia is a prominent clinical feature observed in CRPS-I patients (32). The CPIP model rats developed obvious signs of mechanical allodynia, which well-recapitulates clinical observations in CRPS-I patients. Here, we observed that 100 Hz EA treatment applied right after model establishment significantly alleviated the mechanical allodynia of CPIP model rats. This result is consistent with our recent findings (15). Besides, our recent study also found that 100 Hz EA applied during the maintenance phase of CPIP model can also produce anti-allodynic effects (15). Therefore, these studies all suggest 100 Hz EA can serve as an effective and reliable method for treating CPIP model-related pain behaviors.

The mechanisms through which 100 Hz EA exerts remarkable relieving effect on CPIP-induced mechanical allodynia remain to be examined. In order to explore the genes and molecular pathways which may possibly participate in CPIP model-induced pain and those related to EA-induced pain relieving effect, we carried out RNA-Seq analysis of the ipsilateral DRGs that innervate the hind paw of EA's site of action. Our recent work identified that multiple genes showed expression changes in ipsilateral DRG neurons that innervate the hind limbs of CPIP model rats (26). Our present results are consistent with previous study, showing that CPIP model induced multiple gene expression changes in ipsilateral DRGs (26). We further found that EA treatment partially reversed the up- or down-regulated DEGs in CPIP model rats.

We then focused on the genes that are affected by EA treatment since these genes may be EA-responsive genes and implicated in EA's therapeutic effects. By means of GO analysis, we found that the top enriched biological process of up-regulated DEGs by EA treatment is neutrophil chemotaxis. Neutrophils are usually the first cell type to infiltrate into the tissues when there is noxious stimuli insult (33). Once infiltrated into tissues, they can initiate inflammation, oxidative burst, and clearance of pathogens via recruiting more macrophages (33–35). However, evidence also suggests that infiltrated neutrophils are also capable of driving resolution phase of inflammation via initiation of pro-resolving mechanisms, such as the release of pro-resolving lipid mediators and cytokine scavengers (36–38). These mechanisms can buffer the actions of pro-inflammatory cytokines and exert pro-resolving effects, which enable the conclusion of inflammatory responses (39). Therefore, the exact role of neutrophils in mediating the therapeutic effects of EA on CPIP model rats necessitates further study.

Another mostly enriched biological process exerted by EA we deduced from GO analysis is immune response. Our previous study found that immune response was the mostly enriched biological process in DRGs of CPIP model rats (26). Interestingly, repetitive treatment with EA provoked a further increase in the expression of immune response-related genes in DRGs of CPIP model rats (e.g., Ccl24, Ppbp, Prg2, Ccl3, Ccl9, Cxcr2, Il1rn, and Ccr3, etc.). This finding is consistent with two recent studies exploring the effects of spinal cord stimulation (SCS) on gene expressions in spinal cord of chronic pain model rats (40, 41). SCS is a therapeutic method for attenuating chronic pain (42). Repetitive SCS produced robust analgesic effect on chronic constriction injury (CCI) and paclitaxel-induced chronic pain model rats in these two studies. Further exploration of gene expression changes by RNA-Seq indicated repetitive SCS evoke an obvious increase of expression of many immune response–related genes in the spinal cord of these pain model rats (40, 41). It is usually considered that immune responses in DRGs participate in the initiation and maintenance phase of chronic pain (43). However, more and more evidence suggests that the immune responses can also serve to protect the injured area from further insult by engulfing pathogens, removing cell debris, exerting repair, and neuroprotection mechanisms (44–46). Therefore, EA and SCS may engage similar immune stimulatory effects to exert therapeutic effects on chronic pain conditions. The exact implications of up-regulated expression of immune response-related genes in DRGs after EA treatment of CPIP model rats needs further investigation.

CPIP model rats developed obvious mechanical allodynia in both ipsilateral and contralateral hind limbs, a phenomenon consistent with human (32) CRPS-I patients that exhibit bilateral hypersensitivity to mechanical stimuli (32). This phenomenon, in which trauma or inflammation on ipsilateral side elicits pain on contralateral non-injured side is called mirror-image pain (MIP). At the present, the detailed mechanisms underlying MIP remain largely unknown. It is believed that damage to peripheral or spinal tissues may produce a variety of inflammatory mediators and reactive oxygen species products and initiates inflammatory responses (47, 48). These products may be transported through body or cerebrospinal fluids to the contralateral site, affecting spinal cord, DRGs or peripheral nerves, and producing MIP (49). In addition, glial cell sensitization and their cross-talk in spinal cord are proposed to be involved in MIP (32, 47, 50). Our recent study, together with others, observed significant glial cell activation in bilateral sides of spinal cord dorsal horn of CPIP model rats (13, 18), suggesting glial cell activation may be involved in MIP of CPIP model rats. In the present study, our major focus is to explore the gene changes profile in ipsilateral DRGs after EA treatment. Therefore, we did not examine gene changes in contralateral DRGs of CPIP model rats and examined EA's effects. But these studies will be helpful for further elucidating the mystery underlying MIP of CRPS-I and unraveling EA's therapeutic mechanisms.

In summary, we found that repetitive EA treatment significantly attenuated the mechanical allodynia of the rat CPIP model, an animal model mimicking human CRPS-I. We further identified the genes and gene networks that were affected by EA treatment in DRGs of CPIP model rats. These results suggest that EA may alleviate pain response in CPIP model rats via regulating multiple genes in DRGs. Thus, our work may help to further advance the understandings of the molecular mechanisms underlying EA's therapeutic effects on CRPS-I-related pain.

Data Availability Statement

The original contributions presented in the study are publicly available. The RNA-Seq dataset has been deposited into the National Center for Biotechnology Information's Gene Expression Omnibus repository with accession number GSE158560. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE158560.

Ethics Statement

This animal study was reviewed and approved by the Animal Ethics Committee of Zhejiang Chinese Medical University (permission number ZSLL-2017-183).

Author Contributions

JW, XZ, RC, XL, YL, HN, and DZ performed the experiments and analyzed the data. BoyuL and CY performed bioinformatics analysis of RNA-Seq dataset. XH and YJ coordinated with EA treatment. JF and BoyiL designed and supervised the study. BoyiL wrote the manuscript. All authors reviewed the manuscript.

Funding

This project was supported by Zhejiang Provincial Natural Science Funds for Distinguished Young Scholars (LR17H270001), National Natural Science Foundation of China (81873365 and 81603676), and research funds from Zhejiang Chinese Medical University (Q2019J01 and KC201943).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fneur.2020.580997/full#supplementary-material

References

1. Ott S, Maihofner C. Signs and symptoms in 1,043 patients with complex regional pain syndrome. J Pain. (2018) 19:599–611. doi: 10.1016/j.jpain.2018.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Eldufani J, Elahmer N, Blaise G. A medical mystery of complex regional pain syndrome. Heliyon. (2020) 6:e03329. doi: 10.1016/j.heliyon.2020.e03329

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Freeman R, Edwards R, Baron R, Bruehl S, Cruccu G, Dworkin RH, et al. AAPT diagnostic criteria for peripheral neuropathic pain: focal and segmental disorders. J Pain. (2019) 20:369–93. doi: 10.1016/j.jpain.2018.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Levy RM, Mekhail N, Kramer J, Poree L, Amirdelfan K, Grigsby E, et al. Therapy habituation at 12 months: spinal cord stimulation versus dorsal root ganglion stimulation for complex regional pain syndrome type I and II. J Pain. (2019) 21:399–408. doi: 10.1016/j.jpain.2019.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Urits I, Shen AH, Jones MR, Viswanath O, Kaye AD. Complex regional pain syndrome, current concepts and treatment options. Curr Pain Headache Rep. (2018) 22:10. doi: 10.1007/s11916-018-0667-7

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Birklein F, Ajit SK, Goebel A, Perez R, Sommer C. Complex regional pain syndrome - phenotypic characteristics and potential biomarkers. Nat Rev Neurol. (2018) 14:272–84. doi: 10.1038/nrneurol.2018.20

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Bruehl S. An update on the pathophysiology of complex regional pain syndrome. Anesthesiology. (2010) 113:713–25. doi: 10.1097/ALN.0b013e3181e3db38

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Albazaz R, Wong YT, Homer-Vanniasinkam S. Complex regional pain syndrome: a review. Ann Vasc Surg. (2008) 22:297–306. doi: 10.1016/j.avsg.2007.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Beerthuizen A, Stronks DL, van't Spijker A, Yaksh A, Hanraets BM, Klein J, et al. Demographic and medical parameters in the development of complex regional pain syndrome type 1 (CRPS1): prospective study on 596 patients with a fracture. Pain. (2012) 153:1187–92. doi: 10.1016/j.pain.2012.01.026

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Tran DQ, Duong S, Bertini P, Finlayson RJ. Treatment of complex regional pain syndrome: a review of the evidence. Can J Anaesth. (2010) 57:149–66. doi: 10.1007/s12630-009-9237-0

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Coderre TJ, Xanthos DN, Francis L, Bennett GJ. Chronic post-ischemia pain (CPIP): a novel animal model of complex regional pain syndrome-type I (CRPS-I; reflex sympathetic dystrophy) produced by prolonged hindpaw ischemia and reperfusion in the rat. Pain. (2004) 112:94–105. doi: 10.1016/j.pain.2004.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Hu Q, Zheng X, Chen R, Liu B, Tai Y, Shao X, et al. Chronic post-ischemia pain model for complex regional pain syndrome type-i in rats. J Vis Exp. (2020) 112:94–105. doi: 10.3791/60562

CrossRef Full Text | Google Scholar

13. Hu Q, Wang Q, Wang C, Tai Y, Liu B, Shao X, et al. TRPV1 channel contributes to the behavioral hypersensitivity in a rat model of complex regional pain syndrome type 1. Front Pharmacol. (2019) 10:453. doi: 10.3389/fphar.2019.00453

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Chen R, Yin C, Hu Q, Liu B, Tai Y, Zheng X, et al. Expression profiling of spinal cord dorsal horn in a rat model of complex regional pain syndrome type-I uncovers potential mechanisms mediating pain and neuroinflammation responses. J Neuroinflammation. (2020) 17:162. doi: 10.1186/s12974-020-01834-0

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Hu Q, Zheng X, Li X, Liu B, Yin C, Li Y, et al. electroacupuncture alleviates mechanical allodynia in a rat model of complex regional pain syndrome type-I via suppressing spinal CXCL12/CXCR4 signaling. J Pain. (2020). doi: 10.1016/j.jpain.2020.01.007. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

16. de Logu F, de Pra SD, de David Antoniazzi CT, Kudsi SQ, Ferro PR, Landini L, et al. Macrophages and Schwann cell TRPA1 mediate chronic allodynia in a mouse model of complex regional pain syndrome type I. Brain Behav Immun. (2020 88:535–46. doi: 10.1016/j.bbi.2020.04.037

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Xu X, Tao X, Huang P, Lin F, Liu Q, Xu L, et al. N-methyl-d-aspartate receptor subunit 2B on keratinocyte mediates peripheral and central sensitization in chronic post-ischemic pain in male rats. Brain Behav Immun. (2020) 87:579–90. doi: 10.1016/j.bbi.2020.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Tang Y, Liu L, Xu D, Zhang W, Zhang Y, Zhou J, et al. Interaction between astrocytic colony stimulating factor and its receptor on microglia mediates central sensitization and behavioral hypersensitivity in chronic post ischemic pain model. Brain Behav Immun. (2018) 68:248–60. doi: 10.1016/j.bbi.2017.10.023

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Vickers AJ, Vertosick EA, Lewith G, MacPherson H, Foster NE, Sherman KJ, et al. Acupuncture for chronic pain: update of an individual patient data meta-analysis. J Pain. (2018) 19:455–74. doi: 10.1016/j.jpain.2017.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wei X, He L, Liu J, Ai Y, Liu Y, Yang Y, et al. Electroacupuncture for reflex sympathetic dystrophy after stroke: a meta-analysis. J Stroke Cerebrovasc Dis. (2019) 28:1388–99. doi: 10.1016/j.jstrokecerebrovasdis.2019.02.010

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Berta T, Qadri Y, Tan PH, Ji RR. Targeting dorsal root ganglia and primary sensory neurons for the treatment of chronic pain. Expert Opin Ther Targets. (2017) 21:695–703. doi: 10.1080/14728222.2017.1328057

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Liu B, Linley JE, Du X, Zhang X, Ooi L, Zhang H, et al. The acute nociceptive signals induced by bradykinin in rat sensory neurons are mediated by inhibition of M-type K+ channels and activation of Ca2+-activated Cl- channels. J Clin Invest. (2010) 120:1240–52. doi: 10.1172/JCI41084

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Chai W, Tai Y, Shao X, Liang Y, Zheng GQ, Wang P, et al. Electroacupuncture alleviates pain responses and inflammation in a rat model of acute gout arthritis. Evid Based Complement Alternat Med. (2018) 2018:2598975. doi: 10.1155/2018/2598975

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Li Y, Yin C, Li X, Liu B, Wang J, Zheng X, et al. Electroacupuncture alleviates paclitaxel-induced peripheral neuropathic pain in rats via suppressing TLR4 signaling and TRPV1 upregulation in sensory neurons. Int J Mol Sci. (2019) 20:5917. doi: 10.3390/ijms20235917

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Chaplan SR, Bach FW, Pogrel JW, Chung JM, Yaksh TL. Quantitative assessment of tactile allodynia in the rat paw. J Neurosci Methods. (1994) 53:55–63. doi: 10.1016/0165-0270(94)90144-9

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Yin C, Hu Q, Liu B, Tai Y, Zheng X, Li Y, et al. Transcriptome profiling of dorsal root ganglia in a rat model of complex regional pain syndrome type-I reveals potential mechanisms involved in pain. J Pain Res. (2019) 12:1201–16. doi: 10.2147/JPR.S188758

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Chen W, Ding H, Zhou X, Lin H, Chou KC. iRNA(m6A)-PseDNC: identifying N(6)-methyladenosine sites using pseudo dinucleotide composition. Anal Biochem. (2018) 561–562:59–65. doi: 10.1016/j.ab.2018.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. (2001) 25:402–8. doi: 10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Yin C, Liu B, Wang P, Li X, Li Y, Zheng X, et al. Eucalyptol alleviates inflammation and pain responses in a mouse model of gout arthritis. Br J Pharmacol. (2019) 177:2042–57. doi: 10.1111/bph.14967

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Franceschini A, Szklarczyk D, Frankild S, Kuhn M, Simonovic M, Roth A, et al. STRING v9.1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. (2013) 41:D808–15. doi: 10.1093/nar/gks1094

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Terkelsen AJ, Gierthmuhlen J, Finnerup NB, Hojlund AP, Jensen TS. Bilateral hypersensitivity to capsaicin, thermal, and mechanical stimuli in unilateral complex regional pain syndrome. Anesthesiology. (2014) 120:1225–36. doi: 10.1097/ALN.0000000000000220

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Winterbourn CC, Kettle AJ, Hampton MB. Reactive oxygen species and neutrophil function. Annu Rev Biochem. (2016) 85:765–92. doi: 10.1146/annurev-biochem-060815-014442

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Amulic B, Cazalet C, Hayes GL, Metzler KD, Zychlinsky A. Neutrophil function: from mechanisms to disease. Annu Rev Immunol. (2012) 30:459–89. doi: 10.1146/annurev-immunol-020711-074942

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Yao C, Narumiya S. Prostaglandin-cytokine crosstalk in chronic inflammation. Br J Pharmacol. (2019) 176:337–54. doi: 10.1111/bph.14530

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Jones HR, Robb CT, Perretti M, Rossi AG. The role of neutrophils in inflammation resolution. Semin Immunol. (2016) 28:137–45. doi: 10.1016/j.smim.2016.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Potey PM, Rossi AG, Lucas CD, Dorward DA. Neutrophils in the initiation and resolution of acute pulmonary inflammation: understanding biological function and therapeutic potential. J Pathol. (2019) 247:672–85. doi: 10.1002/path.5221

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Gilroy DW, Bishop-Bailey D. Lipid mediators in immune regulation and resolution. Br J Pharmacol. (2019) 176:1009–23. doi: 10.1111/bph.14587

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Ji RR, Xu ZZ, Strichartz G, Serhan CN. Emerging roles of resolvins in the resolution of inflammation and pain. Trends Neurosci. (2011) 34:599–609. doi: 10.1016/j.tins.2011.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Sivanesan E, Stephens KE, Huang Q, Chen Z, Ford NC, Duan W, et al. Spinal cord stimulation prevents paclitaxel-induced mechanical and cold hypersensitivity and modulates spinal gene expression in rats. Pain Rep. (2019) 4:e785. doi: 10.1097/PR9.0000000000000785

CrossRef Full Text | Google Scholar

41. Stephens KE, Chen Z, Sivanesan E, Raja SN, Linderoth B, Taverna SD, et al. RNA-seq of spinal cord from nerve-injured rats after spinal cord stimulation. Mol Pain. (2018) 14:1744806918817429. doi: 10.1177/1744806918817429

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Mekhail NA, Mathews M, Nageeb F, Guirguis M, Mekhail MN, Cheng J. Retrospective review of 707 cases of spinal cord stimulation: indications and complications. Pain Pract. (2011) 11:148–53. doi: 10.1111/j.1533-2500.2010.00407.x

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Malcangio M. Role of the immune system in neuropathic pain. Scand J Pain. (2019) 20:33–7. doi: 10.1515/sjpain-2019-0138

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Goldstein EZ, Church JS, Hesp ZC, Popovich PG, McTigue DM. A silver lining of neuroinflammation: beneficial effects on myelination. Exp Neurol. (2016) 283:550–9. doi: 10.1016/j.expneurol.2016.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Stephenson J, Nutma E, van der Valk P, Amor S. Inflammation in CNS neurodegenerative diseases. Immunology. (2018) 154:204–19. doi: 10.1111/imm.12922

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Russo MV, McGavern DB. Inflammatory neuroprotection following traumatic brain injury. Science. (2016) 353:783–5. doi: 10.1126/science.aaf6260

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Jancalek R. Signaling mechanisms in mirror image pain pathogenesis. Ann Neurosci. (2011) 18:123–7. doi: 10.5214/ans.0972.7531.11183010

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Xu J, Wei X, Gao F, Zhong X, Guo R, Ji Y, et al. NADPH oxidase 2 derived ROS contributes to LTP of C-fiber evoked field potentials in spinal dorsal horn and persistent mirror-image pain following high frequency stimulus of the sciatic nerve. Pain. (2020) 161:758–72. doi: 10.1097/j.pain.0000000000001761

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Ruohonen S, Jagodi M, Khademi M, Taskinen HS, Ojala P, Olsson T, et al. Contralateral non-operated nerve to transected rat sciatic nerve shows increased expression of IL-1beta, TGF-beta1, TNF-alpha, and IL-10. J Neuroimmunol. (2002) 132:11–7. doi: 10.1016/S0165-5728(02)00281-3

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Choi HS, Roh DH, Yoon SY, Kwon SG, Choi SR, Kang SY, et al. The role of spinal interleukin-1beta and astrocyte connexin 43 in the development of mirror-image pain in an inflammatory pain model. Exp Neurol. (2017) 287:1–13. doi: 10.1016/j.expneurol.2016.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar