Gene Transcript Alterations in the Spinal Cord, Anterior Cingulate Cortex, and Amygdala in Mice Following Peripheral Nerve Injury

Introduction

Neuropathic pain characterized by a broad range of sensory, cognitive, and emotional dysfunction is a complex and debilitating public health problem that affects about 7–10% of the gross population worldwide (van Hecke et al., 2014; Bushnell et al., 2015; Colloca et al., 2017). Clinical and preclinical investigations have observed clusters of behavioral symptoms including spontaneous pain, evoked nociceptive behaviors, pain aversiveness, anxiety, and depression in neuropathic pain rodents and patients (Seminowicz et al., 2009; LaCroix-Fralish et al., 2011). These abnormal behaviors may involve the changes in the activities of nociceptive neurons and the emergence of the new pathological processes and signaling pathways (von Hehn et al., 2012; Guo et al., 2016). However, current treatments have not yielded satisfactory results. Opioids and non-steroidal anti-inflammatory drugs (NSAID) are considered effective approaches to relieve these symptoms, but the efficacy should be re-appraised because of possible safety concerns (Finnerup et al., 2015; Jones et al., 2018). Therefore, identifying the gene expression profiles and gene interactions in pain-associated regions is essential for understanding the pathogenesis under neuropathic pain and developing novel therapeutic strategies to improve the treatment outcomes (Samuel and Farsides, 2017).

Abnormal changes in neural activity and plasticity arising from tissue or nerve injury contribute to pain hypersensitivity (von Hehn et al., 2012). The spinal cord is responsible for receiving information from nociceptors and projecting to the brain and plays a major role in integrating and modulating nociceptive signals (Kuner, 2010). Studies have reported that brain regions anterior cingulate cortex (ACC) and amygdala (AMY) are important areas in pain sensation and involved in the interpretation and assessment of the affective and emotional components of pain (Gao et al., 2004; Phelps and LeDoux, 2005; Cao et al., 2009; Navratilova et al., 2015; Neugebauer, 2015). The lesion of ACC and AMY was documented to inhibit the conditioned place aversion of formalin (Gao et al., 2004; Cao et al., 2009). Increasing evidence indicates that cellular and molecular adaptations within these two regions appear under chronic stress and chronic pain conditions (Simons et al., 2014; Ji et al., 2017; Sellmeijer et al., 2018; Navratilova et al., 2019). However, gene expression patterns in these two areas after nerve injury have not been examined. Moreover, previous studies identified a large amount of differentially expressed genes (DEGs) using gene microarrays and RNA sequencing in neuropathic pain (Jiang et al., 2015; Wu et al., 2016; Descalzi et al., 2017; Zhou et al., 2017), but these studies did not provide a full comparison among distinct pain-associated regions as most of them focused on only one region. Thus, it is imperative to have a comprehensive understanding by sequencing and comparing the gene expression patterns in different pain-associated regions.

To this end, in the present study, we carried out a more thorough analysis of gene expression alterations after nerve injury by examining the DEGs in the spinal cord, ACC, and AMY. A mouse model of L4 spinal nerve ligation (SNL) and the next-generation RNA sequencing with a higher sequencing depth were conducted. Our results revealed the unique transcriptional profiles across three regions responding to peripheral nerve injury and the significant overlapping effects implicating in biological functions and signaling pathways despite some differences. Functional analysis demonstrated that the pain-, anxiety-, and depression-related DEGs were closely associated with neuroinflammation and apoptosis. The conspicuous overlap of pain-, anxiety-, and depression-related DEGs among three regions illustrates some conservative changes in the transcriptome independent of regions. Therefore, our findings may bring some useful information and novel insights into the molecular mechanism that will lead to a new direction for further studies and a potential development of clinical analgesic medications.

Methods

Animal Preparation

Eight-week-old male C57BL/6 mice (25–28 g) were purchased from the Animal Experiment Center (Zhengzhou, Henan, China) and housed in the central animal facility under a standard 12-h light/12-h dark cycle with food and water ad libitum. The mice were kept for at least 7 days before the experiments. The procedures for the care and use of animals were approved by the Animal Care and Use Committee of Zhengzhou University and performed under the guidelines of the International Association for the Study of Pain.

L4 SNL-Induced Neuropathic Pain Model

Mice underwent unilateral L4 SNL (the fourth lumbar L4 spinal nerve) as previously described (Wu et al., 2016; Zhou et al., 2017; Sun et al., 2019). Briefly, animals were anesthetized with isoflurane and the L4 spinal nerve was exposed through the removal of the L5 transverse process. After the exposure and isolation of the L4 spinal nerve, a tight ligation with 7–0 silk thread was made and the nerve was transected distal to the ligature. The surgical procedure for the sham group was identical to that of the SNL group, except that the spinal nerves were not transected or ligated.

Behavioral Testing

Animals were habituated to the testing room with a stable temperature for at least 1 day before behavioral measurements.

Von Frey Filament Testing

Mice were put in individual Plexiglas chambers elevated on the mental mesh screen, and 30 min was allowed for adapting to the environment before the testing. The calibrated von Frey filaments (0.07 g and 0.4 g) were used to stimulate the plantar surface of each hind paw for 1 second. Each application represented one trial, and 10 trials were performed for each hind paw. The times were recorded when the animal exhibited a response (withdrawal, flicking, flinching, or licking) to the stimulation for each set of 10 trials. The paw withdrawal frequencies in 10 trials were calculated to evaluate the mechanical sensitivity (Sun et al., 2019; Li Y. et al., 2020).

Hargreaves Assay

Animals were placed in individual Plexiglas chambers on a glass plate and acclimated for 30 min before the testing. The light beam in a Model 336 Analgesia Meter (UGO BASILE S.R.L., Italy) was applied. The radiant heat stimulus was turned off automatically when the animal displayed paw withdrew. The duration between the light application and paw withdrawing was considered as the paw withdrawal latency (PWL). The test for each paw was repeated three to five times at 5-min intervals. A cutoff of 20 s was set to avoid tissue damage (Sun et al., 2019; Li Y. et al., 2020).

Cold Plate Assay

Animals were placed in a chamber on the cold plate with the temperature at 0°C, which was monitored continuously by a thermometer. The duration between the placement and the first sign of mouse jumping and/or flinching was recorded as the paw withdrawal latency (PWL) to noxious cold stimuli. Each test was repeated three times at 10-min intervals for the ipsilateral hind paw. The cutoff of 20 s was set to avoid tissue damage (Sun et al., 2019; Li Y. et al., 2020).

Conditioned Place Aversion (CPA) Testing

The CPA test was carried out as previously described with minor modifications (Guo et al., 2016; Wu et al., 2019). Briefly, the CPA apparatus consisted of two Plexiglas chambers (15 cm × 15 cm × 15 cm) elevated on the mental mesh screen with a removable board in the middle (15 cm × 15 cm). The experimental process included three distinct sessions: a preconditioning session, a conditioning session, and a postconditioning (testing) session with a duration of 10 min in each session. Before the testing, animals experienced the habituation to the apparatus for at least three consecutive days. Six days post-surgery, animals were allowed to freely explore two chambers for 10 min in the preconditioning session. The amount of the time spent in each chamber was recorded. The mice with a strong initial bias (time spent in one chamber > 500 s) were excluded from the study. In the conditioning session, mice were trained in chambers paired with two von Frey filaments. Low-force (0.07 g) von Frey filament was paired with the non-preferred chamber, in which animals spent less time at the preconditioning session. Medium force (0.4 g) was paired with the preferred chamber in which animals stayed longer at the preconditioning session. During the 10-min training, von Frey filaments were used to stimulate ipsilateral hind paw of sham or SNL mice in the corresponding chamber for 5 min with 10-s intervals. Finally, in the postconditioning (testing) session, mice were allowed free access to both chambers. The duration of time that each mouse spent in each chamber was then recorded for 10 min. Results were presented as “Time in chamber” and “CPA score” that was calculated by the time recorded in pre-condition minus the time recorded in the post-condition.

Tissue Collection and RNA Extraction

Briefly, two groups of mice (SNL and Sham) with three biological replicates were used. Unilateral punches were taken from the SC, ACC, and AMY, respectively. The punches per region were pooled from three mice per sample. A total of nine animals per treatment group were needed. Total RNA was extracted using the miRNeasy kit with on-column digestion of genomic DNA (QIAGEN, Valencia, CA, United States) according to the manufacturer’s instructions. RNA was purified with RNeasy Micro Kit 50 (cat. 74004, Qiagen), and the concentration was measured using the NanoDrop 2000 Spectrophotometer (Thermo Scientific, Wilmington, DE). Sample quality was evaluated with the ratios of A260/280 (1.97∼2.08) and RNA integrity numbers (RIN, 7.5∼8.4) as demonstrated by An Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA).

RNA Sequencing

The total RNA (1.0 μg/sample) was subjected to rRNA depletion by Ribo-Zero rRNA Removal (Human/Mouse/Rat) Kit (Illumina, San Diego, CA, United States). Strand-specific RNA libraries were prepared using TruSeq Stranded Total RNA Sample Preparation Kit (Illumina) without poly-A selection. All assays were performed according to the manufacturer’s instructions. RNA-seq was performed on the Illumina Nova6000 plate High Output Model (Illumina, San Diego, United States) (Hrdlickova et al., 2017), in a 2 × 150-bp paired-end configuration, with a total of more than 2,666 M reads per lane (at least 50 M reads per sample).

Bioinformatics Analysis

The samples from the SC, ACC, and AMY were subjected to multiplexing, sequencing, and differential gene expression analysis such as transcript expression analysis and ncRNA expression analysis. Briefly, Trimmomatic 0.32 was used to trim the sequences (minimal length 50 base pairs, leading and trailing Phred Q 30) for the quality first. The resulting sequencing data were then mapped to the musculus genome sequence version GRCm38.72 downloaded from ENSEMBLE. Gene hit counts and reads per kilobase per million (RPKM) were calculated for each gene to determine the expression levels within the CLCbio software environment (CLC Genomics Workbench 7.0.2, CLC genomics Server). Mapped reads were visualized on the UCSC browser using bigwig files converted from bam files. The significant differentially expressed (DE) mRNAs were defined using a cutoff of P < 0.05 and fold change ≥ 1.74 (log2(± 0.8)) to include more DEGs for the subsequent analyses, such as Go term and KEGG pathway analysis, to get more useful information, especially for the functional analysis and comparisons among three distinct regions. The heatmaps were generated via OmicShare (¹ GENE DENOVO). The function of DE mRNAs was analyzed using the downloaded Gene Cards database² and Comparative Toxicogenomics Database³. The DE mRNAs in the SC were mapped to pain-related genes (refer to pain/itch/touch/thermal/chemical related genes), and the DEGs in the ACC and AMY were aligned to pain- and emotion-related genes (refer to anxiety/depression-related genes). Moreover, the DEGs were compared with the genes related to neuroinflammation (inflammation and immunity) and apoptosis. A Venn diagram was also employed for the DEG comparisons⁴.

Quantitative Real-Time RT-PCR

The RNA-sequence results were verified by q-RT-PCR. Total RNA was extracted from tissues as described above, treated using DNase I (New England Biolabs, Ipswich, MA, United States), and finally reversely transcribed with the Revert Aid First Strand cDNA Synthesis Kit (Thermo) according to the manufacturer’s instructions either oligo (dT) primers or specific RT primers. A template (2 μL) was used for the amplification by real-time PCR with primers as shown in Table 1 (Sangon Biotech, Shanghai, China). Each sample was run in triplicate in a 20-μL volume for the reaction containing 250 nM forward and reverse primers, 10 μL Thermo Scientific Maxima SYBR Green qPCR Master Mix (2×; Thermo Scientific Maxima SYBR Green qPCR Master Mix, Rox solution provided), and 20 ng total cDNA. Reactions were implemented in a 7500 Fast Real-Time PCR Detection System (Applied Biosystems, United States). The cycle parameters were set as follows: an initial 3-min incubation at 95°C, followed by 40 cycles of 95°C for 10 s, 60°C for 30 s, and 72°C for 30 s. Ratios of ipsilateral-side mRNA levels to contralateral-side mRNA levels were calculated using the Δ Ct method (2^−ΔΔCt). All data were normalized to Tuba1α (Wu et al., 2016), which was identified as stable in mice after nerve injury.

