SciELO - Scientific Electronic Library Online

 
vol.51Knockdown of PAICS inhibits malignant proliferation of human breast cancer cell linesBioinformatics analysis of microarray data to reveal the pathogenesis of diffuse intrinsic pontine glioma author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

  • On index processCited by Google
  • Have no similar articlesSimilars in SciELO
  • On index processSimilars in Google

Share


Biological Research

Print version ISSN 0716-9760

Biol. Res. vol.51  Santiago  2018  Epub Aug 29, 2018

http://dx.doi.org/10.1186/s40659-018-0174-7 

Research Article

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

Xiuling Deng1  2 

Dong Wang3 

Shenyuan Wang1 

Haisheng Wang2 

Huanmin Zhou1  * 

1College of life Science, Inner Mongolia Agricultural University, No. 306 Zhaowuda Road, Saihan District, Hohhot 010018, People's Republic of China

2College of Basic Medicine, Inner Mongolia Medical University, Hohhot 010110, People's Republic of China

3Neurology Department, Inner Mongolia People's Hospital, Hohhot 010017, People's Republic of China

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 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 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 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 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 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.

Additional file

Additional file 1: Table S1. GO term analysis of differentially expressed genes between strains. Table S2. GO term analysis of differentially expressed genes only in sheep or goat. Table S3. GO term analysis of codifferentially expressed genes in sheep and goat.

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.

Published online: 17 August 2018

Abbreviations

Abbreviation Expansion

CFA:

complete Freund's adjuvant;

DRG:

dorsal root ganglion;

DEG:

differentially expressed gene;

GO:

gene ontology.

Acknowledgements

Not applicable.

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. [ Links ]

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. [ Links ]

3. Jackson JE. “Camp Pain”: talking with chronic pain patients. Philadelphia: University of Pennsylvania Press, Inc.; 2011. [ Links ]

4. King CD, Keil A, Sibille KT. Chronic pain and perceived stress. 2016. [ Links ]

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. [ Links ]

6. Mogil JS. Pain genetics: past, present and future. Trends Genet. 2012;28(6):258-66. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

10. Krames ES. The dorsal root ganglion in chronic pain and as a target for neuromodulation: a review. Neuromodulation. 2015;18(1):24-32. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

20. Atanasov AT. The linear allometric relationship between total metabolic energy per life span and body mass of mammals. Biosystems. 2007;90(1):224. [ Links ]

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. [ Links ]

22. Dolan S, Evans NP, Richter TA, Nolan AM. Expression of gonadotropinreleasing hormone and gonadotropin-releasing hormone receptor in sheep spinal cord. Neurosci Lett. 2003;346(1-2):120-2. [ Links ]

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. [ Links ]

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. [ Links ]

25. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105-11. [ Links ]

26. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):R25. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

33. De MM, Kraychete DC, Meyer Nascimento RJ. Chronic pain: cytokines, lymphocytes and chemokines. Inflamm Allergy Drug Targets. 2014;13(5):339. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

38. Gwak YS, Hulsebosch CE. GABA and central neuropathic pain following spinal cord injury. Neuropharmacology. 2011;60(5):799-808. [ Links ]

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. [ Links ]

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. [ Links ]

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. [ Links ]

42. Premkumar LS. Targeting TRPV1 as an alternative approach to narcotic analgesics to treat chronic pain conditions. AAPS J. 2010;12(3):361-70. [ Links ]

Received: April 16, 2018; Accepted: August 07, 2018

* Correspondence: huanminzhou@263.net

Xiuling Deng and Dong Wang are co-first authors

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.

Competing interests

The authors declare that they have no competing interests.

Creative Commons License This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License, which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.