Skip to content

Advertisement

  • Research article
  • Open Access

Identification of key genes and pathways involved in response to pain in goat and sheep by transcriptome sequencing

  • 1, 2,
  • 3,
  • 1,
  • 2 and
  • 1Email author
Contributed equally
Biological Research201851:25

https://doi.org/10.1186/s40659-018-0174-7

  • Received: 16 April 2018
  • Accepted: 7 August 2018
  • Published:

Abstract

Purpose

This aim of this study was to investigate the key genes and pathways involved in the response to pain in goat and sheep by transcriptome sequencing.

Methods

Chronic pain was induced with the injection of the complete Freund’s adjuvant (CFA) in sheep and goats. The animals were divided into four groups: CFA-treated sheep, control sheep, CFA-treated goat, and control goat groups (n = 3 in each group). The dorsal root ganglions of these animals were isolated and used for the construction of a cDNA library and transcriptome sequencing. Differentially expressed genes (DEGs) were identified in CFA-induced sheep and goats and gene ontology (GO) enrichment analysis was performed.

Results

In total, 1748 and 2441 DEGs were identified in CFA-treated goat and sheep, respectively. The DEGs identified in CFA-treated goats, such as C-C motif chemokine ligand 27 (CCL27), glutamate receptor 2 (GRIA2), and sodium voltage-gated channel alpha subunit 3 (SCN3A), were mainly enriched in GO functions associated with N-methyl-d-aspartate (NMDA) receptor, inflammatory response, and immune response. The DEGs identified in CFA-treated sheep, such as gamma-aminobutyric acid (GABA)-related DEGs (gamma-aminobutyric acid type A receptor gamma 3 subunit [GABRG3], GABRB2, and GABRB1), SCN9A, and transient receptor potential cation channel subfamily V member 1 (TRPV1), were mainly enriched in GO functions related to neuroactive ligand-receptor interaction, NMDA receptor, and defense response.

Conclusions

Our data indicate that NMDA receptor, inflammatory response, and immune response as well as key DEGs such as CCL27, GRIA2, and SCN3A may regulate the process of pain response during chronic pain in goats. Neuroactive ligand-receptor interaction and NMDA receptor as well as GABA-related DEGs, SCN9A, and TRPV1 may modulate the process of response to pain in sheep. These DEGs may serve as drug targets for preventing chronic pain.

Keywords

  • Chronic pain
  • Transcriptome sequencing
  • Differentially expressed genes
  • Gene ontology
  • Goat
  • Sheep

Background

Chronic pain is considered as a major physical and mental health problem. Millions of people are affected by uncomfortable conditions such as back pain, headache, and arthritis [1, 2]. Chronic pain is often difficult to treat [3, 4], and 60% of patients with chronic pain experience pain after 1 year of treatment [5]. Although extensive research has been conducted on the genetics of chronic pain, the mechanism underlying chronic pain is largely unknown owing to the involvement of several genes [6]. For a better understanding of the etiology and treatment of chronic pain, it is imperative to further elucidate the key mechanism.

Dorsal root ganglion (DRG) contains many primary sensory neurons that are responsible for the transduction and modulation of sensory information and its transmission to the spinal cord, an active participant in the development of chronic neuropathic pain [7]. Accumulating evidence has confirmed that DRG plays a critical role in the induction and maintenance of chronic pain [8]. DRG is involved in the transduction of pain to the central nervous system and exhibits various pathophysiologic changes during chronic pain [9]. Moreover, DRG stimulation in response to peripheral afferent fiber injury may induce multiple abnormal changes that occur within the DRG and stabilize or reduce the hyperexcitability of DRG neurons and consequently decrease chronic neuropathic pain [10]. Furthermore, DRG stimulation is considered as a promising treatment strategy to offer relief from chronic pain, and DRG therapeutics have been implicated for the treatment of chronic pain [11, 12]. KCNQ channels are found expressed in nociceptive DRG neurons and their activation is effective in reducing chronic pain [13]. The differential expression of ATP-gated P2X receptors in DRG is shown to play important roles in the regulation of nociceptive mechanisms in chronic neuropathic pain and visceralgia rat models [14]. The transient receptor potential melastatin 8 in DRG plays a key role in chronic neuropathic pain [15]. Microarray analysis of DRG gene expression has been used for the identification of key cytokines involved in the regulation of chronic pain [16]. Despite these advances, the regulatory mechanisms mediating the development of chronic pain are largely unknown. Given that the spinal DRG is the primary center for the conduction and maintenance of pain, studies on DRG transcriptome may be helpful for the comprehensive and systematic elucidation of the key regulatory mechanisms underlying chronic pain.

In the present study, chronic pain was induced by inoculation of the complete Freund’s adjuvant (CFA) in sheep and goats, as animals injected with CFA are known to quickly develop hyperalgesia and have low latency times on their inflamed paws [17]. A clear advantage of using a large animal model such as sheep is the similarity in the body size, spine dimensions, and cardiac and pulmonary functional parameters to humans [18, 19]. Another benefit of using sheep as a chronic pain model is the longer lifespan of sheep (approximately 20 years) than rats (about 2–4 years) [20]. Previous studies of sheep foot rot have revealed a chronic inflammatory pain condition [21, 22]. Our previous observation on the differences in responses to pain between sheep and goats encouraged us to select sheep and goats as the animal models. After CFA inoculation, we analyzed the pain response mechanism in goats and sheep through transcriptome sequencing using the DRG tissues and performed bioinformatic analysis. Our findings may improve the understanding of the molecular mechanisms underlying chronic pain and provide valuable information for the development of treatment strategy.

Methods

Animal preparation, treatment, and grouping