TABLE 1

Table 1. Primers for RT-qPCR.

Functional Enrichment Analysis of Differentially Expressed Genes (DEGs)

For the function analysis, about 2,230, 1,689, and 1,812 DE mRNAs (P < 0.05, fold change ≥ 1.74), respectively, from the SC, ACC, and AMY were categorized using the Kyoto Encyclopedia of Genes and Genomes (KEGG pathway analysis) and Gene ontology analysis by the database for Annotation, Visualization and Integrated Discovery (DAVID⁵) (Li M. et al., 2020). Likewise, GO Annotations and KEGG Pathways Analysis were applied to predict the role of DE lncRNAs through their target mRNAs. The genetic regulatory networks were clarified by forming hierarchical categories according to the BP, MF, and CC aspects⁶ (Zhou et al., 2017). The significant pathway enrichments of differentially expressed lncRNAs were predicted by the Pathway Analysis⁷ (Zhou et al., 2017).

Protein–Protein Interaction Network Construction

To further analyze and elucidate the functional connection between the differentially expressed encoding genes, the interaction among the significant DEGs in three regions was predicted by the STRING database (version: 11.0)⁸. The top 50 DEGs with the highest correlation degree were screened out to establish the network in the Cytoscape program (version: 3.6.0⁹) (Liu et al., 2019). The connection degree of each node was calculated through the cityscape plugin. The node size was defined by the connection degree. Red and blue colors represented the up- and downregulated genes, respectively.

Co-expression Network Construction

LncRNA–mRNA co-expression networks were established based on the correlation between DE lncRNAs and their neighboring, overlapping, or distant mRNAs in the genome. DEGs with the Pearson correlation coefficients (PCC) > 0.95 or < –0.95 between lncRNAs and mRNAs (FDR < 0.05) were picked out to draw networks by the Cytoscape program (Dong et al., 2014).

Statistical Analysis

All data were collected randomly and expressed as mean ± SEM. The data were statistically analyzed with two-tailed, unpaired Student’s t-test, and one-way and two-way ANOVA with repeated measures. When ANOVA showed a significant difference, pairwise comparisons between means were tested by the post hoc Tukey method. Values of P < 0.05 were considered statistically significant. The data were analyzed by GraphPad Prism 8.0.

Results

SNL Leads to Nociceptive Hypersensitivities and Pain Aversiveness Behaviors

Consistent with early reports (Jiang et al., 2015; Wu et al., 2016), mice exposed to SNL showed significant mechanical allodynia (Figures 1A,B), thermal hyperalgesia (Figure 1C), and cold hyperalgesia (Figure 1D) as indicated by the increases in paw withdrawal frequencies in response to von Frey filament stimuli and the decreases in paw withdrawal latencies in response to heat and cold stimuli, respectively, on day 7 post-surgery on the ipsilateral side as compared to those on the contralateral side. As expected, no nociceptive hypersensitivities were observed on either side of sham-operated mice (Figures 1A–D). Previous studies indicated that peripheral nerve injury led to emotional aversion on day 3 postsurgery (Suzuki et al., 2007; Wu et al., 2019). In line with these studies (Guo et al., 2016; Wu et al., 2019), the time spent in the preferred chamber at initial was sharply reduced in SNL mice on day 7 postsurgery after the repeated stimuli of 0.4 g von Frey filament (Figure 1E; ^∗P < 0.05), demonstrating that SNL mice exhibited emotional aversion. As expected, sham-operated animals did not show any differences when they received the same training (Figure 1E; P > 0.99). The difference in the CPA score between SNL- and sham-operated groups was statistically significant (Figure 1F). Taken together, SNL mice exhibited well-established nociceptive hypersensitivities and emotional aversiveness.

FIGURE 1

Figure 1. Unilateral L4 spinal nerve ligation (SNL) produced nociceptive hypersensitivities and pain aversiveness in mice. (A,B) Paw withdrawal frequencies in response to 0.07 g (A) and 0.4 g (B) von Frey filament stimuli on the ipsilateral and contralateral sides on day 7 following SNL or Sham surgery. n = 18/group. ****P < 0.0001 versus the corresponding sham group or the corresponding contralateral side by two-way ANOVA with repeated measures followed by post hoc Tukey test. (C,D) Paw withdrawal latencies in response to thermal (C) and cold (D) stimuli on the ipsilateral and contralateral sides on day 7 post-SNL or -Sham surgery. n = 18 mice/group for thermal test and 14 mice/group for the cold test. ****P < 0.0001 versus the corresponding sham group or the corresponding contralateral side by one-way (D) or two-way (C) ANOVA with repeated measures followed by post hoc Tukey test. (E,F) Time spent in the corresponding chamber paired with 0.07 g or 0.4 g von Frey filament stimuli (E) and CPA score (F) in the sham-operated and SNL-operated groups. n = 14 mice/group. Pre: Precondition. Post: Postcondition. N.S: not significant, P > 0.99. *P < 0.05 versus the corresponding precondition (E) or sham group (F) by two-way ANOVA with repeated measures followed by post hoc Tukey test (E) or two-tailed unpaired Student’s t-test (F).

RNA-Seq and Genome-Wide Read Mapping in the SC, ACC, and AMY After SNL

More than 50 million (M) reads in each group per region (SC: 55.88 M–95.96 M in sham and 80.78 M–86.10 M in SNL; ACC: 90.74 M-102.62 M in sham and 68.88 M–86.05 M in SNL; AMY: 98.44 M–111.45 M in sham; and 93.65 M–126.80 M in SNL) were achieved. After the trimmed reads were mapped to the reference mouse genome from ENSEMBLE (GRCm38.90), mapped reads were sorted through as exonic, intronic, and intergenic. The proportion of the reads within each category in sham and SNL groups from the SC, ACC, and AMY were illustrated in Supplementary Figure 1A. As expected, many reads were aligned to exonic regions in both groups followed by a considerable proportion of reads mapped to intronic regions (Supplementary Figure 1A). The reads mapped to intergenic regions accounted for a small percentage in both sham and SNL groups among the three regions (Supplementary Figure 1A). Furthermore, Supplementary Figures 1B,D illustrated a robust elevation in the level of reads mapped to the exonic region and a remarkable reduction in the proportion of reads aligning to intronic regions (^∗P < 0.05; ^∗∗P < 0.01) in the SC and AMY, implicating the changes in the functional proteins and signaling pathways after SNL. However, no significant changes were observed in the ACC (Supplementary Figure 1C).

We then analyzed the expression profiles of the DEGs in the SC, ACC, and AMY. Six days after SNL, approximately 38,584, 38,045, and 38,727 genes out of a total of 102,711, 101,218, and 102,589 transcripts, respectively, were identified in the SC, ACC, and AMY. The numbers of the changed genes in three regions were quite similar. In agreement with the changes in the injured DRGs (Wu et al., 2016), the largest transcriptional changes were observed in protein-coding RNAs (49%), followed by other non-coding RNAs (43–44%) and lncRNAs (7–8%) in three regions on day 7 after SNL (Supplementary Figures 1E–G).

Altered Expression Profiles of mRNAs and lncRNAs in the SC, ACC, and AMY After SNL

The robust changes in gene expression of mRNAs and lncRNAs within the spinal cord, ACC, and AMY were observed after nerve injury. About 2,230 (1,616 upregulated, 614 downregulated), 1,689 (1,022 upregulated, 667 downregulated), and 1,812 (1,256 upregulated, 556 downregulated) mRNAs were significantly changed in the SC, ACC, and AMY, respectively (Supplementary Material 1). Besides, approximately 196 (30 upregulated, 136 downregulated), 94 (52 upregulated, 42 downregulated), and 131(50 upregulated, 86 downregulated) lncRNAs were significantly altered in the SC, ACC, and AMY, respectively (Supplementary Material 2). The clustered heatmaps of DE mRNAs (Figure 2A) and DE LncRNAs (Figure 2B) revealed distinct gene expression patterns across three regions after SNL.

FIGURE 2

Figure 2. Differentially expressed mRNAs and lncRNAs in the SC, ACC, and AMY after nerve injury. (A,B) Heatmaps of significantly differentially expressed mRNAs (A) and lncRNAs (B) in the SC, ACC, and AMY from mice on day 7 post-SNL or -sham surgery. n = 9 mice/group. (C,D) Venn diagrams indicate the co-upregulated mRNAs (C) and co-downregulated mRNAs (D) in the SC, ACC, and AMY on day 7 postsurgery. n = 9 mice/group. (E,F) Venn diagrams indicate the co-upregulated lncRNAs (E) and co-downregulated lncRNAs (F) in the SC, ACC, and AMY on day 7 postsurgery. n = 9 mice/group.

Venn diagram was applied to further determine whether these DEGs showed co-expression patterns across three regions. The analyses characterized the co-expression genes by comparing up- and downregulated mRNAs and lncRNAs in the SC, ACC, and AMY. We found the co-regulation of DE mRNAs (Figures 2C,D) and DE lncRNAs (Figures 2E,F) in all three regions, but the robust co-expression patterns were seen in the co-upregulated mRNAs (Figure 2C), and the co-downregulated lncRNAs (Figure 2F), especially in the SC and AMY (Figures 2C,F). The detailed co-expressed DE mRNAs and DE lncRNAs were listed in Supplementary Materials 1, 2, respectively.

Highest Differentially Expressed G Protein-Coupled Receptor mRNAs, Ion Channel mRNAs, and lncRNAs in the SC, ACC, and AMY After SNL

G protein-coupled receptors (GPCRs), ion channels, and lncRNAs are critical in the transmission and modulation of nociceptive information (Zhao et al., 2013; Campbell and Smrcka, 2018; Klein and Oaklander, 2018). Besides 7–8% of the whole transcriptome that are lncRNAs, about 617, 552, and 598 DEGs were identified as GPCR mRNAs and 257, 253, and 250 DEGs were identified as ion channel mRNAs in the SC, ACC, and AMY, respectively. The top 15 up- and downregulated DEGs of GPCR, ion channel, and lncRNA across three regions were displayed in the heatmaps (Figure 3 and Supplementary Figures 2, 3). Consistent with previous reports, the levels of the GPCRs P2ry12 (Horváth et al., 2014), Gpr151 (Jiang et al., 2018), Prokr2 (Maftei et al., 2014), Ccr1 and Ccr5 (Eltayeb et al., 2007; Pevida et al., 2014) in the spinal cord (Figure 3A); D1 receptor (D1R) (Darvish-Ghane et al., 2020) and Cmklr1 (Guo et al., 2012; Doyle et al., 2014) in the ACC (Supplementary Figure 2A); and GPCRs C5ar2 (Carpanini et al., 2019) and Cckbr (Bowers and Ressler, 2015) in the AMY (Supplementary Figure 3A) were remarkably increased on day 7 after SNL. In contrast, the amounts of GPCR histamine receptor H1R (Huang et al., 2017) and G protein-coupled receptor 35 (Gpr35) (Cosi et al., 2011) in the spinal cord (Figure 3A), Grm5 (Ramos-Prats et al., 2019) in the ACC (Supplementary Figure 2A), and Gpr35 (Savitz et al., 2015) and Lpar1 (Pedraza et al., 2014; González de San Román et al., 2019) in the AMY (Supplementary Figure 3A) were dramatically decreased after SNL. For the ion channels, we observed the observably elevated expression of Cacna2d2 (Yu et al., 2019), Orai1 (Dou et al., 2018), Aqp9 (Wu et al., 2020), and Trpc6 (Jin et al., 2017; Wang et al., 2020) in the spinal cord (Figure 3A), Gria1 (Toyoda et al., 2009) and Cacna1c (Jeon et al., 2010) in the ACC (Supplementary Figure 2B), and Cacna1c (Temme and Murphy, 2017), Cacna2d1 (Chen et al., 2018; Young et al., 2016), and Trpc6 (Jin et al., 2017; Wang et al., 2020) in the AMY (Supplementary Figure 3B) after SNL. On the contrary, significant reductions were seen in the levels of ion channel transcripts for Kcnq5 (Manville and Abbott, 2018), Cacna1b (Stevens et al., 2019), and Kcnj6 (Zheng et al., 2015) in the spinal cord (Figure 3B), Gabrb3 (Tripp et al., 2012), Gabra1 (Chandley et al., 2015), and Grik2 (Chandley et al., 2015) in the ACC (Supplementary Figure 2B), as well as the downregulation of the ion channels such as Scn1a, Ano1, Cacna1h (Gangarossa et al., 2014), Cacna1d (McKinney et al., 2009), and Gabra1 (Guilloux et al., 2012) in the AMY (Supplementary Figure 3B) on day 7 after SNL. Interestingly, we detected the upregulation of P2ry12 mRNA in the AMY, which was inconsistent with an earlier report, in which the level of P2ry12 mRNA was unaltered in the AMY post-nerve injury (Barcelon et al., 2019).

