- Open Access
Global gene expression and chromatin accessibility of the peripheral nervous system in animal models of persistent pain
Journal of Neuroinflammation volume 18, Article number: 185 (2021)
Efforts to understand genetic variability involved in an individual’s susceptibility to chronic pain support a role for upstream regulation by epigenetic mechanisms.
To examine the transcriptomic and epigenetic basis of chronic pain that resides in the peripheral nervous system, we used RNA-seq and ATAC-seq of the rat dorsal root ganglion (DRG) to identify novel molecular pathways associated with pain hypersensitivity in two well-studied persistent pain models induced by chronic constriction injury (CCI) of the sciatic nerve and intra-plantar injection of complete Freund’s adjuvant (CFA) in rats.
Our RNA-seq studies identify a variety of biological process related to synapse organization, membrane potential, transmembrane transport, and ion binding. Interestingly, genes that encode transcriptional regulators were disproportionately downregulated in both models. Our ATAC-seq data provide a comprehensive map of chromatin accessibility changes in the DRG. A total of 1123 regions showed changes in chromatin accessibility in one or both models when compared to the naïve and 31 shared differentially accessible regions (DAR)s. Functional annotation of the DARs identified disparate molecular functions enriched for each pain model which suggests that chromatin structure may be altered differently following sciatic nerve injury and hind paw inflammation. Motif analysis identified 17 DNA sequences known to bind transcription factors in the CCI DARs and 33 in the CFA DARs. Two motifs were significantly enriched in both models.
Our improved understanding of the changes in chromatin accessibility that occur in chronic pain states may identify regulatory genomic elements that play essential roles in modulating gene expression in the DRG.
Despite intense research efforts to develop novel analgesic classes, few novel molecular targets have been successfully translated into effective pain therapy . While non-steroidal anti-inflammatory drugs and opioids continue to be the most effective drugs commonly prescribed for the treatment of persistent pain, they have been associated with significant adverse effects. Therefore, increased attention is now focused on genetic and epigenetic mechanisms as an avenue to identify new druggable targets [2, 3].
Current evidence supports the association between changes in gene expression and the transition from acute to chronic pain states in a number of preclinical and clinical models [4, 5]. However, as these models are developed using nerve injury, administration of chemical agents, or evoking a significant inflammatory response, difficulties arise when disentangling gene expression profiles due to the effects of pain behaviors versus the initiating insults. Another factor that complicates the interpretation of transcriptional changes across different pain states is that existing microarray datasets from various preclinical models of persistent pain suffer from poor accuracy for genes expressed in low abundance and at lower coverage. Since these limitations can largely be overcome with newer RNA-seq approaches, the use of RNA-seq to identify shared changes in gene expression in disparate preclinical models is a promising strategy to identify pain-specific gene expression changes.
Chromatin structure is a well-known regulator of gene transcription across eukaryotic organisms [6, 7]. However, chromatin structure in the peripheral nervous system (PNS) remains poorly understood. Intriguingly, recent work has indicated the involvement of epigenomic changes in the regulation of gene expression in the dorsal root ganglion (DRG), which contains the soma of peripheral sensory neurons, in a preclinical model of neuropathic pain . Long-term changes to gene expression patterns likely depend upon modifications to chromatin structure in post-mitotic neurons as well . Thus, changes in chromatin structure at DNA regulatory regions in DRG neurons likely foster long-term changes in membrane potential and excitability, and thus, promote maintenance of persistent pain states. Histone acetylation and DNA methylation have been identified in persistent pain models, and our improved understanding of these and other epigenetic mechanisms through which aberrant gene expression could occur in the PNS unveil early molecular events that underlie the maintenance of chronic pain states and inform novel analgesic treatments.
Here, we used both RNA sequencing and the Assay for Transposase-Accessible Chromatin by sequencing (ATAC-seq) to identify patterns of genome-wide chromatin accessibility and gene transcription in the lumbar DRGs employing two widely used rodent models of neuropathic and inflammatory pain. The differences in chromatin accessibility observed between naïve and pain states in each model allowed us to identify regulatory DNA sequences and putative transcription factors that may drive the changes in gene expression associated with hyperalgesic states. Our integrative approach allowed us a greater understanding of the transcriptional and epigenetic basis of chronic pain in the DRG and identified novel biological processes and regulatory intermediates that may lead to the long-term transcriptional changes associated with persistent pain in the PNS.
Adult male Sprague-Dawley rats (12 weeks old; Harlan Bioproducts for Science, Indianapolis, IN) were allowed to acclimate for a minimum of 48 h prior to any experimental procedures. Animals had access ad libitum to food and water.
Establishment of pain models
CCI of sciatic nerve
Chronic constriction injury (CCI) surgery to the sciatic nerve was performed as previously described . Under 2–3% isoflurane, a small incision was made at the level of the mid-thigh. The sciatic nerve was exposed by blunt dissection through the biceps femoris. The nerve trunk proximal to the distal branching point was loosely ligated with four 4-O silk sutures placed approximately 0.5 mm apart until the epineurium was slightly compressed and minor twitching of the relevant muscles was observed. The muscle layer was closed with 4-O silk suture and the wound closed with metal clips. Pain-related behaviors are highly variable during the first 7 days after CCI, become more consistent, and peak at day 14 before returning to baseline levels 4–6 weeks after surgery .
Intraplantar injection of CFA
Complete Freund’s adjuvant (CFA) (Sigma-Aldrich, St. Louis, MO) was diluted 1:1 with sterile 0.9% saline to produce a 0.5 mg/ml emulsion. Under 2–3% isoflurane, the plantar surface of each hind paw was cleaned and injected with 100 μl of the 50% CFA emulsion using a 27-guage hypodermic needle. Following CFA injection, mechanical hypersensitivity rapidly developed in the ipsilateral hind paw within 2–6 h and persists for 1–2 weeks . Rats were given 48 h from the time of injection for the persistent hypersensitivity to reach its peak level and stabilize. All animals underwent behavior testing to verify presence of hind paw mechanical hypersensitivity using von Frey monofilaments . Any animal that did not show a 50% reduction of paw withdrawal threshold from the pre-injury baseline (i.e., the criteria for the presence of mechanical hypersensitivity) at 14 days for sciatic nerve CCI or 48 h after CFA injection were not used. Absence of mechanical hypersensitivity was verified in all naïve animals.
RNA-sequencing and data processing
Total RNA was extracted from pooled ipsilateral lumbar (L4–6) DRGs from one rat using the QIAGEN RNeasy Mini Prep Kit (QIAGEN, Valencia, CA) with on-column DNase digestion (QIAGEN, Valencia, CA) according to manufacturer’s instructions. RNA concentration was measured using the Nanodrop ND-2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA) and RNA integrity was assessed using RNA Nano Eukaryote chips in an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA).
RNA-seq library construction and sequencing
Five hundred nanograms of total RNA per sample was used to construct sequencing libraries (n = 1 rat/sample). Strand-specific RNA libraries were prepared using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA) after poly(A) selection by the NEBNext poly(A) mRNA Isolation Module (New England Biolabs, Ipswich, MA) according to manufacturer’s instructions. Samples were barcoded using the recommended NEBNext Multiplex Oligos (New England Biolabs). Size range and quality of libraries were verified on the Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). RNA-seq libraries were quantified by qPCR using the KAPA library quantification kit (KAPA Biosystems, Wilmington, MA). Each library was normalized to 2 nM and pooled in equimolar concentrations. Paired-end sequencing was performed in a single lane on an Illumina HiSeq2500 (Illumina, San Diego, CA). Two independent biological replicates of 1 rat per group were run for a total of 6 libraries. Libraries were sequenced to an average depth of 33.9 million reads per sample.
RNA-seq data analysis
Sequencing reads were aligned to annotated RefSeq genes in the rat reference genome (rn6) using HISAT2 , filtered to remove ribosomal RNA, and visualized using the Integrative Genomics Viewer . A gene count matrix that contained raw transcript counts for each annotated gene were generated using the featureCounts function of the Subread package in R  against the Ensemble rn6 transcriptome. This count matrix was then filtered for low count genes so that only those genes with > 0 reads across all samples were retained. We relied on the automatic and independent filtering used by DESeq2 to determine the most appropriate threshold for removing genes with low counts .
To identify genes that were differentially regulated following nerve injury, raw transcript counts were normalized, log2 transformed, and analyzed using the default procedures in DESeq2 . All downstream analyses on RNA-seq data were performed on data obtained from DESeq2. Adjusted p values were corrected using the Benjamini-Hochberg method. An adjusted p value < 0.05 and an absolute log2 fold change > 0.5 were used to define differentially expressed transcripts between naive and each of the chronic pain models. DESeq2 adjusts for multiple testing by implementing the procedures of Benjamini and Hochberg . Genes with differential expression between groups were then included in gene ontology (GO) and pathway analysis to infer their functional roles and relationships. GO analysis for enriched GO biological processes in each set of differentially enriched genes (DEGs) identified by DESeq2 was performed using Metascape . We previously validated our RNA-seq data using qPCR in biological replicates .
ATAC-seq library preparation
Immediately following dissection, the ipsilateral lumbar (L4–L6) DRGs from one rat were transferred directly to cold lysis buffer (0.32 M sucrose, 5 mM calcium chloride, 3 mM magnesium acetate, 10 mM Tris-hydrochloride, pH 8.0, 0.1% Triton X-100, 1 mM dithiothreitol, 5 mM sodium butyrate, 1 mM phenylmethylsulfonyl fluoride). Nuclei were isolated through dounce homogenization of the tissue in lysis buffer followed by sucrose gradient ultracentrifugation (1.8 M sucrose, 3 mM magnesium acetate, 1 mM dithiothreitol, 10 mM Tris-hydrochloride, pH 8.0, 5 mM sodium butyrate, 1 mM phenylmethylsulfonyl fluoride) at 139,800×g at 4 °C for 1 h to remove mitochondrial DNA. The nuclei were resuspended in 1× phosphate-buffered saline and counted 3 times using a Neubauer chamber. Tagmentation by Tn5 was performed using reagents from the Nextera DNA Sample Preparation Kit (FC-121-1030, Illumina; San Diego, CA) as previously described . Each 50 μl reaction contained 50,000 nuclei, 25 μl 2× Tagmentation Buffer, and 2.5 μl TDE1 enzyme and incubated at 37 °C for 30 min. Tagmented DNA was immediately purified using the Clean and Concentrate-5 Kit (Zymo, Irvine, CA) and eluted in 10 μl elution buffer. Tagmented DNA fragments were amplified using Nextera Index adapters, polymerase chain reaction (PCR) primer cocktail, NPM PCR master mix, and 10 cycles of PCR. Each library was purified using Agencourt AMPure XP beads (Beckam Coulter, Atlanta, GA). The fragment distribution of each library was assessed the using the High Sensitivity DNA Kit on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). Libraries were quantified prior to sequencing using the Qubit DNA High Sensitivity kit (ThermoScientific, Waltham, MA) and normalized to 2 nM and pooled in equimolar concentrations. Libraries were sequenced using paired end, dual-index sequencing on an Illumina HiSeq 2500 (Illumina, San Diego, CA) which produced 50 base pair reads. Six independent biological replicates of 1 rat/sample were prepared for each of the CCI and CFA groups and 3 independent biological replicates were prepared for the Naïve group for a total of 15 libraries.
ATAC-seq data processing
The paired-end reads were trimmed using Trimmomatic  to remove adaptors. The trimmed reads were then aligned to rat genome rn6 using Bowtie2  with the following parameters -X2000 --no-mixed --no-discordant. Reads with mapping quality score less than 10 were removed using SAMtools  and duplicated reads were removed using the MarkDuplicates function in Picard (http://broadinstitute.github.io/picard/). Aligned reads were shifted 4 nucleotides upstream for the 5′ end and 5 nucleotides downstream for the 3′ end to remove potential artifacts of Tn5 transposase binding. Tn5 transposase insertion sites were identified by trimming each read to the 5′ end. Bedtools slop was used to extend the insertion site by 75bp upstream and downstream . Reads for each sample were downsampled to 49 million insertion sites to account for differences in sequencing depth. ATAC-seq peaks for each sample were called on the down-sampled bed files using callpeak function of Model-based Analysis of ChIP-seq (MACS2) and the parameters –nomodel –extsize 150 -B –keepdup all –call-summits. Bigwig files were also generated from the down-sampled bed files for visualization in Integrative Genomics Viewer . To improve our confidence in the selection of accessible regions in the dataset, accessible regions were only considered for downstream analysis if the region was called by MACS2 in at least 50% of all samples from that group. A consensus peakset was then determined by the overlap of these regions.
The number of reads that aligned to each peak were counted and differential accessibility at each peak between the CFA and Naïve group and the CCI and Naïve group were determined using the limma package (version 3.38.3) in R (version 3.5.1) . A p value < 0.001 was used to define differentially accessible peaks between naïve and each of the chronic pain models. The genomic feature and the nearest annotated gene were determined using the annotatePeaks.pl function with the Hypergeometric Optimization of Motif Enrichment (HOMER) algorithm (version 4.11.1) . De novo sequence motif discovery was to identify over representation of transcription factor binding sites within differentially accessible regions (DAR)s using the findMotifsGenome.pl function within HOMER .
HEK293 cells were purchased from the American Type Culture Collection (CRL-1573™; Rockville, MD). Cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM, Sigma, St. Louis, MO) supplemented with 10% fetal calf serum and penicillin-streptomycin and incubated at 37 °C in a humidified environment containing 5% CO2. Cells with low passage numbers (i.e., < 20) were used for all experiments.
Luciferase reporter constructs were constructed by cloning a candidate enhancer region into the pGL3 promoter vector (Promega; Madison, WI). Each region was inserted using standard restriction enzyme-based cloning techniques. The regions were obtained by PCR of rat genomic DNA. The 5′ end of the primers were modified to contain BglII (Forward primer) and MluI (reverse primer) restriction sites (Supplemental Table 1). PCR was performed using the Pfu Turbo polymerase (Agilent Technologies; Palo Alto, CA) and touchdown thermocycling. The PCR products were digested and ligated into the BglII (AGATCT) and MluI (ACGCGT) restriction enzyme sites of the pGL3-Promoter luciferase vector (Promega; Madison, WI). The ligated products were transformed into chemically competent DH5α cells using ampicillin (100 μg/ml) to select for the recombinant plasmid-positive colonies. All constructs were verified by restriction enzyme digest and Sanger sequencing.
Transfection and luciferase assays
HEK293 cells were seeded at 20,000 cells/well in 48-well plates in 250 μl of complete media and grown to 60–70% confluence. Cells were then transfected with each reporter construct (450 ng) and 50 ng pGL4.74 Renilla luciferase expression vector (Promega; Madison, WI) using ViaFect Transfection Reagent (Promega; Madison. WI) in 25 μl Opti-MEM (ThermoScientific, Waltham, MA) with a 4:1 ratio in 250 μl complete medium. The transfection efficiency of HEK293 cells was evaluated by transfecting cells with EGFP-N1 (Clontech; Mountain View, CA) in parallel reactions. 48 h post transfection, Firefly and Renilla luciferase activities were measured using the Dual-Glo Luciferase reporter assay system (Promega). Firefly luciferase activities were normalized to the Renilla luciferase activity and expressed as the relative fold difference of the empty pGL3 promoter vector. Each luciferase construct was tested in quadruplicate.
Differential gene expression changes after injury in neuropathic pain and inflammatory pain models
To determine how gene expression is altered in the lumbar DRG following the establishment of two widely used rat models of persistent pain, we compared RNA-seq data obtained 14 days following nerve injury (i.e., CCI model) to naïve rats and 48 h following hind paw inflammation (i.e., CFA) to naïve rats (Fig. 1A). Principal component analysis (PCA) shows clear segregation of the transcriptomes from the CCI and CFA pain models and naïve rats (Fig. 1B). The first two principal components accounted for a total of 91%.
Compared to naïve rats, we identified 2620 (17.8%) DEGs in the DRG following CCI. Of these 2620 DEGs, 1111 (42.4%) genes were upregulated following CCI as compared to naïve and 1509 (57.6%) genes were downregulated (Fig. 1C). Gene ontology analysis of the 1111 upregulated genes indicated that neuronal-activity related biological processes including synapse organization, regulation of ion transport, modulation of chemical synaptic transmission, and sensory perception of pain were among the biological processes that were statistically enriched (Fig. 1D). Following CFA, we identified 2773 (18.8%) genes that were differentially expressed in the DRG when compared to naïve rats. Of these 2773 DEGs, 1119 (40.4%) genes were upregulated and 1654 (59.6%) genes were downregulated (Fig. 1C). As expected, gene ontology analysis of the 1119 upregulated genes revealed neuronal- and pain-related biological processes including synaptic vesicle cycle, regulation of membrane potential, and transmission of nerve impulse were among the biological processes that were statistically enriched (Fig. 1D). CCI- and CFA-specific changes in gene transcription that are important in pre-synaptic activity are shown in Fig. 1E.
Shared differential expression of genes in CCI and CFA models
A total of 1960 genes (13.3% of the 14,742 genes examined) were differentially expressed in both pain models with 752 genes (38.4%) upregulated and 1198 genes (61.1%) downregulated (Fig. 2A; Supplemental Table 2). GO analysis of the 752 upregulated genes showed significant enrichment among a variety of biological processes related to synapse organization and membrane potential (Fig. 2B top). GO analysis of the 1198 downregulated transcripts show significant enrichment among biological processes involved in transmembrane transport and ion binding (Fig. 2B bottom).
The numbers of upregulated genes that encoded GPCRs (Fig. 2C, D) and ion channels (Fig. 2C, E) were markedly higher than those of downregulated genes. Among the 28 GPCR genes upregulated following CCI and CFA, several are known to be involved in abnormal generation of action potentials and altered pain thresholds (e.g., Gal1r, Gabbr2, Grm7, Gpr158, Chrm2, Chrm3, Cnr1, Oprl1, Oprm1, Mrgpre) as well as clinical pain conditions such as neuralgia (e.g., Htr2a, Htr2b, Ptger4, Mrgprx1) (Fig. 2D) [25,26,27,28,29,30,31,32,33,34,35]. Conversely, of the genes identified as transcriptional regulators, the number downregulated was 10-fold higher than the number upregulated (Fig. 2C, F; Supplemental Table 3).
Of note, 10 genes (i.e., Myot, Ca3, Tnnc2, Ankrd23, Eno3, Casq1, Tpm2, Pygm, RT1-Da, Des) were upregulated after CCI and downregulated after CFA. The majority of these genes are involved in the function and maintenance of skeletal muscle and would be expected to be upregulated following our surgical dissection of the biceps femoris, and injury to axons that innervate muscle tissues in the CCI model [36, 37].
ATAC-seq provides a high-resolution chromatin accessibility profile of the rat DRG
Given our findings of gene expression changes in two persistent pain models with different etiologies, we hypothesized that these transcriptional changes were a result of dynamic chromatin occupancy in DRG cells which change the accessibility to cis-regulatory elements for transcriptional machinery. Therefore, we performed ATAC-seq of rat DRGs to determine the genome-wide dynamics of chromatin accessibility of the DRG in naïve rats and in the two pain models. Ninety to 94% of paired-end reads mapped to the rat reference genome and provided an average mapping depth of 17.8 million reads per sample. This represents an average of 35.6 million insertion sites per sample.
We identified a total of 97,485 peaks across all groups (i.e., Naïve, CCI, CFA). Between 53 and 69% of the peaks called in each group were found in only one sample in each group (Supplemental Figure 1). Therefore, to increase our confidence that we were identifying true regions of open accessibility, we retained 56,810 accessible regions that were identified at least half of the samples in each group for use in all downstream analysis. With 33,628 regions (57.4% of consensus regions), we found that the CCI group had the smallest number of accessible sites. The CFA group had 46,238 (81.4%) accessible regions and the naïve group had 45,399 (79.9%) accessible regions.
To determine the reproducibility among replicates, we performed principal component analysis and calculated Pearson correlation coefficients. We found that biological replicates were highly correlated with other samples from the same treatment group and less correlated with samples from either of the other 2 treatment groups. By PCA, samples show separation by group (Fig. 3A). As expected, almost half (48.8%) of these regions were present in all three groups (Fig. 3B). The distance between chromatin accessible regions and the nearest gene transcription start site suggests that these regions are concentrated in cis-regulatory regions (Fig. 3C). Indeed, 54.7%, 28%, and 11.6% of all consensus accessible regions were located in intergenic, intronic, or promoter regions, respectively.
Changes in chromatin accessibility in the DRG after the establishment of persistent pain
We evaluated each of the 56,810 consensus regions for changes in chromatin accessibility in each of the two pain models (p value < 0.001). A total of 1123 (2.0%) of the 56,810 consensus regions showed changes in accessibility in one or both pain models compared to naïve rats. We found 517 DARs in the DRG following CCI compared to naïve with 426 regions having increased accessibility and 91 regions having decreased accessibility (Fig. 4A, B). When comparing the CFA model to naïve, we found 637 DARs with 321 regions with increased accessibility and 316 regions with decreased accessibility following CFA (Fig. 4D, E). Sixty-two percent of all gains and losses in chromatin accessibility after CCI or CFA injection occurred in intergenic regions while 21.3% and 9.8% were located in introns and gene promoters, respectively (Fig. 4G). These changes in genomic distribution of chromatin accessibility after sciatic nerve injury (CCI) or hind paw inflammation (CFA) may facilitate differential gene transcription through chromatin-level regulation at cis-regulatory regions.
Functional annotation of the DARs after CCI shows enrichment for molecular functions involved in phosphatidylinositol 3-kinase regulator activity, and neurotransmitter receptor activity involved in the regulation of postsynaptic membrane potential and enhancer binding (Fig. 4C). These functions converge on mechanisms with the potential to alter neuronal excitability and chromatin structure in DRG cells to produce pain hypersensitivity. Following CFA injection, enrichment in potassium ion transmembrane transport, receptor serine/threonine kinase binding, and structural constituent of synapses are among the molecular functions found in DARs (Fig. 4F). The different molecular functions identified through GO analysis suggest that chromatin structure may be altered to different effect following nerve injury and hind paw inflammation.
A total of 31 DARs were shared following CCI and CFA injection (Table 1). Ten DARs showed decreased accessibility after CCI and CFA compared to naïve. Of these 10 DARs, 2 were associated with upregulation of the nearest annotated gene (i.e., Grik4, Agtpbp1) and 2 were associated with downregulation of their nearest gene (i.e., Wdr60, Dlg3). A total of 21 DARs showed increased accessibility after CCI and CFA compared with Naïve. Of these 21 DARs, 3 were associated with upregulation of the nearest gene (i.e., Rrmb, Arap2, Pcdh9) and 2 were associated with downregulation of the nearest gene (i.e., Kif13b, Lpar1). The ability of these regions to modulate gene transcription were validated by luciferase assay (Fig. 4H).
To determine transcription factors that may have their binding sites within with regions of chromatin accessibility in each of the pain models, we performed motif analyses using HOMER on all of the DARs identified in each pain model (Fig. 5A and Supplemental Tables 4 and 5). We found a total of 17 DNA sequences known to bind transcription factors in the CCI DARs and 33 in the CFA DARs. Binding motifs for Foxo1 and Hlf were significantly enriched following both CCI and CFA (Fig. 5B). CFA was associated with the induction of a wider remodeling of transcription factor binding sites in the DRG than CCI. These enriched sites found after CFA clustered into several important families (e.g., high-mobility group (HMG), basic leucine zipper (bZIP), transcriptional enhancer factor (TEF)). CCI was associated with enriched sites clustered into the ETS-domain (ETS) and zinc finger (Zf) families. The changes in the availability of potential TF binding sites support the long-lasting effects of our pain models on chromatin structure.
In this study, we used both RNA-seq and ATAC-seq on DRG tissues to improve our understanding of the transcriptomic and epigenetic mechanisms in the PNS that may underlie persistent pain. We also aimed to identify novel molecular pathways involved in the development of pain hypersensitivity in two well-studied rodent models of persistent pain with different etiologies. Early research to study gene expression changes in preclinical pain models primarily relied on microarrays to study gene expression changes in the DRG . The decreasing cost of next-generation sequencing has resulted in adoption of RNA-seq as the current standard which has a greater dynamic range for gene expression detection, the ability to measure a larger number of gene transcripts, and can detect differences in sequence and isoforms. Studies that have used RNA-seq to look at DRG changes have primarily examined peripheral nerve injury models [38,39,40,41] with few examining non-nerve injury models (e.g., diabetic neuropathy , ultraviolet-induced inflammation ) which makes it difficult to identify gene expression changes that are specific to pain process versus nerve injury or specific disease-related processes.
In this study, we performed RNA-seq using two widely used persistent pain models which are induced by sciatic nerve CCI and hind paw inflammation to identify pain-related gene expression changes in the DRG. Consistent with prior findings, both nerve injury and inflammatory pain were found to be associated with transcriptional changes of genes involved in cell signaling (i.e., GPCR function, ion channel expression, synaptic transmission). Further, our findings were consistent with previous transcriptome-wide screens that support the upregulation of Reg3b, Vgf, Ccl2, P2rx3, Crh, Scn11a, Drd2, Npy2r, Cacna2d1, and Neto1 in the DRG in various pain models [4, 40, 43, 44]. In addition, our findings are consistent with prior work that found genes downregulated (e.g., Rlbp1, Gja1, Lgr5, Lpar1, Ttyh1) in DRG neurons after CCI . Here, we confirm that these genes are also downregulated in an inflammatory pain model which suggests that these genes have broader roles in pain pathways outside of nerve injury induced neuropathic pain.
Our studies identify several novel genes with previously unknown functions in the development and maintenance of pain. For example, Plxna2 was upregulated in both CCI and CFA models. Plxna2 encodes the plexin-A2 receptor known to be expressed in hippocampal and cortical neurons [45, 46]. Upon binding by its ligands, semaphorin-3A or -6A, plexin-A2 triggers an intracellular signaling cascade which mediates axon repulsion and cell migration during nervous system development [47, 48]. Further research is needed to determine how Plxna2 may participate in pronociceptive mechanisms.
Another salient finding in the current study is that genes involved in epigenetic and transcriptional regulation were largely downregulated following both CCI or CFA injection. Consistent with prior evidence, we found that Dnmt3a, Dnmt3b, Sirt2, Brd3, and Ehmt2 (which encodes G9a) were among the genes downregulated in both pain models, although the magnitude of the log2 fold change did not meet our criteria to be significant in the present study [25, 49, 50]. Mounting evidence supports the roles for epigenetic mechanisms and our findings that a large proportion of transcriptional regulators are downregulated supports this idea. Genes that decrease the regulation of existing transcriptional programs promote the transcriptional changes which are now well established in preclinical and clinical models of persistent pain. These changes in gene transcription and transcriptional regulation may facilitate neuronal hyperexcitability induced remodeling of chromatin structure and neuronal responsivity of cells in the DRG. Our findings are consistent with prior studies which provides indirect evidence of decreased transcriptional control in persistent pain states. Histone deacetyltransferase inhibitors and histone acetyltransferases show analgesic effects in various pain models through non-specific changes in chromatin structure.
One of the limitations in the study of epigenetic mechanisms in the PNS is the small number of cells that comprise each DRG. Traditional chromatin accessibility assays require millions of cells. However, ATAC-seq can assess native chromatin accessibility with much less starting material, and therefore, can detect subtle changes in chromatin accessibility both in homogenous cell lines and in heterogenous mixtures of cell types .
Our study is the first to provide a comprehensive map of changes in chromatin accessibility in the PNS using both neuropathic pain and inflammatory pain models. Chromatin accessibility is necessary for transcription factor binding to cis-regulatory regions and subsequent changes in gene expression [7, 51]. The DRG is a collection of several types of primary sensory neuronal cell bodies and satellite glia cells which acts as the initial point of modulation of action potentials from potentially noxious stimuli. Gene expression changes in the DRG control noxious input to second order neurons in the spinal cord. Therefore, understanding changes in chromatin accessibility that occur in chronic pain states would help to identify regulatory loci in the DRG that play a central role in changing gene expression there. We identified 1123 regions that showed dynamic chromatin accessibility after CCI and/or CFA. Within these regions, we found 48 known DNA binding motifs for transcription factors. Importantly, we identified overrepresentation for the DNA binding motifs of two transcription factors, FOXO1 and HLF, in DARs after both CCI and CFA injection. For example, we found that in both the CCI and CFA models, decreased chromatin accessibility at a 358bp region ~ 130 Kb downstream of Grik4 is associated with increased Grik4 gene expression (Fig. 5C). Grik4 encodes a subunit of a glutamatergic receptor that contributes to excitatory postsynaptic currents and is expressed in the central terminals of nociceptive neurons that synapse in laminae I–III . Over expression of Grik4 in the mouse brain is associated with increased amplitude, greater frequency, and quicker decay of spontaneous EPSCs in CA3 cells [53, 54]. While the effects of increased Grik4 expression in DRG neurons have not been investigated, these alterations in synaptic transmission are consistent with the increased efficacy of synaptic transmission of nociceptive input of central sensitization. The 358 bp DAR we identified upstream of Grik4 in the rat maps to chr11:120511578-120511909 on hg19 and is located within a previously identified regulatory region for Grik4 (i.e., GeneHancer #GH11J120511). Our findings suggest that this site is a transcriptional repressor binding site and that decreased accessibility inhibits this repression to increase Grik4 gene expression. Indeed, we found a potential Foxo1 binding site within this DAR (Fig. 5C). Foxo1 is well established as a transcriptional repressor in neural tissues. Further research is needed to evaluate this region for regulatory potential and determine if Foxo1 binds to this region in the rat DRG and can modulate pain behaviors.
DRG neurons differ in soma size, neurochemical properties, and functional specificity. Different subtypes of DRG neurons may show variable genetic and epigenetic changes after injury. Our bulk ATAC-seq provides an average of chromatin accessibility at a specific genomic locus. Change in chromatin accessibility may be due to an increase in the number of cell types in which a specific region is accessible or an increase in accessibility in the existing population of cells . As low-cell number methodologies mature, future research at the single-cell level may identify the epigenetic mechanisms responsible for the changes in chromatin structure in specific neuronal subtypes. Data from of chromatin confirmation assays in combination with multimodal epigenome and transcriptome studies will further our understanding of how altered transcriptional regulation may promote chronic pain.
In our study, we identified genes that are similarly differentially expressed in two disparate pain models as well as cis-regulatory regions with changes in chromatin accessibility that may drive these changes in gene expression and then validated the regulatory potential of these putative regulatory regions were indeed capable of altering gene expression. Our improved understanding of these regulatory genomic elements that may play pain-specific roles in modulating gene expression in the DRG provides new molecular insights into important changes the PNS that influence sensitivity and susceptibility to chronic pain.
Availability of data and materials
Raw and processed sequencing data for all ATAC-seq data and the RNA-seq data for the CFA samples have been deposited in the NCBI GEO database under accession #GSE167369. The RNA-seq data files for the naïve and CCI groups were previously published  and are available under accession #GEO100122.
Assay for Transposase-Accessible Chromatin by sequencing
Chronic constriction injury
Complete Freund’s adjuvant
Differentially accessible region
Differentially expressed gene
Dorsal root ganglion
Hypergeometric Optimization of Motif Enrichment
Polymerase chain reaction
Peripheral nervous system
Kissin I. The development of new analgesics over the past 50 years: a lack of real breakthrough drugs. Anesth Analg. 2010;110(3):780–9. https://doi.org/10.1213/ANE.0b013e3181cde882.
Louwies T, Ligon CO, Johnson AC, Greenwood-Van MB. Targeting epigenetic mechanisms for chronic visceral pain: a valid approach for the development of novel therapeutics. Neurogastroenterol Motil. 2019;31(3):e13500. https://doi.org/10.1111/nmo.13500.
Niederberger E, Resch E, Parnham MJ, Geisslinger G. Drugging the pain epigenome. Nat Rev Neurol. 2017;13(7):434–47. https://doi.org/10.1038/nrneurol.2017.68.
LaCroix-Fralish ML, Austin JS, Zheng FY, Levitin DJ, Mogil JS. Patterns of pain: meta-analysis of microarray studies of pain. Pain. 2011;152(8):1888–98. https://doi.org/10.1016/j.pain.2011.04.014.
Stephens KE, Zhou W, Ji Z, Chen Z, He S, Ji H, et al. Sex differences in gene regulation in the dorsal root ganglion after nerve injury. BMC Genomics. 2019;20(1):147. https://doi.org/10.1186/s12864-019-5512-9.
Li B, Carey M, Workman JL. The role of chromatin during transcription. Cell. 2007;128(4):707–19. https://doi.org/10.1016/j.cell.2007.01.015.
Zaret KS, Carroll JS. Pioneer transcription factors: establishing competence for gene expression. Genes Dev. 2011;25(21):2227–41. https://doi.org/10.1101/gad.176826.111.
Palmisano I, Danzi MC, Hutson TH, Zhou L, McLachlan E, Serger E, et al. Epigenomic signatures underpin the axonal regenerative ability of dorsal root ganglia sensory neurons. Nat Neurosci. 2019;22(11):1913–24. https://doi.org/10.1038/s41593-019-0490-4.
Sweatt JD. The emerging field of neuroepigenetics. Neuron. 2013;80(3):624–32. https://doi.org/10.1016/j.neuron.2013.10.023.
Bennett GJ, Xie YK. A peripheral mononeuropathy in rat that produces disorders of pain sensation like those seen in man. Pain. 1988;33(1):87–107. https://doi.org/10.1016/0304-3959(88)90209-6.
Ren K, Dubner R. Inflammatory models of pain and hyperalgesia. ILAR J. 1999;40(3):111–8. https://doi.org/10.1093/ilar.40.3.111.
Chaplan SR, Bach FW, Pogrel JW, Chung JM, Yaksh TL. Quantitative assessment of tactile allodynia in the rat paw. J Neurosci Methods. 1994;53(1):55–63. https://doi.org/10.1016/0165-0270(94)90144-9.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60. https://doi.org/10.1038/nmeth.3317.
Robinson JT, Thorvaldsdottir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29(1):24–6. https://doi.org/10.1038/nbt.1754.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30. https://doi.org/10.1093/bioinformatics/btt656.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. https://doi.org/10.1186/s13059-014-0550-8.
Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10(1):1523. https://doi.org/10.1038/s41467-019-09234-6.
Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nat Methods. 2013;10(12):1213–8. https://doi.org/10.1038/nmeth.2688.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://doi.org/10.1093/bioinformatics/btu170.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9. https://doi.org/10.1038/nmeth.1923.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9. https://doi.org/10.1093/bioinformatics/btp352.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2. https://doi.org/10.1093/bioinformatics/btq033.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38(4):576–89. https://doi.org/10.1016/j.molcel.2010.05.004.
Zhang Y, Chi D. Overexpression of SIRT2 Alleviates neuropathic pain and neuroinflammation through deacetylation of transcription factor nuclear factor-kappa B. Inflammation. 2018;41(2):569–78. https://doi.org/10.1007/s10753-017-0713-3.
Cetereisi D, Kramvis I, Gebuis T, van der Loo RJ, Gouwenberg Y, Mansvelder HD, et al. Gpr158 deficiency impacts hippocampal CA1 neuronal excitability, dendritic architecture, and affects spatial learning. Front Cell Neurosci. 2019;13:465. https://doi.org/10.3389/fncel.2019.00465.
Gasparini CF, Smith RA, Griffiths LR. Genetic and biochemical changes of the serotonergic system in migraine pathobiology. J Headache Pain. 2017;18(1):20. https://doi.org/10.1186/s10194-016-0711-0.
Lembo PM, Grazzini E, Groblewski T, O'Donnell D, Roy MO, Zhang J, et al. Proenkephalin A gene products activate a new family of sensory neuron—specific GPCRs. Nat Neurosci. 2002;5(3):201–9. https://doi.org/10.1038/nn815.
Liu F, Yajima T, Wang M, Shen JF, Ichikawa H, Sato T. Effects of trigeminal nerve injury on the expression of galanin and its receptors in the rat trigeminal ganglion. Neuropeptides. 2020;84:102098. https://doi.org/10.1016/j.npep.2020.102098.
Okumura Y, Yamagishi T, Nukui S, Nakao K. Discovery of AAT-008, a novel, potent, and selective prostaglandin EP4 receptor antagonist. Bioorg Med Chem Lett. 2017;27(5):1186–92. https://doi.org/10.1016/j.bmcl.2017.01.067.
Popova D, Desai N, Blendy JA, Pang ZP. Synaptic regulation by OPRM1 variants in reward neurocircuitry. J Neurosci. 2019;39(29):5685–96. https://doi.org/10.1523/JNEUROSCI.2317-18.2019.
Ramesh D, D'Agata A, Starkweather AR, Young EE. Contribution of endocannabinoid gene expression and genotype on low back pain susceptibility and chronicity. Clin J Pain. 2018;34(1):8–14. https://doi.org/10.1097/AJP.0000000000000508.
Sachau J, Bruckmueller H, Gierthmuhlen J, Magerl W, May D, Binder A, et al. The serotonin receptor 2A (HTR2A) rs6313 variant is associated with higher ongoing pain and signs of central sensitization in neuropathic pain patients. Eur J Pain. 2021;25(3):595–611. https://doi.org/10.1002/ejp.1696.
Zhao C, Quan X, He J, Zhao R, Zhang Y, Li X, et al. Identification of significant gene biomarkers of low back pain caused by changes in the osmotic pressure of nucleus pulposus cells. Sci Rep. 2020;10(1):3708. https://doi.org/10.1038/s41598-020-60714-y.
Zhou LL, Zhu YM, Qian FY, Yuan CC, Yuan DP, Zhou XP. MicroRNA1433p contributes to the regulation of pain responses in collageninduced arthritis. Mol Med Rep. 2018;18(3):3219–28. https://doi.org/10.3892/mmr.2018.9322.
Komori N, Takemori N, Kim HK, Singh A, Hwang SH, Foreman RD, et al. Proteomics study of neuropathic and nonneuropathic dorsal root ganglia: altered protein regulation following segmental spinal nerve ligation injury. Physiol Genomics. 2007;29(2):215–30. https://doi.org/10.1152/physiolgenomics.00255.2006.
Mayeux V, Valmier J. Skeletal muscle contraction modulates carbonic anhydrase phenotype in adult mouse dorsal root ganglion neurons. Brain Res. 1995;694(1-2):191–9. https://doi.org/10.1016/0006-8993(95)00698-P.
Hammer P, Banck MS, Amberg R, Wang C, Petznick G, Luo S, et al. mRNA-seq with agnostic splice site discovery for nervous system transcriptomics tested in chronic pain. Genome Res. 2010;20(6):847–60. https://doi.org/10.1101/gr.101204.109.
Hu G, Huang K, Hu Y, Du G, Xue Z, Zhu X, et al. Single-cell RNA-seq reveals distinct injury responses in different types of DRG sensory neurons. Sci Rep. 2016;6(1):31851. https://doi.org/10.1038/srep31851.
Pokhilko A, Nash A, Cader MZ. Common transcriptional signatures of neuropathic pain. Pain. 2020;161(7):1542–54. https://doi.org/10.1097/j.pain.0000000000001847.
Sun W, Kou D, Yu Z, Yang S, Jiang C, Xiong D, et al. A transcriptomic analysis of neuropathic pain in rat dorsal root ganglia following peripheral nerve injury. Neuromolecular Med. 2020;22(2):250–63. https://doi.org/10.1007/s12017-019-08581-3.
Athie MCP, Vieira AS, Teixeira JM, Dos Santos GG, Dias EV, Tambeli CH, et al. Transcriptome analysis of dorsal root ganglia's diabetic neuropathy reveals mechanisms involved in pain and regeneration. Life Sci. 2018;205:54–62. https://doi.org/10.1016/j.lfs.2018.05.016.
Dawes JM, Antunes-Martins A, Perkins JR, Paterson KJ, Sisignano M, Schmid R, et al. Genome-wide transcriptional profiling of skin and dorsal root ganglia after ultraviolet-B-induced inflammation. PLoS One. 2014;9(4):e93338. https://doi.org/10.1371/journal.pone.0093338.
Reinhold AK, Batti L, Bilbao D, Buness A, Rittner HL, Heppenstall PA. Differential transcriptional profiling of damaged and intact adjacent dorsal root ganglia neurons in neuropathic pain. PLoS One. 2015;10(4):e0123342. https://doi.org/10.1371/journal.pone.0123342.
Hatanaka Y, Kawasaki T, Abe T, Shioi G, Kohno T, Hattori M, et al. Semaphorin 6A-Plexin A2/A4 interactions with Radial glia regulate migration termination of superficial layer cortical neurons. iScience. 2019;21:359–74.
Suto F, Tsuboi M, Kamiya H, Mizuno H, Kiyama Y, Komai S, et al. Interactions between plexin-A2, plexin-A4, and semaphorin 6A control lamina-restricted projection of hippocampal mossy fibers. Neuron. 2007;53(4):535–47. https://doi.org/10.1016/j.neuron.2007.01.028.
Nogi T, Yasui N, Mihara E, Matsunaga Y, Noda M, Yamashita N, et al. Structural basis for semaphorin signalling through the plexin receptor. Nature. 2010;467(7319):1123–7. https://doi.org/10.1038/nature09473.
Renaud J, Kerjan G, Sumita I, Zagar Y, Georget V, Kim D, et al. Plexin-A2 and its ligand, Sema6A, control nucleus-centrosome coupling in migrating granule cells. Nat Neurosci. 2008;11(4):440–9. https://doi.org/10.1038/nn2064.
Gallagher SJ, Mijatov B, Gunatilake D, Gowrishankar K, Tiffen J, James W, et al. Control of NF-kB activity in human melanoma by bromodomain and extra-terminal protein inhibitor I-BET151. Pigment Cell Melanoma Res. 2014;27(6):1126–37. https://doi.org/10.1111/pcmr.12282.
Jiang BC, He LN, Wu XB, Shi H, Zhang WW, Zhang ZJ, et al. Promoted interaction of C/EBPalpha with demethylated Cxcr3 gene promoter contributes to neuropathic pain in mice. J Neurosci. 2017;37(3):685–700. https://doi.org/10.1523/JNEUROSCI.2262-16.2016.
Klemm SL, Shipony Z, Greenleaf WJ. Chromatin accessibility and the regulatory epigenome. Nat Rev Genet. 2019;20(4):207–20. https://doi.org/10.1038/s41576-018-0089-8.
Fernández-Montoya J, Avendaño C, Negredo P. The glutamatergic system in primary somatosensory neurons and is involvement in sensory input-dependent plasticity. Int J Mol Sci. 2018;19:69. https://doi.org/10.3390/ijms19010069.
Aller MI, Pecoraro V, Paternain AV, Canals S, Lerma J. Increased dosage of high-affinity Kainate receptor gene grik4 alters synaptic transmission and reproduces autism spectrum disorders features. J Neurosci. 2015;35(40):13619–28. https://doi.org/10.1523/JNEUROSCI.2217-15.2015.
Zhang Y, Liu T, Meyer CA, et al. Model-based Analysis of ChIP-Seq (MACS). Genome Biol. 2008;9:R137. https://doi.org/10.1186/gb-2008-9-9-r137.
We thank Rakel Tryggvadóttir and Dr. Jaclyn Daniels for their technical assistance.
This study was supported by grants from the National Institutes of Health (Bethesda, Maryland, USA) F32NR015728 (KES), KL2TR003108 (KES), NS110598 (YG), NS117761 (YG), R01GM118760 (SDT), the Arkansas Children’s Research Institute (KES), the Arkansas Breast Cancer Research Program (KES) as well as a seed grant from the Johns Hopkins Blaustein Pain Research Fund (SDT). The funding bodies had no role in the design of the study and collection, analysis, interpretation of data, and in the writing of the manuscript.
Ethics approval and consent to participate
All procedures were reviewed and approved by the Johns Hopkins Animal Care and Use Committee and were performed in accordance with the National Institutes of Health’s Guide for the Care and Use of Laboratory Animals.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Consistency of peaks among replicates. The graph shows the number of peaks identified in one or more samples. Due to the large number of accessible regions identified in only one sample, we included only regions that were identified in at least 50% of the samples within the study group.
Genomic locations and PCR primers for luciferase assay constructs.
Genes differentially expressed in the CCI and CFA models.
Genes encoding transcriptional regulators differentially expressed in both CCI and CFA models.
Transcriptional factor motifs significantly enriched in differentially accessible regions after CCI.
Transcriptional factor motifs significantly enriched in differentially accessible regions after CFA.
About this article
Cite this article
Stephens, K.E., Zhou, W., Renfro, Z. et al. Global gene expression and chromatin accessibility of the peripheral nervous system in animal models of persistent pain. J Neuroinflammation 18, 185 (2021). https://doi.org/10.1186/s12974-021-02228-6
- Chromatin accessibility
- Dorsal root ganglion
- Nerve injury