All experiments were approved by the local animal care and use committee. Six sheep and six goats weighing 30–40 kg and 2- to 3-year old were obtained from the key laboratory of Inner Mongolia Autonomous Region (Hohhot, Inner Mongolia, China) and housed in the central housing facility under a standard condition. Food and water were available ad libitum. The CFA-induced arthritic rat model has been used to study chronic pain [23]; hence, chronic pain was induced by CFA inoculation in the sheep and goats. CFA (Sigma, USA) is a mixture of heat-killed Mycobacterium tuberculosis and paraffin oil at a concentration of 1 mg/mL. Prior to experiments, CFA was diluted 1:1 in 0.9% sterile saline. Each sheep or goat was injected in the left footpad (hind paw) with 1 mL of CFA under isoflurane anesthesia. This dose was chosen for the following reasons: After the injection of 1 mL of CFA, pain responses such as redness, swelling, and hyperalgesia were markedly observed at the injection site in both sheep and goats, indicating that the inflammatory pain was successfully induced by CFA; however, the induction effect was not satisfactory under or above this dose of CFA. Sheep or goats injected with 1 mL of saline served as controls. The animals were divided into four groups, namely, CFA-treated sheep (st1-3), control sheep (sc1-3), CFA-treated goat (gt1-3), and control goat (gc1-3) groups (n = 3 in each group).

Isolation of DRG

After 72 h from CFA inoculation, the animals were sequentially slaughtered under the anesthesia of Deer Sleeping Ling. This time point was selected for the following reasons: (1) From 6 to 72 h after the injection of CFA, the redness and swelling at the injection sites in the left footpads of goats and sheep gradually aggravated. The animals showed signs of severe pain with different degrees of functional abnormalities such as tiptoe, lameness, and inability to stand at 72 h of injection, and these symptoms began to ease after 96 h of injection. (2) The time point of 72 h is long enough for the development of changes in the transcription/translation of proteins that are involved in the sensory neuron response to inflammation [24]. The lumbar spinal bone was placed on the ice and both sides of the muscles were removed under sterile conditions. The spinous process and transverse process of spines were fully exposed and cut-off for the isolation of the spinal cord. At the spinal canal of the posterior root of the spinal nerve, DRG was stripped-off and the L3-L5 DRG was collected. After rinsing in the RNA protection solution, the L3-L5 DRG was quickly stored in liquid nitrogen.

Construction of the cDNA library of DRGs and high-throughput sequencing

Total RNA was extracted from DRGs using TRIzol regents and treated with RQ1 DNase (Promega) to remove DNA. The concentration and purity of the extracted RNA was determined by measuring the ratio of the absorbance at 260 and 280 nm (A260/A280) using SmartSpec plus (Bio-Rad Laboratories, Hercules, Calif). The integrity of the extracted RNA was confirmed with 1.5% agarose gel electrophoresis. A total of 10 μg of total RNA was used for RNA-seq library preparation. For directional RNA-seq library preparation, polyadenylated mRNAs were purified with oligo(dT)-conjugated magnetic beads (Invitrogen) and iron fragmented at 95 °C. After end repair and 5′ adaptor ligation, the cDNAs were reverse transcribed, purified, and amplified. The PCR products (200–500 bp) were finally purified and stored at − 80°C for subsequent high-throughput sequencing. The cDNA libraries were prepared and applied onto Illumina Hiseq 2000 system for 100 nucleotide pair-end sequencing by Majorbio. Inc (Shanghai, China).

Preprocessing and quality control of sequencing data

Pre-processing for the retrieval of the clean reads was performed as follows: The reads with two N were excluded; the adapter sequences were removed according to the joint information; reads with low Q < 16 were removed. The quality assessment of the raw reads was conducted using FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

Read mapping and transcription annotation