FIGURE 3

Figure 3. Heatmaps of the representative differentially expressed genes (DEGs) in the spinal cord after nerve injury. (A-C), Top 15 up- and downregulated genes of G protein-coupled receptors (A), ion channels (B), and lncRNAs (C) in the fourth lumbar spinal cord on day 7 after SNL. Colors in the heatmaps indicate the Row Z-score among the different datasets. The up- and downregulated genes are colored in red and green, respectively. n = 9 mice/group.

Additionally, top 15 up- and downregulated lncRNAs in three regions were shown in heatmaps (Figure 3 and Supplementary Figures 2, 3C). As expected, some DEGs for lncRNAs that were identified in the present study (Table 2) have been previously reported to implicate in pain (Li et al., 2018; Tang et al., 2018; Che et al., 2019; Meng et al., 2019; Han et al., 2020; Wen et al., 2020; Wu et al., 2020). Consistently, we detected the increased expression of lncRNA Dancr and H19, and the decreased expression of lncRNA Meg3 in the spinal cord (Li et al., 2018; Tang et al., 2018; Che et al., 2019; Han et al., 2020; Wen et al., 2020). The downregulation of lncRNA Malat1 in the spinal cord and the upregulation of lncRNA Meg3 and lncRNA H19 in the ACC were observed on day 7 after SNL, but the functions of these lncRNAs under neuropathic pain conditions need to be further determined (Meng et al., 2019; Wu et al., 2020).

TABLE 2

Table 2. Spinal nerve ligation-induced differentially expressed lncRNAs that have been previously reported implicated in pain.

Validation of the DEGs for lncRNAs and mRNAs in the SC, ACC, and AMY

We next conducted a quantitative real-time RT-PCR assay to validate the reliability of RNA sequencing results by analyzing the expression of significant DE lncRNAs and mRNAs on day 7 after SNL in three regions. The expression of three lncRNAs (Pantr1, Mir9-3hg, and Miat), two ion channel mRNAs (Kcnk18 and Kcnma1), and three G protein-coupled receptor mRNAs (P2ry12, Cmklr1, and Adrb3) was measured in the SC (Figure 4A), ACC (Figure 4B), and AMY (Figure 4C), respectively. As expected, the levels of the selected lncRNAs and mRNAs were concomitant with the sequencing results (Figures 4A–C). It was noteworthy that the amount of Kcnma1 was elevated in the AMY (Figure 4C) but reduced in the ACC (Figure 4B).

FIGURE 4

Figure 4. Validations of differentially expressed lncRNAs and mRNAs in the SC, ACC, and AMY after nerve injury. (A) Levels of lncRNA Pantr1, G protein-coupled receptor P2ry12, and ion channel Kcnk18 in the spinal cord on day 7 after SNL. n = 12 mice/group. *P < 0.05; **P < 0.01, and ***P < 0.001 versus the corresponding sham group by two-tailed unpaired Student’s t-test. (B) Amounts of lncRNA Mir9-3hg, G protein-coupled receptor Cmklr1, and ion channel Kcnma1 in the ACC on day 7 after SNL. n = 12 mice/group. *P < 0.05; **P < 0.01, versus the corresponding sham group by two-tailed unpaired Student’s t-test. (C) Levels of lncRNA Miat, G protein-coupled receptor Adrb3, and ion channel Kcnma1 in the AMY on day 7 after SNL. n = 12 mice/group. *P < 0.05; **P < 0.01, and ***P < 0.001 versus the corresponding sham group by two-tailed unpaired Student’s t-test.

Functional Enrichment Analysis of the Differentially Expressed Genes After SNL

To explore the functional enrichments of these DEGs, we performed Gene Ontology and KEGG pathway analyses to categorize the up- and downregulated mRNAs based on the distinct processes using the DAVID bioinformatics database. The top 10 analyzed results of biological processes from the up- (red panels on the left) and downregulated (blue panels on the right) mRNAs in three regions were displayed in Figure 5 and Supplementary Material 3. The most significant enriched biological processes of upregulated genes in the spinal cord were immune system process, apoptotic process, innate immune process, inflammatory response, defense response, regulation of RNA transcription, and cell proliferation, while the downregulated genes in the spinal cord were mainly involved in protein phosphorylation, covalent chromatin modification, negative regulation of NF-kappaB transcription factor activity, and cytokine production after nerve injury (Figure 5A). The upregulated DEGs in the ACC were highly enriched in transcription, regulation of transcription, signal transduction, locomotory behavior, and sensory perception of pain and G protein-coupled receptor pathways, in contrast to the prominent enrichments in transport, cell differentiation, neuron migration, and positive regulation of synapse assembly for downregulated genes (Figure 5B). The upregulated genes in the AMY were related to cell proliferation, apoptotic process, nerve system development, and Histone H4 acetylation besides the enrichments in transcription-related processes, whereas the downregulated genes in the AMY were markedly enriched in regulation of membrane potential, response to wounding, and neuron projection extensive as well as transport processes (Figure 5C). For the molecular function enrichments, we observed the striking enrichments in protein binding, DNA binding, and protein homodimerization activity for the upregulated genes and metal ion binding, ATP binding, and transferase activity for downregulated genes in the spinal cord (Supplementary Material 3). In the ACC, the upregulated genes were prominently implicated in protein binding, action binding, ion channel binding, and protein kinase binding, while the downregulated genes were involved in protein binding, transcription factor binding, and DNA binding (Supplementary Material 3). In the AMY, the upregulated genes were distinctly enriched in protein binding, protein N-terminal binding, and chromatin binding, while the downregulated genes were mainly enriched in lipid binding, endopeptidase inhibitor activity, and protein homodimerization activity (Supplementary Material 3). Within the category of “cellular component,” the DEGs in three regions were robustly enriched in membrane, cytoplasm, and nucleus (Supplementary Material 3).

FIGURE 5

Figure 5. Biological process analysis of the differentially expressed up- and downregulated mRNAs in the SC, ACC, and AMY after nerve injury. (A–C) Analysis of the Gene Ontology database showed top 10 biological processes from upregulated mRNAs (red panels on the left) and downregulated mRNAs (blue panels on the right) in the SC (A), ACC (B), and AMY (C) on day 7 post-SNL according to the P-value. The DAVID database was used to do the GO enrichment analysis. Red and blue bars represent up- and downregulated mRNA enrichments, respectively.

Pathway analyses showed that most significant pathway enrichments in the spinal cord contained chemokine signaling pathway, tumor necrosis factor (TNF) signaling pathway, and Fc gamma R-mediated phagocytosis for the upregulated genes (red panels on the left) and axon guidance, hippo signaling pathway, and T cell receptor signaling pathway for the downregulated genes (blue panels on the right) (Figure 6A). In the ACC, the obvious enrichments were seen in the cAMP signaling pathway, neurotrophin signaling pathway and morphine addiction for the upregulated genes and the insulin signaling pathway, Ras signaling pathway, inflammatory mediator regulation of TRP channels, and type II diabetes mellitus for the downregulated genes (Figure 6B). In the AMY, the dramatic enrichments were detected in the MAPK signaling pathway, osteoclast differentiation, TNF signaling pathway for the upregulated genes and the metabolic pathways, chemical carcinogenesis, and complement and coagulation cascades for the downregulated genes (Figure 6C). These findings indicated the overlapped function in biological processes and pathways among three regions, which was further verified by the more detailed function analyses by comparing the function enrichments of the overall DE mRNAs in three regions (Supplementary Figure 4).

FIGURE 6

Figure 6. KEGG pathway analysis of the differentially up- and downregulated mRNAs in the SC, ACC, and AMY after nerve injury. (A–C) The top 10 enrichments of KEGG pathways from upregulated mRNAs (red panels on the left) and downregulated mRNAs (blue panels on the right) in the SC (A), ACC (B), and AMY (C) on day 7 post-SNL according to the P-value. The DAVID database was used to do the KEGG pathway analysis. Red and blue bars represent up- and downregulated mRNA enrichments, respectively.

PPI Network Establishment to Analyze Protein–Protein Interactions in Three Regions After SNL

To gain insight into the functional connection among the DE mRNAs and their potential role in neuropathic pain, a PPI network was conducted using the STRING database. The top 50 protein-coding DEGs with the highest correlation degree in each region were screened out and used to generate the network. As shown in Figure 7A, the increased DEGs, such as Itgam, Itgax, Tyrobp, Ptprc, Cd14, Fcgr3, and Cd44, were the crucial molecules among the hub genes in the network of SC, whereas the increased DEGs, such as Trp53, Mapk1, Mapk3, Fn1, Cnb3, Ube2c, Ube2d3, Fbxl19, Cdc34, Keap1, and Lmo7, and the decreased DEGs, such as Atg7, Socs1, Lnx1, Nedd4l, and Uba7, played a major role in the network of ACC (Figure 7C). The hub genes in the network of AMY (Figure 7E) revealed the vital position of Trp53, Mapk14, Tnf, Icam1, Itgan, Casp3, Myc, and Cd44. In addition, considering that neuroinflammation and apoptosis participated in many pathological processes including neurological and psychiatric disorders (Chelyshev et al., 2001; Ji et al., 2018; Matsuda et al., 2019), we defined the function of top 50 DEGs in the network of SC, ACC, and AMY by comparing them with a total number of 3,773 genes related to neuroinflammation (inflammation and immunity) and apoptosis (1,543, 3,348, and 1,279 related genes for inflammation, apoptosis, and immunity, respectively). Venn diagrams showed that approximately 82%, 86%, and 82% of DEGs, respectively, from the SC, ACC, and AMY were mapped to neuroinflammation- and apoptosis-related genes (Figures 7B,D,F).

FIGURE 7

Figure 7. PPI network establishments to analyze protein–protein interactions. (A,C,E) Top 50 differentially expressed genes (DEGs) were picked out based on the connection degree of genes and constructed the network in the SC (A), ACC (C), and AMY (E). The size of the node represents the connection degree that indicates the importance of the gene in the network. Red and blue colors represent the up- and downregulated genes, respectively. (B,D,F) Venn diagrams indicated the number and proportion of the selected top 50 DEGs mapped to apoptosis-, inflammation-, and immunity-related genes in the SC (B), ACC (D), and AMY (E), respectively.

Differentially Expressed mRNAs Implicated in Pain, Anxiety, and Depression Disorders Across Three Regions After SNL

We next used the Gene Cards database and CTD database to characterize the DEGs involved in pain and emotional disorders in the SC, ACC, and AMY. Based on the relevance score, about 230 (196 upregulated, 34 downregulated), 157 (100 upregulated, 57 downregulated), and 149 (100 upregulated, 49 downregulated) pain-related DEGs in the SC, ACC, and AMY, respectively, were observed (Supplementary Figure 5 and Supplementary Material 4). The ACC and AMY contained approximately 220 (140 upregulated, 80 downregulated) and 201 (146 upregulated, 55 downregulated) anxiety-related DEGs as well as about 278 (176 upregulated, 102 downregulated) and 287 (207 upregulated, 80 downregulated) depression-related DEGs, respectively (Supplementary Figure 5 and Supplementary Material 5). The top 20 highest up- and downregulated pain-related DEGs were displayed in heatmaps (Figures 8A–C).

FIGURE 8

Figure 8. Comparisons of pain-, anxiety-, and depression-related genes in the SC, ACC, and AMY. (A–C) Heatmaps of representative top 20 up- and downregulated pain-related DEGs in the SC (A), ACC (B), and AMY (C) on day 7 after SNL. (D) Venn diagram indicated the number of the overlapped pain-related genes among the SC, ACC, and AMY on day 7 after SNL. Detailed information was displayed in Supplementary Table 3. (E–G) Venn diagrams indicated the number of overlapped pain-, anxiety-, and depression-related DEGs in the ACC and AMY (detailed information in Supplementary Table 4). Colors in the heatmaps indicate the Row Z-score among the different datasets. Up- and downregulated genes are colored in red and green, respectively.

Pain-, Anxiety-, and Depression-Related DEGs Displayed a High Correlation to Neuroinflammation and Apoptosis in Three Regions After SNL

To further determine the function of DEGs in SNL-induced neuropathic pain pathogenesis, we investigated the correlation between the DEGs implicated in pain, depression, and anxiety disorders and the genes related to neuroinflammation (inflammation and immunity) and apoptosis. About 230 pain-related genes in the SC, 157 pain-related genes, 220 anxiety-related genes, and 278 depression-related genes in the ACC and 149 pain-related genes, 201 anxiety-related genes, and 287 depression-related genes in the AMY were mapped to the neuroinflammation- (inflammation and immunity) and apoptosis-related genes. There were about 51%, 60%, and 71% of pain-related DEGs in the SC; 42%, 52%, and 58% of pain-related DEGs in the ACC; and 47%, 59%, and 68% of pain-related DEGs in the AMY mapping to the datasets of apoptosis, inflammation, and immunity, respectively (Supplementary Table 1). Approximately 179, 99, and 110 pain-related overlapping DEGs were seen in three distinct regions, respectively (Supplementary Table 1). About 26%, 30%, and 36% of anxiety-related DEGs in the ACC; 31%, 36%, and 49% of anxiety-related DEGs in the AMY; 28%, 30%, and 41% of depression-related DEGs in the ACC; and 27%, 30%, and 44% of depression-related DEGs in the AMY were mapped to the genes associated with apoptosis, inflammation, and immunity, respectively (Supplementary Table 2). Furthermore, we found about 88 and 106 overlapping anxiety-related DEGs and 125 and 139 overlapping depression-related DEGs in the ACC and AMY, respectively (Supplementary Table 2).

Comparisons of Pain-, Anxiety-, and Depression-Related DEGs Among the SC, ACC, and AMY After SNL

To obtain more information about gene adaptations in response to SNL, we compared pain-related DEGs among three distinct regions in the present study (Figure 8D). Strikingly, about 103 overlapping genes were seen among the SC, ACC, and AMY (Supplementary Table 3). Of them, we observed the increased expression of C3, Cacna1c, Casp3, Cfh, Crem, Fn1, Men1, Oprl1, Sparc, Trp53, and Txnip and the decreased expression of Scn1a in all three regions. Interestingly, Scn1a was also reported downregulated in injured DRGs (Wu et al., 2016). In addition to the concordant expression patterns, we saw the discrepancy in the expression of the genes such as Alad, Ank1, Aurka, Ccnd3, Hgf, Mtm1, Nrf1, Tcf4, Wnk1, and Tardbp among the SC, ACC, and AMY. Consistent with the earlier findings in the injured DRGs (Wu et al., 2016), the expression of Gabra1 was reduced and the expression of Cacna2d1 was elevated, both in ACC and AMY after SNL. We also detected the upregulation of Acadvl, Acp5, Aifm1, Crh, Crlf1, Cxcl10, Deaf1, Eif4g1, Elane, F13a, Fgfr1, Fgfr3, Fhit, Gnas, Icam1, Il15, Il31ra, Itgam, Lgals1, Socs3, Stim1, Trappc2, Tsc22d3, Ube3a, Vgf, Vim, and Vip and the downregulation of Kif1a and Lipe in both SC and AMY. Finally, there were also many overlapped genes between the SC and ACC including the upregulated expression of Ikzf1, Pdyn, Plaur, Ppp1r1b, Pygl, and Sh2b3 and the downregulated expression of Anxa4, Cacna1b, Hrh1, Nf2, and Sgk1, as well as the inconsistent expression patterns in Arnt, Capn3, Ccr6, Cnbp, Igf1, Mapk8, Nfkb1, Id2, Litaf, Runx2, Tpm1, and Trps1.

Moreover, pain-, anxiety-, and depression-related DEGs in the ACC and AMY were compared (Figures 8E–G), A total of 338 pain-, anxiety-, and depression-related DEGs were observed in the ACC and AMY (Supplementary Table 4). Besides pain-related DEGs mentioned above, we saw the upregulated expression of Camk2b, Ccnd3, Ccng2, Cyp26b1, Fkbp5, Fgfr3, Gcnt2, Gnas, Gria1, Keap1, Map2, Mapk1, Mapk9, Nr4a2, Nr4a3, Oprk1, Oprl1, Pde4b, Pde7b, Ppfibp1, Psrc1, Tpm1, Tpm3, and Ube2c mRNAs and the downregulated expression of Tcf4, Pln, and Mtm1 mRNAs in both ACC and AMY. Additionally, the genes such as Cryab, Dclk1, Eif5a, Hgf, Magi1, Mef2c, Myc, Nr3c1, Plec Rai1, and Tardbp were downregulated in the ACC but upregulated in the AMY. The genes including Arrdc3, Diablo, Fmo5, Fnbp1, Ggt1, Nos1, Postn, Tcf7l2, and Xiap were significantly elevated in the ACC but reduced in the AMY. These findings suggested some shared pathogenesis mechanisms in neuropathic pain and emotional disorders, probably providing more evidence and basis for future researches.

Functional Prediction of DE lncRNAs in SNL

We next determined the function of lncRNAs through their related mRNAs by selecting the genes with the absolute value of correlation > 0.95 and the co-localization within 100 kb at the upstream and downstream (Supplementary Material 6).

GO enrichment analysis results were graphically displayed in directed acyclic graphs (DAGs), in which the branch represented the relationship of the inclusion that defined the smaller and smaller scales from top to bottom. The top 10 GO enrichments were selected as the master nodes of DAGs. They were shown together with the GO terms of containment relationships and systematically GO terms. DAGs were plotted from the biological process, molecular function, and cellular component aspects in the SC, ACC, and AMY, respectively (Figure 9 and Supplementary Figures 6, 7A-C). According to the distribution of predicted target genes in the Gene Ontology, the function of DE lncRNAs was clarified and displayed in the form of histograms by the -log (P-value) of each GO term. Apparently, the most significant biological process enrichments were synapse assembly, cell–cell adhesion via plasma membrane adhesion molecules, synapse organization, cellular macromolecule metabolic process, nucleic acid metabolic process, and regulation of RNA metabolic process in three regions. The noteworthy cellular component enrichments were seen in intracellular, MHC class I protein complex, nuclear lumen, membrane-enclosed lumen, and intracellular organelle lumen in the SC, ACC, and AMY. The most robust molecular functions were enriched in binding, nucleic acid binding, heterocyclic compound binding, organic cyclic compound binding, and TAP binding among three regions (Figure 9 and Supplementary Figures 6, 7D).

FIGURE 9

Figure 9. Functional prediction of DE lncRNAs by GO and KEGG analyses in the spinal cord of SNL mice. (A–C) Directed acyclic graphs (DAGs) graphically display the significant GO enrichment results with the candidate targeted genes in biological process (A), molecular function (B), and cellular component (C). (D) Significant molecular function, biological process, and cellular component enrichment analysis of DE lncRNA-related mRNAs. The enrichment scores (−log10 (P-value)) of the GO term were shown in the histogram. (E,F) The DE lncRNA-related mRNA-enriched KEGG pathways represented by enrichment scores (−log10 (P-value)) (E) and the scatterplot showing statistics of pathway enrichment (F), respectively. The color of pathway terms is defined by the enrichment P-value.

Similarly, the top 20 KEGG enrichments were shown in the histograms by the -log (P-value) of each pathway (Figure 9 and Supplementary Figures 6, 7E) and together with the enriched distribution maps (Figure 9 and Supplementary Figures 6, 7F), in which the degree of KEGG enrichment was assessed by the Rich factor, P-value, and the number of genes. The enrichment was more significant with the greater rich factor and the larger number of genes but the smaller P-value. The most significantly enriched pathways were related to cell adhesion molecules (CAMs), graft-versus-host disease, type 1 diabetes mellitus, antigen processing and presentation, allograft rejection autoimmune thyroid disease, and cellular senescence among three regions after SNL (Figure 9 and Supplementary Figures 6, 7F). Overall, these data demonstrated the overlapped effects in DE lncRNA function among three regions, similar to the patterns of the DE mRNAs.

lncRNA–mRNA Co-expression Network Analysis

To observe the potential interaction between lncRNAs and mRNAs in three regions after SNL, gene co-expression networks were constructed based on the correlation analysis. The networks were established by numbers of DE lncRNAs and the most potential top 50 DE mRNAs targets (PCC > 0.95 or <–0.95, and FDR < 0.05) (Supplementary Material 7). The cis-acting regulatory networks were constructed with 136 relationships between 88 lncRNAs (61 known, 26 predicted) and top 50 mRNAs in the SC (Figure 10A), 137 relationships between 87 lncRNAs (56 known, 31 predicted) and top 50 mRNAs in the ACC (Figure 10C), and 137 relationships between 89 lncRNAs (59 known, 30 predicted) and top 50 mRNAs in the AMY (Figure 10E). In contrast, the co-expression networks for trans-acting regulation consisted of 80 relationships between 17 lncRNAs (13 known, 4 predicted) and top 50 mRNAs in the SC (Figure 10B), 45 relationships between 13 lncRNAs (10 known, 3 predicted) and top 50 mRNAs in the ACC (Figure 10D), and 51 relationships between 17 lncRNAs (11 known, 6 predicted) and top 50 mRNAs in the AMY (Figure 10F).

FIGURE 10

Figure 10. Analysis of lncRNA–mRNA co-expression network in the SC, ACC, and AMY of SNL mice. (A,B) Co-expression network of cis-acting regulatory elements (A) and trans-acting regulatory elements (B) in the SC. (C,D) Co-expression network of cis-acting regulatory elements (C) and trans-acting regulatory elements (D) in the ACC. (E,F) The co-expression network of cis-acting regulatory elements (E) and trans-acting regulatory elements (F) in the AMY. LncRNA–mRNA co-expression network was constructed based on the Pearson correlation coefficients. Pink color nodes represent lncRNA, while blue nodes represent the co-expression mRNAs.

Discussion