The clean reads were mapped to the reference genome of sheep (ftp://ftp.ncbi.nlm.nih.gov/genomes/Ovis_aries/) and goat (http://goat.kiz.ac.cn/GGD) using TopHat2.0.1 [25] with default parameters. Bowtie was used during mapping, and no more than one mismatched read was allowed. For calculating the gene expression, the mRNA length and sequencing depth were homogenized using reads per kilo base of a gene per million reads (RPKM) algorithm.

Correlation analysis for the gene expression levels in different samples

Correlation analysis for gene expression levels between any two samples was performed for the three duplicate samples of sheep/goat under the same treatment condition to avoid any individual differences. If the correlation coefficient (R) was close to 1, the expression patterns between different samples were considered to have higher similarity, showing that the data had high degree of homogenization and were reliable for the subsequent analysis. In our study, the sample was deleted if R2 > 0.8. To further avoid any individual difference, clustering heatmap analysis for the RPKM data was conducted using the pheatmap version 1.0.8 (http://cran.r-project.org/web/packages/pheatmap/index.html) in R3.4.1.

Identification of differentially expressed genes (DEGs)

We identified DEGs in the CFA-treated sheep group and the control sheep group using the likelihood ratio tests in edgeR package [26] in R language. The cut-off thresholds were fold change (FC) ≥ 2 or < 0.5 and P ≤ 0.01. Moreover, all gene sequences of goats were annotated by blasting to the sheep genome sequences. The DEGs in the CFA-treated goat group were compared with those in the control goat group using the edgeR package with the same thresholds. Moreover, the common DEGs (co-DEGs) between the CFA-treated goat and sheep groups were identified and visualized by Venn diagram.

Gene ontology (GO) enrichment analysis

GO (http://www.geneontology.org) [27] is used for exploring the functions of large-scale genomic or transcriptomic data, and mainly includes biological process (BP), molecular function (MF), and cellular component (CC). GO enrichment analyses for DEGs in different groups were performed using an online tool database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov/) [28]. To further obtain the significant GO functions, the enriched GO terms were clustered using the functional annotation clustering tool in DAVID, wherein the enrichment score for each cluster was calculated as the negative log of the geometric mean of P values in the cluster. Each functional cluster contained GO functions with similar biological meaning, owing to the presence of similar gene members. The more the enrichment score, the more significant was the cluster. Enrichment score > 0.5 was set as the cut-off value.

Results

High-throughput sequencing data

The results showed that the A260/A280 of extracted total RNA from all DRGs samples was 1.8–2.1, indicating that total RNA from all DRGs samples could be used for construction of the cDNA library of DRGs and subsequent high-throughput sequencing. Using the Illumina sequencing technology, more than 41,000,000 clean reads were obtained from each sample. The GC content of these 12 samples was all less than 50 and Q30(%) of these samples was about 90%, indicating that the obtained data were good and reliable for the subsequent analysis.

According to the sequencing data, a total of 20,951 genes were identified from six goat samples, accounting for 97.95% (20,951/21,389) of all goat genes. In addition, 22,066 genes were identified from six sheep samples, accounting for 88.15% (22,066/25,033) of all sheep genes. These data indicate that most of genes were expressed in DRG tissues.

Correlation analysis for the gene expression levels in different samples

To avoid the individual difference, correlation analysis for gene expression levels between any two samples was performed (data not shown). The results revealed the high correlation between any two samples in the control goat or/and CFA-treated goat groups (all R2 > 0.92), indicating that the individual difference was small. Furthermore, high correlation was observed between any two samples in the control sheep or/and CFA-treated sheep groups (all R2 > 0.8). High correlations were also observed between control goat and sheep samples as well as between CFA-treated goat and sheep samples (all R2 > 0.89), indicating that the expression patters of common genes in goat and sheep were similar. Heatmap analysis results showed that the control and CFA-treated goat samples were clustered (Fig. 1a), indicating that CFA treatment induced changes in the expression of goat genes. In addition, sc3 were clustered together with three CFA-treated sheep samples (Fig. 1b), suggesting that sc3 samples had larger individual differences with other two control sheep samples. The sc3 sample was detected and was not used for subsequent analyses. After the detection of the sc3 sample, the results of the hierarchical clustering analysis showed that the control and CFA-treated sheep samples were clustered (Fig. 1c).
Fig. 1
Fig. 1

Heatmap analysis of goat (a) and sheep (b, c) samples. The below abscissa axis represents samples. Red shows upregulated genes, while green shows downregulated genes. In (c) sc3 sample was detected

Identification of DEGs

Using the edgeR package, a total of 1748 DEGs (741 upregulated and 1007 downregulated) were identified in the CFA-treated goat group as compared with the control goat group. In addition, 2441 DEGs (993 upregulated and 1508 downregulated) were identified in the CFA-treated sheep group as compared with the control sheep group. A total of 471 common DEGs (132 upregulated and 339 downregulated) were screened out from the CFA-treated goat and sheep groups as compared with the control sheep and goat groups (Fig. 2). Furthermore, 2953 DEGs (1762 upregulated and 1173 downregulated) were identified between the control sheep and goat groups.
Fig. 2
Fig. 2

Venn diagram showed the common differentially expressed genes (co-DEGs) between the CFA-treated goat and sheep groups as compared with the control sheep and goat groups

Functional analysis for DEGs

To better understand the function of DEGs, GO functional enrichment analysis was performed. DEGs identified in the CFA-treated goats were mainly enriched in GO functions associated with gated channel activity, N-methyl-d-aspartate (NMDA) receptor, inflammatory response, immune response, and defense response (Table 1 and Fig. 3). Key DEGs such as C-C motif chemokine ligand 27 (CCL27), glutamate receptor 2 (GRIA2), glutamate ionotropic receptor NMDA type subunit 2A (GRIN2A), calcium voltage-gated channel subunit alpha 1E (CACNA1E), sodium voltage-gated channel alpha subunit 3 (SCN3A), and potassium two pore domain channel subfamily K member 1 (KCNK1) were enriched in these GO functions and may play key roles in the process of response to pain in goats (Table 1). Moreover, the DEGs identified in the CFA-treated sheep were mainly enriched in GO functions related to gated channel activity, neuroactive ligand-receptor interaction, NMDA receptor, voltage-gated potassium channel complex, and defense response (Table 2 and Fig. 4). Many gamma-aminobutyric acid (GABA)-related DEGs such as gamma-aminobutyric acid type A receptor gamma 3 subunit (GABRG3), gamma-aminobutyric acid type A receptor beta 2 subunit (GABRB2), and gamma-aminobutyric acid type A receptor beta 1 subunit (GABRB1) were markedly upregulated. In addition, SCN9A and transient receptor potential cation channel subfamily V member 1 (TRPV1) were found markedly downregulated. These DEGs may be involved in the process of response to pain in sheep. GO clusters enriched by DEGs identified between CFA-treated animals and control animals are shown in Additional file 1: Table S1.
Table 1

The significantly enriched GO terms associated with pain in goats

Go term

Count

P value

Gene

Up-regulated

 GO:0006811 ~ ion transport

33

1.69E−04

SCN3A, SLC38A8, GABRB1, GRIK3, CNGB1, SLCO1A2, GRID2, SLC4A4, KCNG3, CHRNA1, CLCA1, GABRA4, TRPC5, SLC12A3, GRIN2A, SLC22A20, KCNK1, CACNA2D4, ITPR2, SLCO1B3, GRIA2, SLC17A4, SLC26A7, CACNA1E, KCNH8, ATP7B, KCNH5

 GO:0022836 ~ gated channel activity

18

3.76E−04

GABRA4, SCN3A, TRPC5, GABRB1, GRIK3, GRIN2A, CNGB1, KCNK1, ITPR2, CACNA2D4, GRIA2, GRID2, KCNH8, CACNA1E, KCNG3, CHRNA1, FGF2, KCNH5

 GO:0022834 ~ ligand-gated channel activity

11

3.98E−04

GRIA2, GABRA4, TRPC5, GABRB1, GRIK3, GRID2, GRIN2A, CNGB1, CHRNA1, KCNK1, ITPR2

 GO:0005216 ~ ion channel activity

20

6.64E−04

CLCA1, GABRA4, SCN3A, TRPC5, GRIK3, GABRB1, GRIN2A, CNGB1, KCNK1, ITPR2, CACNA2D4, GRIA2, SLC26A7, GRID2, KCNH8, CACNA1E, KCNG3, CHRNA1, FGF2, KCNH5

 GO:0022832 ~ voltage-gated channel activity

9

0.058679

SCN3A, GRIN2A, CACNA1E, KCNH8, KCNG3, KCNK1, FGF2, KCNH5, CACNA2D4

 GO:0008066 ~ glutamate receptor activity

5

0.004105

GRM3, GRIA2, GRIK3, GRID2, GRIN2A

 IPR001508:NMDA receptor

4

0.005749

GRIA2, GRIK3, GRID2, GRIN2A

 GO:0006816 ~ calcium ion transport

8

0.03123

CLCA1, NPY, TRPC5, GRIN2A, CACNA1E, JPH1, ITPR2, CACNA2D4

 GO:0006955 ~ immune response

18

0.291463

CXCL1, PSMB10, AQP9, S100A7, C6, CXCL9, CD180, PDCD1LG2, CCL27, TNFSF8, APOA4, C8A, CD96, OASL, CFHR5, VSIG4, IL1RAPL1, IFI6

 GO:0009611 ~ response to wounding

13

0.432142

CXCL1, SELP, C6, MAP1B, CXCL9, GRIN2A, CD180, C8A, FGA, SERPINB2, VSIG4, CFHR5, FGF2

 GO:0006954 ~ inflammatory response

8

0.526807

C8A, CXCL1, SELP, C6, CXCL9, CFHR5, VSIG4, CD180

 GO:0006952 ~ defense response

14

0.540829

CXCL1, SELP, S100A7, C6, CXCL9, FOXN1, VARS, CD180, CLEC1A, APOA4, C8A, VSIG4, CFHR5, IL1RAPL1

Down-regulated

 IPR001073:Complement C1q protein

8

5.07E−04

CBLN4, C1QTNF3, C1QTNF1, C1QTNF2, EMILIN2, MMRN1, ADIPOQ, COL8A2

 GO:0009611 ~ response to wounding

46

8.77E−05

GGCX, ERBB3, F13A1, TGFB3, AFAP1L2, MMRN1, SRF, TIMP3, SHH, AZU1, CCL24, IGSF10, S1PR3, HRH1, CASP3, CCL21, LTB4R, LOX, THBS1, FN1, RAB27A, INA, KLF6, CR2, LIPA, GATM, CYP1A1, PDPN, SAAL1, ANXA1, CHST3, COL5A1, NOTCH3, TNFAIP6, SDC1, LYVE1, THBD, TFRC, ADM, SERPINF2, ITGA5, PDGFRA, HBEGF, VCAN, IGFBP4, AOC3

 GO:0006954 ~ inflammatory response

19

0.270975

LIPA, CR2, PDPN, SAAL1, ANXA1, AFAP1L2, CCL24, AZU1, S1PR3, TNFAIP6, HRH1, TFRC, LTB4R, CCL21, SERPINF2, THBS1, IGFBP4, AOC3, FN1

 GO:0006952 ~ defense response

25

0.855432

AFAP1L2, CCL24, AZU1, S1PR3, HRH1, LTB4R, CCL21, THBS1, FN1, RAB27A, CR2, LIPA, PDPN, SAAL1, ANXA1, COLEC12, INHBB, TNFAIP6, INHBA, BPI, TFRC, SERPINF2, IGFBP4, ICOSLG, AOC3

 GO:0045202 ~ synapse

27

0.015927

CDK5R1, SNCAIP, SYT4, ERBB3, SYT6, SV2B, PRKACA, SV2A, APBA1, GABRG2, EFNB1, SYT11, RIMBP2, SYT12, NLGN1, LIN7C, PPP1CC, RPH3A, SSPN, CBLN4, EGFLAM, COLQ, MAP1S, GRM8, CHRNB3, VAMP3, DOC2B

 GO:0044456 ~ synapse part

16

0.180734

GABRG2, SNCAIP, SYT4, ERBB3, SYT11, NLGN1, LIN7C, SYT6, RPH3A, SSPN, GRM8, CHRNB3, SV2B, SV2A, DOC2B, APBA1

 GO:0045596 ~ negative regulation of cell differentiation

19

0.01342

TWSG1, SIX2, SOCS5, JAG1, LIG4, GLI2, SOX9, ADIPOQ, SHH, HES1, NOTCH3, INHBA, RNF6, HOXA2, RHOA, TGFBR3, TOB2, BMPR1A, TOB1

 IPR001565:Synaptotagmin

5

0.005111

SYT4, SYT11, SYT6, RPH3A, DOC2B

 GO:0005544 ~ calcium-dependent phospholipid binding

4

0.115606

SYT11, ANXA1, DOC2B, RBM12

 GO:0045664 ~ regulation of neuron differentiation

12

0.048765

NOTCH3, AMIGO1, HES1, HOXA2, CDK5R1, RNF6, NLGN1, RHOA, SEMA4D, GLI2, SHH, KALRN

Fig. 3
Fig. 3

The significantly enriched gene ontology terms associated with pain in goats

Table 2

The significantly enriched GO terms associated with pain in sheep

Go term

Count

P value

Gene

Up-regulated

 GO:0015276 ~ ligand-gated ion channel activity

14

2.48E−05

GABRG3, SHROOM2, GRIK1, GABRB2, GLRA3, GABRB1, GRIN2A, GRIA4, CNGB3, KCNJ6, GRID2, CHRNA5, KCNH7, SCNN1B

 GO:0006811 ~ ion transport

34

2.10E−03

SLC5A5, GRIK1, GABRB2, ATOX1, GABRB1, GLRA3, SLC39A12, ATP5G1, CNGB3, CHRNA5, GRID2, SLC30A3, SCNN1B, ATP5I, OCA2, SLC30A7, TRPM3, GABRG3, KCND2, KCNB2, SLC12A3, GRIN2A, GRIA4, CACNG2, TCN1, ATP13A4, SLCO1B3, KCNJ6, CATSPER3, KCNH7, HEPH, KCTD16, SLC13A4, CLCN7

 GO:0022836 ~ gated channel activity

20

4.61E−04

GABRG3, SHROOM2, KCND2, GRIK1, GABRB2, KCNB2, GLRA3, GABRB1, GRIN2A, CACNG2, GRIA4, CNGB3, KCNJ6, CATSPER3, GRID2, CHRNA5, KCNH7, KCTD16, SCNN1B, CLCN7

 GO:0005216 ~ ion channel activity

21

0.002572

TRPM3, GABRG3, SHROOM2, KCND2, GRIK1, GABRB2, KCNB2, GLRA3, GABRB1, GRIN2A, CACNG2, GRIA4, CNGB3, KCNJ6, CATSPER3, GRID2, CHRNA5, KCNH7, KCTD16, SCNN1B, CLCN7

 GO:0060284 ~ regulation of cell development

4

0.895878

MUSK, MAP1B, NLGN1, SOX5

 GO:0007268 ~ synaptic transmission

15

0.020093

GABRG3, KCND2, MYO6, GRIK1, GABRB2, GLRA3, PCDHB4, NLGN1, GRIN2A, GRIA4, MUSK, GRID2, CHRNA5, GAD1, NPFF

 GO:0019226 ~ transmission of nerve impulse

16

0.033585

GABRG3, KCND2, MYO6, GRIK1, GABRB2, GLRA3, PCDHB4, NLGN1, GRIN2A, CACNG2, GRIA4, MUSK, GRID2, CHRNA5, GAD1, NPFF

 GO:0006812 ~ cation transport

22

0.042508

TRPM3, SLC5A5, KCND2, KCNB2, SLC12A3, ATOX1, SLC39A12, GRIN2A, CACNG2, ATP5G1, TCN1, ATP13A4, KCNJ6, CATSPER3, KCNH7, HEPH, KCTD16, SLC30A3, SLC13A4, SCNN1B, ATP5I, SLC30A7

 IPR001508:NMDA receptor

4

0.009544

GRIK1, GRID2, GRIN2A, GRIA4

 GO:0005254 ~ chloride channel activity

5

0.116044

GABRG3, GABRB2, GABRB1, GLRA3, CLCN7

 GO:0005244 ~ voltage-gated ion channel activity

9

0.134055

KCND2, KCNJ6, CATSPER3, KCNB2, GRIN2A, KCNH7, CACNG2, KCTD16, CLCN7

 GO:0031404 ~ chloride ion binding

5

0.138619

GABRG3, GABRB2, GABRB1, GLRA3, CLCN7

Down-regulated

 GO:0005248 ~ voltage-gated sodium channel activity

9

2.43E−06

SCN2B, SCN3B, PKD2, SCN9A, SCN4B, SCN11A, SCN8A, SCN7A, SCN10A

 GO:0001518 ~ voltage-gated sodium channel complex

7

6.00E−05

SCN2B, SCN9A, SCN4B, SCN11A, SCN8A, SCN7A, SCN10A

 GO:0031402 ~ sodium ion binding

21

1.81E−04

SLC8A3, SLC13A5, PDXK, SCN2B, SLC12A2, SCN3B, ATP1B2, SLC9A3, KIAA1919, SLC34A2, SLC17A6, SLC6A8, SCN9A, SLC4A8, SCN4B, SLC13A3, SCN11A, SCN7A, SCN8A, HCN4, SCN10A

 GO:0022843 ~ voltage-gated cation channel activity

19

0.015225

KCNJ16, KCNAB2, SCN2B, SCN3B, KCNA2, KCTD20, KCNJ5, KCTD9, KCNK6, SCN9A, PKD2, SCN4B, SCN11A, SCN7A, SCN8A, KCNQ2, HCN4, KCTD12, SCN10A

 hsa04080:neuroactive ligand-receptor interaction

10

0.978027

GABRG2, GABRE, HTR1B, PTGER3, TRPV1, P2RY1, LTB4R2, AVPR1A, GABBR2, PTAFR

 GO:0005083 ~ small GTPase regulator activity

36

4.43E−04

ALS2, TBC1D9, CYTH1, PREX1, MYO9B, CYTH3, ARFGEF2, RIMS1, PLEKHG3, PLEKHG2, DMXL2, PLEKHG7, SOS1, GOPC, ARHGAP1, TBC1D13, PLEKHG5, RAPGEF2, FGD5, ERRFI1, IQSEC1, TBC1D2B, ABR, TRIO, RPH3A, VAV2, ARHGEF10, MAP4K4, TBC1D24, CDC42BPG, ADAP2, SGSM1, SGSM2, CIT, CDC42BPB, RANBP10

 GO:0008076 ~ voltage-gated potassium channel complex

7

0.526697

KCNJ5, KCTD9, KCNK6, KCNA2, KCTD20, KCNQ2, KCTD12

 GO:0006952 ~ defense response

33

0.971848

ALS2, YWHAZ, C3, TRPV1, TLR2, AFAP1L2, CX3CL1, ITGB1, TLR7, IL31RA, FCN3, SCN9A, THBS1, MX2, FN1, SPP1, PTGER3, LIPA, IL1RL1, SAAL1, SOCS6, SMAD1, CD1D, DDX58, SIGLEC1, ABCC9, SARM1, CCR5, TFRC, CCR2, ICOSLG, CD14, PTAFR

Fig. 4
Fig. 4

The significantly enriched gene ontology terms associated with pain in sheep

We found that the DEGs only identified in goat, such as immunoglobulin superfamily member 6 (IGSF6), C-X-C motif chemokine ligand 11 (CXCL11), pentraxin 3 (PTX3), C-C motif chemokine ligand 17 (CCL17), and chemokine (C-C motif) receptor-like 1 (CCRL1) were significantly enriched in chemokine signaling pathway, chemotaxis, inflammatory response, cytokine–cytokine receptor interaction, and immune response (Additional file 1: Table S2), while the DEGs only identified in sheep, including colony-stimulating factor 2 (CSF2), CCAAT/enhancer binding protein gamma (CEBPG), RELB proto-oncogene, nuclear factor kappa B (NF-KB) subunit (RELB), B cell CLL/lymphoma 3 (BCL3), and integrin subunit beta 1 (ITGB1), were markedly enriched in intermediate filament, leukocyte differentiation, B cell activation, cellular response to stress, cytokine–cytokine receptor interaction, and cell–cell signaling (Additional file 1: Table S2).

We also performed GO analysis of co-DEGs between goat and sheep. The results showed that the upregulated co-DEGs were significantly enriched in embryonic organ development, transmembrane protein, cytoskeleton, neuroactive ligand-receptor interaction, and synapse, while the downregulated co-DEGs were enriched in cell adhesion, extracellular matrix, extracellular matrix-receptor interaction, cell motion, and signal peptide (Additional file 1: Table S3).

Discussion

In the present study, we investigated the mechanism underlying the response to pain in goat and sheep through transcriptome sequencing using DRG tissues of CFA-induced sheep and goats and subsequent bioinformatics analysis. The results showed that the DEGs identified in CFA-treated goats, such as CCL27, GRIA2, and SCN3A, were mainly enriched in GO functions associated with NMDA receptor, inflammatory response, and immune response and may play key roles in the process of response to pain in goats. The DEGs identified in CFA-treated sheep, such as GABA-related DEGs (GABRG3, GABRB2, and GABRB1), SCN9A, and TRPV1, were mainly enriched in GO functions related to neuroactive ligand-receptor interaction, NMDA receptor, and defense response and may be involved in the process of response to pain in sheep. These data elucidated the possible mechanism involved in the process of response to pain in sheep and goat and may suggest potential drug targets for the treatment of chronic pain.

In a previous study, NMDA receptor was shown to play a key role in the neurophysiology of pain transmission and modulation [29]. In the insular cortex, the increase in synaptic NMDA receptors contributes to the development of neuropathic pain [30]. Wang et al. demonstrated that Ht31 peptide repressed inflammatory pain in mice by blocking the NMDA receptor-mediated nociceptive transmission [31]. Moreover, inflammatory response is also considered as the key mechanism involved in the regulation of neuropathic pain, and modulation of this process may serve as a treatment strategy for neuropathic pain [32]. The key DEGs associated with pain in goat were identified in the present study. CCL27 is a cytokine involved in chronic pain [33]. In a CFA-induced chronic inflammatory pain model, mice with null mutation of GRIA2 have increased thermal and mechanical hyperalgesia [34]. SCN3A is a voltage-gated sodium channel gene that regulates the important role of miR-30b in spinal nerve ligation-induced neuropathic pain in rats [35]. Given the key role of these GO functions and their enriched genes, we speculate that NMDA receptor and inflammatory response as well as the above mentioned enriched DEGs may be the key mechanisms that mediate the response to pain in goat.

We found that the NMDA receptor is a key GO function enriched by DEGs identified in CFA-treated sheep, indicative of its involvement in the process of response to pain in sheep. In addition, DEGs identified in CFA-treated sheep were mainly enriched in GO functions related to neuroactive ligand-receptor interaction. Sprangers et al. confirmed the association of neuroactive ligand-receptor interaction with disease-related pain and pain severity, duration, and relief [35], suggestive of the key role of this pathway in the process of response to pain. Moreover, GABA-related DEGs such as GABRG3, GABRB2, and GABRB1 were markedly upregulated. The inhibitory transmitter GABA has been shown to play a crucial role in the modulation of pain transmission [36, 37]. Gwak and Hulsebosch revealed that the spinal cord injury-induced hypofunction of GABAergic tone was involved in enhanced synaptic transmission, which results in neuronal hyperexcitability in dorsal horn neurons and subsequent chronic neuropathic pain [38]. These findings imply that the inhibitory transmitter GABA may serve as a key mediator involved in the response to pain in sheep. We also observed downregulation in SCN9A and TRPV1 expressions. Previous studies have confirmed that the mutations in SCN9A (NaV1.7) may cause pain disorders and pain sensitivity [39, 40]. In addition, TRPV1 is a molecular sensor of heat and capsaicin that was shown to regulate chronic pain by modulating central terminal sensitization [41], and targeting TRPV1 is considered as a promising strategy for the treatment of chronic pain [42]. Given the key role of SCN9A and TRPV1 in chronic pain modulation, we speculate that the downregulation of SCN9A and TRPV1 may contribute to the process of response to pain in sheep.

Some limitations of the present study include the small sample size owing to the limitation of experimental cost that may influence the stability of statistical power. In addition, the identified DEGs and pathways were not verified by functional experiments. Moreover, CFA is used as a co-adjuvant in antibody production. Whether CFA acts as an immunopotentiator and participates in the immune response to pain is unclear. Therefore, further investigations with additional experiments and high throughput data are warranted to confirm the findings of the present study.

Conclusion

In conclusion, our data indicate that NMDA receptor, inflammatory response, and immune response as well as DEGs such as CCL27, GRIA2, and SCN3A may be the key factors that mediate the process of response to chronic pain in goats. Neuroactive ligand-receptor interaction and NMDA receptor as well as GABA-related DEGs (GABRG3, GABRB2, and GABRB1), SCN9A, and TRPV1 may be involved in the process of response to pain in sheep. Our findings may provide some drug targets to prevent chronic pain.

Notes

Abbreviations

CFA: 

complete Freund’s adjuvant

DRG: 

dorsal root ganglion

DEG: 

differentially expressed gene

GO: 

gene ontology

Declarations

Authors’ contributions

HZ carried out the conception and design of the research. DW participated in obtaining funding. HW participated in the acquisition of data and revised the manuscript. SW carried out the analysis and interpretation of data. XD participated in the acquisition, analysis and interpretation of data, drafting the manuscript and revision of manuscript for important intellectual content. All authors read and approved the final manuscript.

Acknowledgements

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

Not applicable. This study was a primary research and further study is in progress.

Consent for publication

Not applicable.

Ethics approval and consent to participate

This study was approved by Ethics Committee of Inner Mongolia Agricultural University, Inner Mongolia Medical University and Inner Mongolia People’s Hospital Hohhot.

Funding

Not applicable.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
College of life Science, Inner Mongolia Agricultural University, No. 306 Zhaowuda Road, Saihan District, Hohhot, 010018, People’s Republic of China
(2)
College of Basic Medicine, Inner Mongolia Medical University, Hohhot, 010110, People’s Republic of China
(3)
Neurology Department, Inner Mongolia People’s Hospital, Hohhot, 010017, People’s Republic of China

References

  1. Gatchel RJ, Peng YB, Peters ML, Fuchs PN, Turk DC. The biopsychosocial approach to chronic pain: scientific advances and future directions. Psychol Bull. 2007;133(4):581.View ArticlePubMedGoogle Scholar
  2. Bushnell MC, Čeko M, Low LA. Cognitive and emotional control of pain and its disruption in chronic pain. Nat Rev Neurosci. 2013;14(7):502.View ArticlePubMedPubMed CentralGoogle Scholar
  3. Jackson JE. “Camp Pain”: talking with chronic pain patients. Philadelphia: University of Pennsylvania Press, Inc.; 2011.Google Scholar
  4. King CD, Keil A, Sibille KT. Chronic pain and perceived stress. 2016.Google Scholar
  5. Costa LD, Maher CG, McAuley JH, Hancock MJ, Herbert RD, Refshauge KM, Henschke N. Prognosis for patients with chronic low back pain: inception cohort study. BMJ. 2009;339:b3829.View ArticlePubMed CentralGoogle Scholar
  6. Mogil JS. Pain genetics: past, present and future. Trends Genet. 2012;28(6):258–66.View ArticlePubMedGoogle Scholar
  7. Krames ES. The dorsal root ganglion in chronic pain and as a target for neuromodulation: a review. Neuromodul J Int Neuromodul Soc. 2015;18(1):24.View ArticleGoogle Scholar
  8. Liem L, Van DE, Huygen FJ, Staats P, Kramer J. The dorsal root ganglion as a therapeutic target for chronic pain. Reg Anesth Pain Med. 2016;41(4):1.View ArticleGoogle Scholar
  9. Mccallum JB, Kwok WM, Sapunar D, Fuchs A, Hogan QH. Painful peripheral nerve injury decreases calcium current in axotomized sensory neurons. Anesthesiology. 2006;105(1):160.View ArticlePubMedPubMed CentralGoogle Scholar
  10. Krames ES. The dorsal root ganglion in chronic pain and as a target for neuromodulation: a review. Neuromodulation. 2015;18(1):24–32.View ArticlePubMedGoogle Scholar
  11. Pope JE, Deer TR, Kramer J. A systematic review: current and future directions of dorsal root ganglion therapeutics to treat chronic pain. Pain Med. 2013;14(10):1477.View ArticlePubMedGoogle Scholar
  12. Deer TR, Grigsby E, Weiner RL, Wilcosky B, Kramer JM. A prospective study of dorsal root ganglion stimulation for the relief of chronic pain. Neuromodul J Int Neuromodul Soc. 2013;16(1):67–72.View ArticleGoogle Scholar
  13. Wu Z, Li L, Xie F, Du J, Zuo Y, Frost JA, Carlton SPD, Walters ET, Yang Q. Activation of KCNQ channels suppresses spontaneous activity in DRG neurons and reduces chronic pain after spinal cord injury. J Neurotrauma. 2017;34(6):1260.View ArticlePubMedPubMed CentralGoogle Scholar
  14. Chen L, Liu Y, Yue K, Ru Q, Xiong Q, Ma B, Tian X, Li C. Differential expression of ATP-gated P2X receptors in DRG between chronic neuropathic pain and visceralgia rat models. Purinergic Signal. 2016;12(1):79–87.View ArticlePubMedGoogle Scholar
  15. Su L, Wang C, Yu YH, Ren YY, Xie KL, Wang GL. Role of TRPM8 in dorsal root ganglion in nerve injury-induced chronic pain. BMC Neurosci. 2011;12(1):120.View ArticlePubMedPubMed CentralGoogle Scholar
  16. Strong JA, Xie W, Coyle DE, Zhang JM. Microarray analysis of rat sensory ganglia after local inflammation implicates novel cytokines in pain. PLoS ONE. 2012;7(7):e40779.View ArticlePubMedPubMed CentralGoogle Scholar
  17. Bertorelli R, Corradini L, Rafiq K, Tupper J, Calò G, Ongini E. Nociceptin and the ORL-1 ligand [Phe1ψ (CH2-NH)Gly2]nociceptin(1–13)NH2 exert anti-opioid effects in the Freund’s adjuvant-induced arthritic rat model of chronic pain. Br J Pharmacol. 1999;128(6):1252–8.View ArticlePubMedPubMed CentralGoogle Scholar
  18. Hassenbusch SJ, Satterfield WC, Gradert TL. A sheep model for continuous intrathecal infusion of test substances. Hum Exp Toxicol. 1999;18(2):82–7.View ArticlePubMedGoogle Scholar
  19. Scheerlinck JP, Snibson KJ, Bowles VM, Sutton P. Biomedical applications of sheep models: from asthma to vaccines. Trends Biotechnol. 2008;26(5):259–66.View ArticlePubMedGoogle Scholar
  20. Atanasov AT. The linear allometric relationship between total metabolic energy per life span and body mass of mammals. Biosystems. 2007;90(1):224.View ArticlePubMedGoogle Scholar
  21. Ley SJ, Waterman AE, Livingston A. A field study of the effect of lameness on mechanical nociceptive thresholds in sheep. Vet Rec. 1995;137(4):85–7.View ArticlePubMedGoogle Scholar
  22. Dolan S, Evans NP, Richter TA, Nolan AM. Expression of gonadotropin-releasing hormone and gonadotropin-releasing hormone receptor in sheep spinal cord. Neurosci Lett. 2003;346(1–2):120–2.View ArticlePubMedGoogle Scholar
  23. Nagakura Y, Okada M, Kohara A, Kiso T, Toya T, Iwai A, Wanibuchi F, Yamaguchi T. Allodynia and hyperalgesia in adjuvant-induced arthritic rats: time course of progression and efficacy of analgesics. J Pharmacol Exp Ther. 2003;306(2):490–7.View ArticlePubMedGoogle Scholar
  24. Miki K, Zhou QQ, Guo W, Guan Y, Terayama R, Dubner R, Ren K. Changes in gene expression and neuronal phenotype in brain stem pain modulatory circuitry after inflammation. J Neurophysiol. 2002;87(2):750.View ArticlePubMedGoogle Scholar
  25. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.View ArticlePubMedPubMed CentralGoogle Scholar
  26. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):R25.View ArticlePubMedPubMed CentralGoogle Scholar
  27. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25–9.View ArticlePubMedPubMed CentralGoogle Scholar
  28. Da Wei Huang BTS, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2008;4(1):44–57.View ArticleGoogle Scholar
  29. Bennett GJ. Update on the neurophysiology of pain transmission and modulation: focus on the NMDA-receptor. J Pain Symptom Manage. 2000;19(1 Suppl):S2.View ArticlePubMedGoogle Scholar
  30. Qiu S, Chen T, Koga K, Guo YY, Xu H, Song Q, Wang JJ, Descalzi G, Kaang BK, Luo JH. An increase in synaptic NMDA receptors in the insular cortex contributes to neuropathic pain. Sci Signal. 2013;6(275):ra34.View ArticlePubMedGoogle Scholar
  31. Wang WT, Pan GQ, Zhang ZY, Suo ZW, Yang X, Hu XD. Ht31 peptide inhibited inflammatory pain by blocking NMDA receptor-mediated nociceptive transmission in spinal dorsal horn of mice. Neuropharmacology. 2015;89:290.View ArticlePubMedGoogle Scholar
  32. Zhang J, Echeverry S, Lim TK, Lee SH, Shi XQ, Huang H. Can modulating inflammatory response be a good strategy to treat neuropathic pain? Curr Pharm Des. 2015;21(7):831–9.View ArticlePubMedGoogle Scholar
  33. De MM, Kraychete DC, Meyer Nascimento RJ. Chronic pain: cytokines, lymphocytes and chemokines. Inflamm Allergy Drug Targets. 2014;13(5):339.Google Scholar
  34. Liu XJ, Salter MW. Glutamate receptor phosphorylation and trafficking in pain plasticity in spinal cord dorsal horn. Eur J Neurosci. 2010;32(2):278–89.View ArticlePubMedPubMed CentralGoogle Scholar
  35. Su S, Shao J, Zhao Q, Ren X, Cai W, Li L, Bai Q, Chen X, Xu B, Wang J. MiR-30b attenuates neuropathic pain by regulating voltage-gated sodium channel Nav1.3 in rats. Front Mol Neurosci. 2017;10:126.View ArticlePubMedPubMed CentralGoogle Scholar
  36. Hossaini M, Duraku LS, Saraç C, Jongen JL, Holstege JC. Differential distribution of activated spinal neurons containing glycine and/or GABA and expressing c-fos in acute and chronic pain models. Pain. 2010;151(2):356.View ArticlePubMedGoogle Scholar
  37. Zeilhofer HU, Benke D, Yevenes GE. Chronic pain states: pharmacological strategies to restore diminished inhibitory spinal pain control. Annu Rev Pharmacol. 2012;52(52):111–33.View ArticleGoogle Scholar
  38. Gwak YS, Hulsebosch CE. GABA and central neuropathic pain following spinal cord injury. Neuropharmacology. 2011;60(5):799–808.View ArticlePubMedPubMed CentralGoogle Scholar
  39. Holliday KL, Thomson W, Neogi T, Felson DT, Wang K, Wu FC, Huhtaniemi IT, Bartfai G, Casanueva F, Forti G. The non-synonymous SNP, R1150W, in SCN9A is not associated with chronic widespread pain susceptibility. Mol Pain. 2012;8(1):72.PubMedPubMed CentralGoogle Scholar
  40. Duan G, Xiang G, Zhang X, Yuan R, Zhan H, Qi D. A single-nucleotide polymorphism in SCN9A may decrease postoperative pain sensitivity in the general population. Anesthesiology. 2013;118(2):436–42.View ArticlePubMedGoogle Scholar
  41. Kim YS, Chu Y, Han L, Li M, Li Z, Lavinka PC, Sun S, Tang Z, Park K, Caterina MJ. Central terminal sensitization of TRPV1 by descending serotonergic facilitation modulates chronic pain. Neuron. 2014;81(4):873–87.View ArticlePubMedPubMed CentralGoogle Scholar
  42. Premkumar LS. Targeting TRPV1 as an alternative approach to narcotic analgesics to treat chronic pain conditions. AAPS J. 2010;12(3):361–70.View ArticlePubMedPubMed CentralGoogle Scholar

Copyright

Advertisement