Neuropathic pain is a somatosensory disorder resulting from nerve injury or diseases affecting the peripheral and central nervous systems (Colloca et al., 2017). With the high incidence and poor management in the clinic, it is a major public health problem. Over the past decades, the potential mechanisms underlying neuropathic pain have been extensively studied. However, the effective treatments are still limited due to the largely unknown molecular mechanisms (Finnerup et al., 2015; Jones et al., 2018). Evidence demonstrates that the alteration in gene expression profiles at different levels of the nervous system plays an important role in the development and maintenance of neuropathic pain. In the present study, we reported gene transcript alterations in pain- and emotion-associated regions in the central nervous system following peripheral nerve injury in mice using the next-generation RNA sequencing assay. Bioinformatics and pathway analyses revealed that particular differentially expressed gene patterns and biological networks in the SC, ACC, and AMY were in correspondence with SNL-induced nociceptive hypersensitivities and pain-related aversion.

In the present study, numerous DEGs belonging to GPCR and ion channel genes were identified. Many of these DEGs have identified function in pain or emotion dysfunction. Ccr1 and Ccr5 were shown to be involved in heat hyperalgesia in mice (Eltayeb et al., 2007; Pevida et al., 2014). The modification of Cmklr1 was shown to implicate depression in the prefrontal cortex and hippocampus responding to chronic restraint stress (CRS) (Guo et al., 2012; Doyle et al., 2014). Genetic knockout of Orai1 nearly eliminated the second phase of formalin-induced pain and attenuated carrageenan-induced pain hypersensitivity and neuronal excitability (Dou et al., 2018). TRPC6 inhibition in the spinal cord blocked the induction of morphine tolerance and hyperalgesia in rats (Jin et al., 2017; Wang et al., 2020). Likewise, we identified the changes in lncRNA expression after SNL. The expression of lncRNA Malat1 was controversial. Following our sequencing result, Meng C et al. demonstrated that the inhibition of spinal Malat1 expression contributed to neuropathic pain after brachial plexus avulsion (Meng et al., 2019). On the contrary, another study reported that the inhibition of spinal Malat1 reduced the incidence of CCI-induced neuropathic pain (Wu et al., 2020). These findings suggest that the molecular mechanism underlying neuropathic pain may vary with different etiologies and courses. LncRNA H19 was upregulated in the injured DRGs and hippocampus neurons following peripheral nerve injury (Han et al., 2020; Wen et al., 2020). Consistently, our sequencing data showed that SNL increased its expression in the spinal cord and ACC. However, whether the increased H19 in these two regions contributes to neuropathic pain needs to be confirmed.

Furthermore, lots of pain-, anxiety-, and depression-related genes were identified in the SC, ACC, and AMY following SNL. Consistent with previous studies, Itgam (CD11b) (Ghasemlou et al., 2015), Tlr3 (Liu et al., 2012), Bdnf (Sapio et al., 2019), and Stim1 (Gao et al., 2016) were significantly elevated, while Kif1a (Wang et al., 2018) and Rbl2 (Chen et al., 2019) were robustly decreased in the spinal cord. However, there was a discrepancy in the expression of Gsk3b. As reported, Gsk3b was downregulated at day 3 and upregulated at day 10 in the spinal cord after partial sciatic nerve ligation (Weng et al., 2014). Unexpectedly, we saw the decreased expression of spinal Gsk3b on day 7 after SNL. These results may imply that spinal Gsk3b expression is time-dependent post-nerve injury. In the ACC and AMY, we observed the upregulation of Drd1 (Darvish-Ghane et al., 2020), Tnf (Akter et al., 2020), and Il33 (Fairlie-Clarke et al., 2018), as well as the increased expression of Runx1 and Cd68 among pain-related DEGs. However, the function of Runx1 and Cd68 in these two regions is still unknown and remains to be investigated.

According to the sequencing data from the injured DRGs (Wu et al., 2016), we observed similar expression changes of some genes such as the upregulation of Atf3, Ccr1, Ccr5, and Gal, and the downregulation of Gria2 in the spinal cord as well as the elevated expression of Cacna2d1 and the reduced expression of Gabra1 in the ACC and AMY after SNL. Besides, we found the consistent expression patterns of some genes such as the upregulation of Atf3, Bdnf, C1qa, C3, Cck, Ccl2, Cd68, Csf1, Cx3cr1, Gch1, Itgam, Ngfr, Pdyn, and Tlr2 in the spinal cord of SNL mice, when compared to the sequencing results in CCI rats (Du et al., 2018; Korczeniewska et al., 2020). These data suggest the common gene expression patterns independent of regions, models, or species. The ACC and AMY are important brain areas in pain and emotion modulation, but the gene expression profiles in these two regions after nerve injury have not been identified before. We achieved more anxiety- and depression-related DEGs than pain-related DEGs in these two regions, consistent with their significant function in emotion processing. Interestingly, among the anxiety- and depression-related DEGs, we found the upregulation of Bcl3, C3, Gpat3, and Tnf in the AMY as well as the increased expression of Crhr2, Nant, Sar1a, and Tgif1 and the decreased expression of Ier2, Il12a, Nrep, and Tnfrsf25 in the ACC. These alternations were in step with the sequencing results in 12- and 24-month-old mice (Li M. et al., 2020), suggesting some common gene expression alterations in different pathological processes, at least for anxiety- and depression-related genes.

Consistent with previous reports (Jiang et al., 2015; Wu et al., 2016), GO term and KEGG pathway enrichment analyses in three regions showed notable enrichments in apoptotic, inflammation, immunity, cytokine production and defense response, behavior, and sensory perception of pain as well as the enrichments in chemokine signaling pathway, MAPK signaling pathway, TNF signaling pathway, cAMP signaling pathway, Type II diabetes mellitus, and T cell receptor signaling pathway. Consistently, functional analyses observed that large percentages of pain-, anxiety-, and depression-related DEGs were highly related to neuroinflammation and apoptosis that were considered to occupy an important position in pain states (Chelyshev et al., 2001; Ji et al., 2018; Matsuda et al., 2019). Among the overlapped pain-related DEGs, the amounts of Cacna1c, Casp3, C3, and TXNIP were sharply elevated in all three regions following SNL. ACC-conditional deletion of Cav1.2 channels impaired observational fear learning and reduced behavioral pain responses, while neuronal deletion of Cav1.2 led to significant deficits in the extinction of conditioned fear and altered sIPSC and sEPSC activity within the amygdala (Jeon et al., 2010; Temme and Murphy, 2017). Evidence indicated that Casp3 was upregulated in neuropathic pain and that the activation of Casp3 was required in long-term depression (Li et al., 2010; Yang et al., 2018). The deletion of complement C3 was shown to reduce pain-, anxiety-, and depression-like behaviors and to improve learning and spatial memory in aged mice (Shi et al., 2015; Crider et al., 2018). Recent reports suggested that activation of TXNIP/NLRP3 axis was positively associated with pain and emotion disorders and the neuroprotective properties by pharmacological inhibition or genetic deletion of TXNIP following cerebrovascular and neurodegenerative diseases (Nasoohi et al., 2018; Pan et al., 2018). Taken together, the remarkable overlapped DEGs might be the most potential candidates for the researches on pain-, anxiety-, or depression-disorders.

Despite that we reported the unique transcriptome profiles and conducted a series of functional analyses, the present study still has some limitations. Firstly, it should be noted that the SC, ACC, and AMY contain a variety of cell populations including different types of neurons, astrocytes, and microglia cells. However, all bioinformatics analyses presented in this work were obtained from all cell populations. Thus, future studies on the cell-type-specific changes in gene expression following peripheral nerve injury should be performed using single-cell sequencing analysis. Secondly, we used the database of anxiety- and depression-related genes to analyze the DEGs in the ACC and AMY after SNL. However, six days post-SNL, anxiety- and depression-like behaviors were not completely developed even if SNL-induced aversion was detected (Suzuki et al., 2007; Wu et al., 2019). However, the evidence demonstrated that the analyzed anxiety- and depression-related DEGs such as Gria1 (Rivera et al., 2020), Mapk1 (Sierra-Fonseca et al., 2019), Mapk9 (Thomson et al., 2020), and Fkbp5 (Zannas et al., 2019) contributed to anxiety or depression symptoms in rodents. The further study on the gene expression profiles and emotion-related behaviors including anxiety- and depression-like behaviors at the later stage of neuropathic pain should be carried out. Thirdly, investigations on the gender-specific and age-specific pain mechanisms should also be included due to the more frequent prevalence of pain in women and aged patients (Bouhassira et al., 2008). Fourthly, other brain regions such as the medial prefrontal cortex, nucleus accumbens, and periaqueductal gray were reported to participate in pain- and emotion-related behaviors under the conditions of chronic stress and/or chronic pain as well (Bouhassira et al., 2008; Descalzi et al., 2017; Smith et al., 2021). To obtain more valuable information, RNA sequencing analysis at these brain regions needs to be considered in the future. Finally, although the present study demonstrated gene transcript alternations and their functional analyses in the SC, ACC, and AMY, whether these changes contribute to the induction and maintenance of neuropathic pain and whether they can serve as new targets remain to be further determined.

Conclusion

In summary, we for the first time provided the unique gene expression profiles of lncRNAs and mRNAs in three pain-related regions and revealed the implication of neuroinflammation and apoptosis in the pathogenesis of neuropathic pain using different bioinformatics analyses. The comparisons of RNA sequencing results provide a more thorough analysis of gene expression alterations in three distinct pain-related regions. Overall, our findings present comprehensive information that may facilitate the discovery of novel analgesic strategies.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: National Center for Biotechnology Information (NCBI) BioProject, https://www.ncbi.nlm.nih.gov/bioproject/, PRJNA705299.

Ethics Statement

The animal study was reviewed and approved by the Animal Care and Use Committee of Zhengzhou University.

Author Contributions

WZ conceived the project, supervised all the experiments, and edited the manuscript. SS, ML, DW, JC, XR, Y-XT, and WZ assisted with experimental design. SS and ML carried out behavioral tests, surgery, and tissue collection. SS performed the RT-PCR assay and wrote the draft of the manuscript. SS, ML, and DW analyzed the data. Y-XT and WZ edited the manuscript. All authors read and discussed the manuscript.

Funding

This work was supported by the start-up fund from the Zhengzhou University and the National Natural Science Foundation of China (Grant Numbers 81471144 and 81771195).

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.

Acknowledgments

The authors thank Dr. Xiang Gao and Mrs. Tingting Wang for their supports in RNA sequencing. The authors also thank Mrs. Peijun Jia and Mrs. Yake Wang for their assistance in editing the manuscript.

Supplementary Material

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

Supplementary Figure 1 | Transcriptome profiling in the SC, ACC and AMY in mice after nerve injury. (A) The mapped proportions of exonic, intronic, and intergenic reads in the SC, ACC and AMY on day 7 post-SNL or sham surgery. (B-D) The proportions of reads that align to exonic, intronic and intergenic regions from the SC (B), ACC (C), and AMY (D) on day 7 post-surgery. n = 3 biological repeats. ^∗P < 0.05; ^∗∗P < 0.01 versus the corresponding sham group by two-tailed unpaired Student’s t-test. (E-G) Distribution of differentially expressed RNAs in SC (E), ACC (F), and AMY (G) on day 6 after SNL.

Supplementary Figure 2 | Heatmaps of the representative differentially expressed genes (DEGs) in the ACC after nerve injury. Top 15 up-and down-regulated G protein-coupled receptor mRNAs (A), ion channel mRNAs (B), and lncRNAs (C) in the ACC on day 7 after SNL. Colors in the heatmaps indicate the Row Z-score among the different datasets. The up- and down-regulated genes are colored in red and green, respectively.

Supplementary Figure 3 | Heatmaps of the representative differentially expressed genes (DEGs) in the AMY after nerve injury. Top 15 up- and down-regulated G protein-coupled receptor mRNAs (A), ion channel mRNAs (B), and lncRNAs (C) in the AMY on day 7 after SNL. Colors in the heatmaps indicate the Row Z-score among the different datasets. The up- and down-regulated genes are colored in red and green, respectively.

Supplementary Figure 4 | The common biological processes and KEGG pathways in the SC, ACC and AMY after nerve injury. (A-D) The top 10 common enrichments in biological processes (blue panels on the left) and KEGG pathways (red panels on the right) between SC and ACC (A), SC and AMY (B), or ACC and AMY (C). The top 10 common function enrichments in biological processes (blue panels on the left) and KEGG pathways (red panels on the right) among SC, ACC and AMY (D) according to the P-value. Blue and red bars represent the biological process and pathway enrichments, respectively.

Supplementary Figure 5 | Heatmaps of the identified pain-, anxiety- and depression-related genes in the SC, ACC and AMY after nerve injury. (A) Heatmaps of 157 (100 upregulated and 57 downregulated) pain-, 220 (140 upregulated and 80 downregulated) anxiety- and 278 (176 upregulated and 102 downregulated) depression-related DEGs in the ACC. (B) Heatmaps of 149 (100 upregulated and 49 downregulated) pain-, 201 (146 upregulated and 55 downregulated) anxiety- and 287 (207 upregulated and 80 downregulated) depression-related DEGs in the AMY; C, Heatmaps of 230 (196 upregulated and 34 downregulated) pain-related DEGs in the spinal cord. Colors in the heatmaps indicate the Row Z-score among the different datasets. The up- and down-regulated genes are colored in red and green, respectively.

Supplementary Figure 6 | The functional prediction of DE lncRNAs by GO and KEGG analyses in the ACC from SNL mice. (A-C) Directed Acyclic Graphs (DAGs) graphically display the significant GO enrichment results with the candidate targeted genes in biological process (A), molecular function (B), and cellular component (C). (D) The significant molecular function, biological process and cellular component enrichment analysis of DE lncRNA-related mRNAs. (E,F) The DE lncRNA-related mRNAs enriched KEGG pathways represented by enrichment scores (−log10 (P-value)) (E) and the scatterplot showing statistics of pathway enrichment (F), respectively. The color of pathway terms is defined by the enrichment P-value.

Supplementary Figure 7 | The functional prediction of DE lncRNAs by GO and KEGG analyses in the AMY from SNL mice. (A-C) Directed Acyclic Graphs (DAGs) graphically display the significant GO enrichment results with the candidate targeted genes in biological process (A), molecular function (B), and cellular component (C). (D) The significant molecular function, biological process and cellular component enrichment analysis of DE lncRNA-related mRNAs. (E,F) The DE lncRNA-related mRNAs enriched KEGG pathways represented by enrichment scores (−log10 (P-value)) (E) and the scatterplot showing statistics of pathway enrichment (F), respectively. The color of pathway terms is defined by the enrichment P-value.

Footnotes

1. ^ http://www.omicshare.com/tools
2. ^ https://www.genecards.org/
3. ^ http://ctdbase.org/
4. ^ https://bioinfogp.cnb.csic.es/tools/venny/index.html
5. ^ http://david.abcc.ncifcrf.gov/
6. ^ http://www.geneontology.org
7. ^ http://www.genome.jp/kegg/
8. ^ https://string-db.org/cgi/
9. ^ www.cytoscape.org/

References

Akter, S., Uddin, K. R., Sasaki, H., and Shibata, S. (2020). Gamma oryzanol alleviates high-fat diet-induced anxiety-like behaviors through downregulation of dopamine and inflammation in the amygdala of mice. Front. Pharmacol. 11:330. doi: 10.3389/fphar.2020.00330

PubMed Abstract | CrossRef Full Text | Google Scholar

Barcelon, E. E., Cho, W. H., Jun, S. B., and Lee, S. J. (2019). Brain microglial activation in chronic pain-associated affective disorder. Front. Neurosci. 13:213. doi: 10.3389/fnins.2019.00213

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouhassira, D., Lantéri-Minet, M., Attal, N., Laurent, B., and Touboul, C. (2008). Prevalence of chronic pain with neuropathic characteristics in the general population. Pain 136, 380–387. doi: 10.1016/j.pain.2007.08.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowers, M. E., and Ressler, K. J. (2015). Interaction between the cholecystokinin and endogenous cannabinoid systems in cued fear expression and extinction retention. Neuropsychopharmacology 40, 688–700. doi: 10.1038/npp.2014.225

PubMed Abstract | CrossRef Full Text | Google Scholar

Bushnell, M. C., Case, L., Ceko, M., Cotton, V. A., Gracely, J. L., Low, L. A., et al. (2015). Effect of environment on the long- term consequences of chronic pain. Pain 156, S42–S49. doi: 10.1097/01.j.pain.0000460347.77341.bd

CrossRef Full Text | Google Scholar

Campbell, A. P., and Smrcka, A. V. (2018). Targeting G protein-coupled receptor signalling by blocking G proteins. Nat. Rev. Drug Discovery 17, 789–803. doi: 10.1038/nrd.2018.135

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, H., Gao, Y. J., Ren, W. H., Li, T. T., Duan, K. Z., Cui, Y. H., et al. (2009). Activation of extracellular signal-regulated kinase in the anterior cingulate cortex contributes to the induction and expression of affective pain. J. Neurosci. 29, 3307–3321. doi: 10.1523/JNEUROSCI.4300-08.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Carpanini, S. M., Torvell, M., and Morgan, B. P. (2019). Therapeutic inhibition of the complement system in diseases of the central nervous system. Front. Immunol. 10:362. doi: 10.3389/fimmu.2019.00362

PubMed Abstract | CrossRef Full Text | Google Scholar

Chandley, M. J., Crawford, J. D., Szebeni, A., Szebeni, K., and Ordway, G. A. (2015). NTRK2 expression levels are reduced in laser captured pyramidal neurons from the anterior cingulate cortex in males with autism spectrum disorder. Mol. Autism 6:28. doi: 10.1186/s13229-015-0023-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Che, X., Deng, X., Xie, K., Wang, Q., Yan, J., Shao, X., et al. (2019). Long noncoding RNA MEG3 suppresses podocyte injury in diabetic nephropathy by inactivating Wnt/β-catenin signaling. PeerJ 7:e8016. doi: 10.7717/peerj.8016

PubMed Abstract | CrossRef Full Text | Google Scholar

Chelyshev, I., Cherepnev, G. V., and Saĭtkulov, K. I. (2001). Apoptoz v nervnoĭ sisteme [Apoptosis in the nervous system]. Ontogenez 32, 118–129.

Google Scholar

Chen, J., Li, L., Chen, S. R., Chen, H., Xie, J. D., Sirrieh, R. E., et al. (2018). The α2δ-1-NMDA receptor complex is critically involved in neuropathic pain development and gabapentin therapeutic actions. Cell Rep. 22, 2307–2321. doi: 10.1016/j.celrep.2018.02.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Q. B., Li, Z. H., Fu, Y., Lv, N. N., Tian, N., Han, L., et al. (2019). Downregulated long non-coding RNA LINC00899 inhibits invasion and migration of spinal ependymoma cells via RBL2-dependent FoxO pathway. Cell Cycle (Georgetown, Tex.) 18, 2566–2579. doi: 10.1080/15384101.2019.1652046

PubMed Abstract | CrossRef Full Text | Google Scholar

Colloca, L., Ludman, T., Bouhassira, D., Baron, R., Dickenson, A. H., Yarnitsky, D., et al. (2017). Neuropathic pain. Nat. Rev. Dis. Primers 3:17002. doi: 10.1038/nrdp.2017.2

PubMed Abstract | CrossRef Full Text | Google Scholar

Cosi, C., Mannaioni, G., Cozzi, A., Carlà, V., Sili, M., and Cavone, L. (2011). G-protein coupled receptor 35 (GPR35) activation and inflammatory pain: studies on the antinociceptive effects of kynurenic acid and zaprinast. Neuropharmacology 60, 1227–1231. doi: 10.1016/j.neuropharm.2010.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Crider, A., Feng, T., Pandya, C. D., Davis, T., Nair, A., Ahmed, A. O., et al. (2018). Complement component 3a receptor deficiency attenuates chronic stress-induced monocyte infiltration and depressive-like behavior. Brain Behav. Immun. 70, 246–256. doi: 10.1016/j.bbi.2018.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Darvish-Ghane, S., Quintana, C., Beaulieu, J. M., and Darvish (2020). D1 receptors in the anterior cingulate cortex modulate basal mechanical sensitivity threshold and glutamatergic synaptic transmission. Mol. Brain 13:121. doi: 10.1186/s13041-020-00661-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Descalzi, G., Mitsi, V., Purushothaman, I., Gaspari, S., Avrampou, K., Loh, Y. E., et al. (2017). Neuropathic pain promotes adaptive changes in gene expression in brain networks involved in stress and depression. Sci. Signal. 10:eaaj1549. doi: 10.1126/scisignal.aaj1549

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, R., Jia, D., Xue, P., Cui, X., Li, K., Zheng, S., et al. (2014). Genome-wide analysis of long noncoding RNA (lncRNA) expression in hepatoblastoma tissues. PLoS One 9:e85599. doi: 10.1371/journal.pone.0085599

PubMed Abstract | CrossRef Full Text | Google Scholar

Dou, Y., Xia, J., Gao, R., Gao, X., Munoz, F. M., Wei, D., et al. (2018). Orai1 plays a crucial role in central sensitization by modulating neuronal excitability. J. Neurosci. 38, 887–900. doi: 10.1523/JNEUROSCI.3007-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Doyle, J. R., Krishnaji, S. T., Zhu, G., Xu, Z. Z., Heller, D., Ji, R. R., et al. (2014). Development of a membrane-anchored chemerin receptor agonist as a novel modulator of allergic airway inflammation and neuropathic pain. J. Biol. Chem. 289, 13385–13396. doi: 10.1074/jbc.M113.522680

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, H., Shi, J., Wang, M., An, S., Guo, X., Wang, Z., et al. (2018). Analyses of gene expression profiles in the rat dorsal horn of the spinal cord using RNA sequencing in chronic constriction injury rats. J. Neuroinflamm. 15:280. doi: 10.1186/s12974-018-1316-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Eltayeb, S., Berg, A. L., Lassmann, H., Wallström, E., Nilsson, M., Olsson, T., et al. (2007). Temporal expression and cellular origin of CC chemokine receptors CCR1, CCR2 and CCR5 in the central nervous system: insight into mechanisms of MOG-induced EAE. J. Neuroinflamm. 4:14. doi: 10.1186/1742-2094-4-14

PubMed Abstract | CrossRef Full Text | Google Scholar

Fairlie-Clarke, K., Barbour, M., Wilson, C., Hridi, S. U., Allan, D., and Jiang, H. R. (2018). Expression and function of IL-33/ST2 axis in the central nervous system under normal and diseased conditions. Front. Immunol. 9:2596. doi: 10.3389/fimmu.2018.02596

PubMed Abstract | CrossRef Full Text | Google Scholar

Finnerup, N. B., Attal, N., Haroutounian, S., McNicol, E., Baron, R., Dworkin, R. H., et al. (2015). Pharmacotherapy for neuropathic pain in adults: a systematic review and meta-analysis. Lancet. Neurol. 14, 162–173. doi: 10.1016/S1474-4422(14)70251-0

CrossRef Full Text | Google Scholar

Gangarossa, G., Laffray, S., Bourinet, E., and Valjent, E. (2014). T-type calcium channel Cav3.2 deficient mice show elevated anxiety, impaired memory and reduced sensitivity to psychostimulants. Front. Behav. Neurosci. 8:92. doi: 10.3389/fnbeh.2014.00092

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, X., Xia, J., Munoz, F. M., Manners, M. T., Pan, R., Meucci, O., et al. (2016). STIMs and Orai1 regulate cytokine production in spinal astrocytes. J. Neuroinflamm. 13:126. doi: 10.1186/s12974-016-0594-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y. J., Ren, W. H., Zhang, Y. Q., and Zhao, Z. Q. (2004). Contributions of the anterior cingulate cortex and amygdala to pain- and fear-conditioned place avoidance in rats. Pain 110, 343–353. doi: 10.1016/j.pain.2004.04.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghasemlou, N., Chiu, I. M., Julien, J. P., and Woolf, C. J. (2015). CD11b+Ly6G- myeloid cells mediate mechanical inflammatory pain hypersensitivity. Proc. Natl. Acad. Sci. U S A. 112, E6808–E6817. doi: 10.1073/pnas.1501372112

PubMed Abstract | CrossRef Full Text | Google Scholar

González de San Román, E., Manuel, I., Ledent, C., Chun, J., and Rodr et al., (2019). CB1 and LPA1 receptors relationship in the mouse central nervous system. Front. Mol. Neurosci. 12:223. doi: 10.3389/fnmol.2019.00223

PubMed Abstract | CrossRef Full Text | Google Scholar

Guilloux, J. P., Douillard-Guilloux, G., Kota, R., Wang, X., Gardier, A. M., Martinowich, K., et al. (2012). Molecular evidence for BDNF- and GABA-related dysfunctions in the amygdala of female subjects with major depression. Mol. Psychiatry 17, 1130–1142. doi: 10.1038/mp.2011.113

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, W., Chu, Y.-X., Imai, S., Yang, J. L., Zou, S., Mohammad, Z., et al. (2016). Further observations on the behavioral and neural effects of bone marrow stromal cells in rodent pain models. Mol. Pain 12:1744806916658043. doi: 10.1177/1744806916658043

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, X., Fu, Y., Xu, Y., Weng, S., Liu, D., Cui, D., et al. (2012). Chronic mild restraint stress rats decreased CMKLR1 expression in distinct brain region. Neurosci. Lett. 524, 25–29. doi: 10.1016/j.neulet.2012.06.075

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, C. L., Liu, Y. P., Guo, C. J., Du, T. T., Jiang, Y., and Wang, K. L. (2020). The lncRNA H19 binding to let-7b promotes hippocampal glial cell activation and epileptic seizures by targeting Stat3 in a rat model of temporal lobe epilepsy. Cell Proliferation 53:e12856. doi: 10.1111/cpr.12856

PubMed Abstract | CrossRef Full Text | Google Scholar

Horváth, G., Gölöncsér, F., Csölle, C., Király, K., Andó, R. D., Baranyi, M., et al. (2014). Central P2Y12 receptor blockade alleviates inflammatory and neuropathic pain and cytokine production in rodents. Neurobiol. Dis. 70, 162–178. doi: 10.1016/j.nbd.2014.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Hrdlickova, R., Toloue, M., and Tian, B. (2017). RNA-Seq methods for transcriptome analysis. Wiley interdisciplinary reviews. RNA 8:10.1002/wrna.1364. doi: 10.1002/wrna.1364

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C., Lu, F., Li, P., Cao, C., and Liu, Z. (2017). Tlx3 function in the dorsal root ganglion is pivotal to itch and pain sensations. Front. Mol. Neurosci. 10:205. doi: 10.3389/fnmol.2017.00205

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeon, D., Kim, S., Chetana, M., Jo, D., Ruley, H. E., Lin, S. Y., et al. (2010). Observational fear learning involves affective pain system and Cav1.2 Ca2+ channels in ACC. Nat. Neurosci. 13, 482–488. doi: 10.1038/nn.2504

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, G., Zhang, W., Mahimainathan, L., Narasimhan, M., Kiritoshi, T., Fan, X., et al. (2017). 5-HT2C receptor knockdown in the amygdala inhibits neuropathic-pain-related plasticity and behaviors. J. Neurosci. 37, 1378–1393. doi: 10.1523/JNEUROSCI.2468-16.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, R. R., Nackley, A., Huh, Y., Terrando, N., and Maixner, W. (2018). Neuroinflammation and central sensitization in chronic and widespread pain. Anesthesiology 129, 343–366. doi: 10.1097/ALN.0000000000002130

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, B. C., Sun, W. X., He, L. N., Cao, D. L., Zhang, Z. J., Gao, Y. J., et al. (2015). Identification of lncRNA expression profile in the spinal cord of mice following spinal nerve ligation-induced neuropathic pain. Mol. Pain 11:43. doi: 10.1186/s12990-015-0047-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, B. C., Zhang, W. W., Yang, T., Guo, C. Y., Cao, D. L., Zhang, Z. J., et al. (2018). Demethylation of G-Protein-Coupled Receptor 151 promoter facilitates the binding of krüppel-like factor 5 and enhances neuropathic pain after nerve injury in mice. J. Neurosci. 38, 10535–10551. doi: 10.1523/JNEUROSCI.0702-18.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, H., Sun, Y. T., Guo, G. Q., Chen, D. L., Li, Y. J., Xiao, G. P., et al. (2017). Spinal TRPC6 channels contributes to morphine-induced antinociceptive tolerance and hyperalgesia in rats. Neurosci. Lett. 639, 138–145. doi: 10.1016/j.neulet.2016.12.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, J. M. R., Viswanath, O., Peck, J., Kaye, A. D., Gill, J. S., Simopoulos, T. T., et al. (2018). A brief history of the opioid epidemic and strategies for pain medicine. Pain Therapy 7, 13–21. doi: 10.1007/s40122-018-0097-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Klein, M. C., and Oaklander, A. L. (2018). Ion channels and neuropathic pain. eLife 7:e42849. doi: 10.7554/eLife.42849

PubMed Abstract | CrossRef Full Text | Google Scholar

Korczeniewska, O. A., Katzmann Rider, G., Gajra, S., Narra, V., Ramavajla, V., Chang, Y. J., et al. (2020). Differential gene expression changes in the dorsal root versus trigeminal ganglia following peripheral nerve injury in rats. Eur J. Pain (London, England) 24, 967–982. doi: 10.1002/ejp.1546

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuner, R. (2010). Central mechanisms of pathological pain. Nat. Med. 16, 1258–1266. doi: 10.1038/nm.2231

PubMed Abstract | CrossRef Full Text | Google Scholar

LaCroix-Fralish, M. L., Austin, J. S., Zheng, F. Y., Levitin, D. J., and Mogil, J. S. (2011). Patterns of pain: meta-analysis of microarray studies of pain. Pain 152, 1888–1898. doi: 10.1016/j.pain.2011.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., Su, S., Cai, W., Cao, J., Miao, X., Zang, W., et al. (2020). Differentially expressed genes in the brain of aging mice with cognitive alteration and depression- and anxiety-like behaviors. Front. Cell Dev. Biol. 8:814. doi: 10.3389/fcell.2020.00814

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Tang, C., Wang, J., Guo, P., Wang, C., Wang, Y., et al. (2018). Methylene blue relieves the development of osteoarthritis by upregulating lncRNA MEG3. Exp. Therapeutic Med. 15, 3856–3864. doi: 10.3892/etm.2018.5918

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Guo, X., Sun, L., Xiao, J., Su, S., Du, S., et al. (2020). N6-Methyladenosine demethylase FTO contributes to neuropathic pain by stabilizing g9a expression in primary sensory neurons. Adv. Sci. (Weinheim, Baden-Wurttemberg, Germany) 7:1902402. doi: 10.1002/advs.201902402

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Jo, J., Jia, J. M., Lo, S. C., Whitcomb, D. J., Jiao, S., et al. (2010). Caspase-3 activation via mitochondria is required for long-term depression and AMPA receptor internalization. Cell 141, 859–871. doi: 10.1016/j.cell.2010.03.053

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, K., Fu, Q., Liu, Y., Cao, J., Miao, X., Zang, W., et al. (2019). An integrative bioinformatics analysis of microarray data for identifying hub genes as diagnostic biomarkers of preeclampsia. Biosci. Rep. 39:BSR20190187. doi: 10.1042/BSR20190187

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, T., Berta, T., Xu, Z. Z., Park, C. K., Zhang, L., Lü, N., et al. (2012). TLR3 deficiency impairs spinal cord synaptic transmission, central sensitization, and pruritus in mice. J. Clin. Invest. 122, 2195–2207. doi: 10.1172/JCI45414

PubMed Abstract | CrossRef Full Text | Google Scholar

Maftei, D., Marconi, V., Florenzano, F., Giancotti, L. A., Castelli, M., Moretti, S., et al. (2014). Controlling the activation of the Bv8/prokineticin system reduces neuroinflammation and abolishes thermal and tactile hyperalgesia in neuropathic animals. Br. J. Pharmacol. 171, 4850–4865. doi: 10.1111/bph.12793

PubMed Abstract | CrossRef Full Text | Google Scholar

Manville, R. W., and Abbott, G. W. (2018). Gabapentin is a potent activator of KCNQ3 and KCNQ5 potassium channels. Mol. Pharmacol. 94, 1155–1163. doi: 10.1124/mol.118.112953

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsuda, M., Huh, Y., and Ji, R. R. (2019). Roles of inflammation, neurogenic inflammation, and neuroinflammation in pain. J. Anesthesia 33, 131–139. doi: 10.1007/s00540-018-2579-4

PubMed Abstract | CrossRef Full Text | Google Scholar

McKinney, B. C., Sze, W., Lee, B., and Murphy, G. G. (2009). Impaired long-term potentiation and enhanced neuronal excitability in the amygdala of Ca(V)1.3 knockout mice. Neurobiol. Learn. Memory 92, 519–528. doi: 10.1016/j.nlm.2009.06.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, C., Yang, X., Liu, Y., Zhou, Y., Rui, J., Li, S., et al. (2019). Decreased expression of lncRNA Malat1 in rat spinal cord contributes to neuropathic pain by increasing neuron excitability after brachial plexus avulsion. J. Pain Res. 12, 1297–1310. doi: 10.2147/JPR.S195117

PubMed Abstract | CrossRef Full Text | Google Scholar

Nasoohi, S., Ismael, S., and Ishrat, T. (2018). Thioredoxin-Interacting protein (TXNIP) in cerebrovascular and neurodegenerative diseases: regulation and implication. Mol. Neurobiol. 55, 7900–7920. doi: 10.1007/s12035-018-0917-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Navratilova, E., Ji, G., Phelps, C., Qu, C., Hein, M., Yakhnitsa, V., et al. (2019). Kappa opioid signaling in the central nucleus of the amygdala promotes disinhibition and aversiveness of chronic neuropathic pain. Pain 160, 824–832. doi: 10.1097/j.pain.0000000000001458

PubMed Abstract | CrossRef Full Text | Google Scholar

Navratilova, E., Xie, J. Y., Meske, D., Qu, C., Morimura, K., Okun, A., et al. (2015). Endogenous opioid activity in the anterior cingulate cortex is required for relief of pain. J. Neurosci. 35, 7264–7271. doi: 10.1523/JNEUROSCI.3862-14.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Neugebauer, V. (2015). Amygdala pain mechanisms. Handb. Exp. Pharmacol. 227, 261–284. doi: 10.1007/978-3-662-46450-2_13

CrossRef Full Text | Google Scholar

Pan, Z., Shan, Q., Gu, P., Wang, X. M., Tai, L. W., Sun, M., et al. (2018). miRNA-23a/CXCR4 regulates neuropathic pain via directly targeting TXNIP/NLRP3 inflammasome axis. J. Neuroinflamm. 15:29. doi: 10.1186/s12974-018-1073-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedraza, C., Sánchez-López, J., Castilla-Ortega, E., Rosell-Valle, C., Zambrana-Infantes, E., García-Fernández, M., et al. (2014). Fear extinction and acute stress reactivity reveal a role of LPA (1) receptor in regulating emotional-like behaviors. Brain Struct. Funct. 219, 1659–1672. doi: 10.1007/s00429-013-0592-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Pevida, M., Lastra, A., Meana, Á, Hidalgo, A., Baamonde, A., Menéndez, L., et al. (2014). The chemokine CCL5 induces CCR1-mediated hyperalgesia in mice inoculated with NCTC 2472 tumoral cells. Neuroscience 259, 113–125. doi: 10.1016/j.neuroscience.2013.11.055

PubMed Abstract | CrossRef Full Text | Google Scholar

Phelps, E. A., and LeDoux, J. E. (2005). Contributions of the amygdala to emotion processing: from animal models to human behavior. Neuron 48, 175–187. doi: 10.1016/j.neuron.2005.09.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramos-Prats, A., Kölldorfer, J., Paolo, E., Zeidler, M., Schmid, G., Ferraguti, F., et al. (2019). An appraisal of the influence of the metabotropic glutamate 5 (mGlu5) receptor on sociability and anxiety. Front. Mol. Neurosci. 12:30. doi: 10.3389/fnmol.2019.00030

PubMed Abstract | CrossRef Full Text | Google Scholar

Rivera, P., Tovar, R., Ramírez-López, M. T., Navarro, J. A., Vargas, A., Suárez, J., et al. (2020). Sex-Specific anxiety and prefrontal cortex glutamatergic dysregulation are long-term consequences of pre-and postnatal exposure to hypercaloric diet in a rat model. Nutrients 12:1829. doi: 10.3390/nu12061829

PubMed Abstract | CrossRef Full Text | Google Scholar

Samuel, G. N., and Farsides, B. (2017). The UK’s 100,000 Genomes Project: manifesting policymakers’ expectations. New Genet. Soc. 36, 336–353. doi: 10.1080/14636778.2017.1370671

PubMed Abstract | CrossRef Full Text | Google Scholar

Sapio, M. R., Iadarola, M. J., LaPaglia, D. M., Lehky, T., Thurm, A. E., Danley, K. M., et al. (2019). Haploinsufficiency of the brain-derived neurotrophic factor gene is associated with reduced pain sensitivity. Pain 160, 1070–1081. doi: 10.1097/j.pain.0000000000001485

PubMed Abstract | CrossRef Full Text | Google Scholar

Savitz, J., Dantzer, R., Wurfel, B. E., Victor, T. A., Ford, B. N., Bodurka, J., et al. (2015). Neuroprotective kynurenine metabolite indices are abnormally reduced and positively associated with hippocampal and amygdalar volume in bipolar disorder. Psychoneuroendocrinology 52, 200–211. doi: 10.1016/j.psyneuen.2014.11.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Sellmeijer, J., Mathis, V., Hugel, S., Li, X. H., Song, Q., Chen, Q. Y., et al. (2018). Hyperactivity of anterior cingulate cortex areas 24a/24b drives chronic pain-induced anxiodepressive-like consequences. J. Neurosci. 38, 3102–3115. doi: 10.1523/JNEUROSCI.3195-17.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Seminowicz, D. A., Laferriere, A. L., Millecamps, M., Yu, J. S., Coderre, T. J., and Bushnell, M. C. (2009). MRI structural brain changes associated with sensory and emotional function in a ratmodel of long-term neuropathic pain. NeuroImage 47:10071014. doi: 10.1016/j.neuroimage.2009.05.068

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, Q., Colodner, K. J., Matousek, S. B., Merry, K., Hong, S., Kenison, J. E., et al. (2015). Complement C3-Deficient mice fail to display age-related hippocampal decline. J. Neurosci. 35, 13029–13042. doi: 10.1523/JNEUROSCI.1698-15.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Sierra-Fonseca, J. A., Parise, L. F., Flores-Ramirez, F. J., Robles, E. H., and Iñiguez, S. D. (2019). Dorsal hippocampus ERK2 signaling mediates anxiolytic-related behavior in male rats. Chronic stress (Thousand Oaks, Calif.) 3:2470547019897030. doi: 10.1177/2470547019897030

PubMed Abstract | CrossRef Full Text | Google Scholar

Simons, L. E., Moulton, E. A., Linnman, C., Carpino, E., Becerra, L., Borsook, D., et al. (2014). The human amygdala and pain: evidence from neuroimaging. Hum. Brain Mapp. 35, 527–538. doi: 10.1002/hbm.22199

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, M. L., Asada, N., and Malenka, R. C. (2021). Anterior cingulate inputs to nucleus accumbens control the social transfer of pain and analgesia. Science (New York, N.Y.) 371, 153–159. doi: 10.1126/science.abe3040

PubMed Abstract | CrossRef Full Text | Google Scholar

Stevens, A. M., Liu, L., Bertovich, D., Janjic, J. M., and Pollock, J. A. (2019). Differential expression of neuroinflammatory mRNAs in the rat sciatic nerve following chronic constriction injury and pain-relieving nanoemulsion NSAID delivery to infiltrating macrophages. Int. J. Mol. Sci. 20:5269. doi: 10.3390/ijms20215269

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, L., Gu, X., Pan, Z., Guo, X., Liu, J., Atianjoh, F. E., et al. (2019). Contribution of DNMT1 to neuropathic pain genesis partially through epigenetically repressing kcna2 in primary afferent neurons. J. Neurosci. 39, 6595–6607. doi: 10.1523/JNEUROSCI.0695-19.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Suzuki, T., Amata, M., Sakaue, G., Nishimura, S., Inoue, T., Shibata, M., et al. (2007). Experimental neuropathy in mice is associated with delayed behavioral changes related to anxiety and depression. Anesthesia Analgesia 104, 1570–1577. doi: 10.1213/01.ane.0000261514.19946.66

CrossRef Full Text | Google Scholar

Tang, Z., Gong, Z., and Sun, X. (2018). LncRNA DANCR involved osteolysis after total hip arthroplasty by regulating FOXO1 expression to inhibit osteoblast differentiation. J. Biomed. Sci. 25:4. doi: 10.1186/s12929-018-0406-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Temme, S. J., and Murphy, G. G. (2017). The L-type voltage-gated calcium channel CaV1.2 mediates fear extinction and modulates synaptic tone in the lateral amygdala. Learn. Memory (Cold Spring Harbor, N.Y.) 24, 580–588. doi: 10.1101/lm.045773.117

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomson, A. C., Kenis, G., Tielens, S., de Graaf, T. A., Schuhmann, T., Rutten, B. P. F., et al. (2020). Transcranial magnetic stimulation-induced plasticity mechanisms: TMS-Related gene expression and morphology changes in a human neuron-like cell model. Front. Mol. Neurosci. 13:528396. doi: 10.3389/fnmol.2020.528396

PubMed Abstract | CrossRef Full Text | Google Scholar

Toyoda, H., Zhao, M. G., Ulzhöfer, B., Wu, L. J., and Xu, H. (2009). Roles of the AMPA receptor subunit GluA1 but not GluA2 in synaptic potentiation and activation of ERK in the anterior cingulate cortex. Mol. Pain 5:46. doi: 10.1186/1744-8069-5-46

PubMed Abstract | CrossRef Full Text | Google Scholar

Tripp, A., Oh, H., Guilloux, J. P., Lewis, D. A., and Sibille, E. (2012). Brain-derived neurotrophic factor signaling and subgenual anterior cingulate cortex dysfunction in major depressive disorder. Am. J. Psychiatry 169, 1194–1202. doi: 10.1176/appi.ajp.2012.12020248

PubMed Abstract | CrossRef Full Text | Google Scholar

van Hecke, O., Austin, S. K., Khan, R. A., Smith, B. H., and Torrance, N. (2014). Neuropathic pain in the general population: a systematic review of epidemiological studies. Pain 155, 654–662. doi: 10.1016/j.pain.2013.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

von Hehn, C. A., Baron, R., and Woolf, C. J. (2012). Deconstructing the neuropathic pain phenotype to reveal neural mechanisms. Neuron 73, 638–652. doi: 10.1016/j.neuron.2012.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Zhai, W., Yu, Z., Sun, L., Li, H., Shen, H., et al. (2018). Neuroprotection EXerted by Netrin-1 and kinesin motor KIF1A in secondary brain injury following experimental intracerebral hemorrhage in rats. Front. Cell. Neurosci. 11:432. doi: 10.3389/fncel.2017.00432

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Zhao, M., Jia, P., Liu, F. F., Chen, K., Meng, F. Y., et al. (2020). The analgesic action of larixyl acetate, a potent TRPC6 inhibitor, in rat neuropathic pain model induced by spared nerve injury. J. Neuroinflamm. 17:118.

Google Scholar

Wen, J., Yang, Y., Wu, S., Wei, G., Jia, S., Hannaford, S., et al. (2020). Long noncoding RNA H19 in the injured dorsal root ganglion contributes to peripheral nerve injury-induced pain hypersensitivity. Transl. Perioperat. Pain Med. 7, 176–184.

Google Scholar

Weng, H. R., Gao, M., and Maixner, D. W. (2014). Glycogen synthase kinase 3 beta regulates glial glutamate transporter protein expression in the spinal dorsal horn in rats with neuropathic pain. Exp. Neurol. 252, 18–27. doi: 10.1016/j.expneurol.2013.11.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, J., Wang, C., and Ding, H. (2020). LncRNA MALAT1 promotes neuropathic pain progression through the miR-154-5p/AQP9 axis in CCI rat models. Mol. Med. Rep. 21, 291–303.

Google Scholar

Wu, S., Marie Lutz, B., Miao, X., Liang, L., Mo, K., Chang, Y. J., et al. (2016). Dorsal root ganglion transcriptome analysis following peripheral nerve injury in mice. Mol. Pain 12:1744806916629048. doi: 10.1177/1744806916629048

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, X. B., He, L. N., Jiang, B. C., Wang, X., Lu, Y., Gao, Y. J., et al. (2019). Increased CXCL13 and CXCR5 in anterior cingulate cortex contributes to neuropathic pain-related conditioned place aversion. Neurosci. Bull. 35, 613–623. doi: 10.1007/s12264-019-00377-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, M., Chen, W., Zhang, Y., Yang, R., Wang, Y., Yuan, H., et al. (2018). EphrinB/EphB signaling contributes to spinal nociceptive processing via calpain-1 and caspase-3. Mol. Med. Rep. 18, 268–278. doi: 10.3892/mmr.2018.8996

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, E. E., Bryant, C. D., Lee, S. E., Peng, X., Cook, B., Nair, H. K., et al. (2016). Systems genetic and pharmacological analysis identifies candidate genes underlying mechanosensation in the von Frey test. Genes Brain Behav. 15, 604–615. doi: 10.1111/gbb.12302

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, H., Zhang, P., Chen, Y. R., Wang, Y. J., Lin, X. Y., Li, X. Y., et al. (2019). Temporal changes of spinal transcriptomic profiles in mice with spinal nerve ligation. Front. Neurosci. 13:1357. doi: 10.3389/fnins.2019.01357

PubMed Abstract | CrossRef Full Text | Google Scholar

Zannas, A. S., Jia, M., Hafner, K., Baumert, J., Wiechmann, T., Pape, J. C., et al. (2019). Epigenetic upregulation of FKBP5 by aging and stress contributes to NF-κB-driven inflammation and cardiovascular risk. Proc. Natl. Acad. Sci. U S A. 116, 11370–11379. doi: 10.1073/pnas.1816847116

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, X., Tang, Z., Zhang, H., Atianjoh, F. E., Zhao, J. Y., Liang, L., et al. (2013). A long noncoding RNA contributes to neuropathic pain by silencing Kcna2 in primary afferent neurons. Nat. Neurosci. 16, 1024–1031. doi: 10.1038/nn.3438

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, Y., Xu, H., Zhan, L., Zhou, X., Chen, X., Gao, Z., et al. (2015). Activation of peripheral KCNQ channels relieves gout pain. Pain 156, 1025–1035. doi: 10.1097/j.pain.0000000000000122

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, J., Xiong, Q., Chen, H., Yang, C., and Fan, Y. (2017). Identification of the spinal expression profile of non-coding rnas involved in neuropathic pain following spared nerve injury by sequence analysis. Front. Mol. Neurosci. 10:91. doi: 10.3389/fnmol.2017.00091

PubMed Abstract | CrossRef Full Text | Google Scholar