Gene Network Dysregulation in the Trigeminal Ganglia and Nucleus Accumbens of a Model of Chronic Migraine-Associated Hyperalgesia

The pharmacological agent nitroglycerin (NTG) elicits hyperalgesia and allodynia in mice. This model has been used to study the neurological disorder of trigeminovascular pain or migraine, a debilitating form of hyperalgesia. The present study validates hyperalgesia in an established mouse model of chronic migraine triggered by NTG and advances the understanding of the associated molecular mechanisms. The RNA-seq profiles of two nervous system regions associated with pain, the trigeminal ganglia (TG) and the nucleus accumbens (NAc), were compared in mice receiving chronic NTG treatment relative to control (CON) mice. Among the 109 genes that exhibited an NTG treatment-by-region interaction, solute carrier family 32 (GABA vesicular transporter) member 1 (Slc32a1) and preproenkephalin (Penk) exhibited reversal of expression patterns between the NTG and CON groups. Erb-b2 receptor tyrosine kinase 4 (Erbb4) and solute carrier family 1 (glial high affinity glutamate transporter) member 2 (Slc1a2) exhibited consistent differential expression between treatments across regions albeit at different magnitude. Period circadian clock 1 (Per1) was among the 165 genes that exhibited significant NTG treatment effect. Biological processes disrupted by NTG in a region-specific manner included adaptive and innate immune responses; whereas glutamatergic and dopaminergic synapses and rhythmic process were disrupted in both regions. Regulatory network reconstruction highlighted the widespread role of several transcription factors (including Snrnp70, Smad1, Pax6, Cebpa, and Smpx) among the NTG-disrupted target genes. These results advance the understanding of the molecular mechanisms of hyperalgesia that can be applied to therapies to ameliorate chronic pain and migraine.


INTRODUCTION
Migraine affects approximately 15% of the world's population and the World Health Organization considers migraine a top ten most disabling conditions with migraineurs experiencing enhanced intensity of pain or hyperalgesia. The therapeutic options for chronic migraine sufferers are limited and provide incomplete symptom relief (Bigal et al., 2008). The nitric oxide (NO) donor nitroglycerin (NTG) has been used extensively to understand the nociceptive system and pain processing (Cury et al., 2011). Nitroglycerin provokes pain in migraine-susceptible patients (Christiansen et al., 1999;Afridi et al., 2005;Olesen, 2008), and has been shown to evoke hyperalgesia in rodents (Di et al., 2016;Ferrari et al., 2016;Tipton et al., 2016;Demartini et al., 2017).
The trigeminal ganglia (TG) and nucleus accumbens (NAc) are two nervous system regions associated with chronic migraine and pain perception (Akerman et al., 2011;Schwartz et al., 2014). In response to NO, TG neurons trigger vasodilation and neurogenic inflammation, which elicit hyperalgesia, sensitization, and allodynia (Bellamy et al., 2006;Akerman et al., 2011). Chronic pain and migraine signal other nervous system regions including disruption of the nucleus accumbens (NAc) processes that play a role in migraine and hyperalgesia comorbidities including depression, irritability, fatigue, sleepiness, and loss of appetite (Burstein and Jakubowski, 2005;De Felice et al., 2013;Yuan et al., 2013;Elman and Borsook, 2016). Alterations in neural circuits in these nervous system regions have also been shown to play a role in medication overuse headache (Calabresi and Cupini, 2005;Torta et al., 2016) that can subsequently result in chronic migraine symptoms.
Studies of NTG-treated rodents can offer important insights into the mechanisms of hyperalgesia that can shed light into solutions for chronic pain such as migraine (Bates et al., 2010;Ferrari et al., 2015;Pedersen et al., 2016;Sufka et al., 2016). A study of the effect of one acute NTG infusion in the TG of male rats uncovered 15 differentially expressed genes relative to controls including glutamine synthetase (Glul), period circadian clock 1 (Per1), and genes that participate in immune responses including TAP binding protein (Tapbp), RT1 class Ia, locus A2 (RT1-A2) and RT1 class I, locus A3 (RT1-A3) (Pedersen et al., 2016). The molecular mechanisms disrupted by chronic NTG treatment evoking hyperalgesia in mice are incompletely understood. A comprehensive study of these pathways can offer insights into effective therapies to alleviate chronic pain and migraine (Bates et al., 2010;Ferrari et al., 2015;Pedersen et al., 2016;Sufka et al., 2016).
The objective of this study is to advance the understanding of the molecular disruptions that occur in a chronic hyperalgesic state associated with migraine. A transcriptomic analysis was undertaken to identify genes, biological processes and regulatory networks impacted by chronic NTG exposure in the TG and NAc of mice. The characterization of gene and pathway dysregulation can offer insights into the molecular mechanisms disrupted by NTG-evoked hyperalgesia in mice. These findings can shed light into solutions for chronic pain conditions such as migraine (Bates et al., 2010;Ferrari et al., 2015;Pedersen et al., 2016;Sufka et al., 2016).

Animals Experiments
Male C57BL6/J mice (Jackson Laboratories, Bar Harbor, ME) between 9 and 12 weeks old were studied. Mice were group housed in a 12-12 light-dark cycle, and food was available ad libitum. All experimental procedures were approved by the University of Illinois at Chicago Office of Animal Care and Institutional Biosafety Committee, in accordance with AALAC guidelines and the Animal Care Policies of the University of Illinois at Chicago, as well as with the European Union directive on the subject of animal rights. Mice were weighed daily during treatment, and no adverse effects of treatment were observed on body weight.
Chronic migraine-associated pain was modeled using intraperitoneal injections of NTG (10 mg/kg, IP) every second day for 9 days totaling 5 test days (Pradhan et al., 2014a;Tipton et al., 2016). The NTG injection was prepared daily from a stock solution of 5.0 mg/mL NTG and diluted in 0.9% saline. Control mice received intraperitoneal injections of 0.9% saline and were treated in parallel to NTG-treated mice. Mice were randomly allocated to either chronic NTG (NTG) or control (CON) groups. Male mice were studied to remove the effect of the estrus cycle stage on the molecular disruptions triggered by NTG.
The hind paw and cephalic mechanical threshold experiments used separate groups of mice from the same population. For all behavioral experiments, mice were counterbalanced into treatment groups following the first basal test for mechanical sensitivity on naïve mice. The experimenter was blinded to the drug condition being tested, basal responses are assessed before daily injection and injection volume was 10 ml/kg. No adverse effects of injection were observed in any of the experiments. All mice were tested in a separate behavior room with low-light (∼35-50 lux) and low-noise conditions, between 09:00 and 16:00. For all behavioral tests, mice were habituated to the testing rack for 2 days prior to the first test day, and on each test day for 20 minutes prior to the first measurement. For cephalic testing, mice were tested in 118 ml paper cups, to which they had been previously habituated for 1 h over 2 days.
The plantar surface of the mouse hind paw and the periorbital region caudal to the eyes and near the midline were tested for the peripheral and cephalic mechanical threshold measurements. To assess mechanical sensitivity, the threshold for responses to punctate mechanical stimuli (mechanical allodynia) was tested according to the up-and-down method (Chaplan et al., 1994). The region of interest was stimulated with a series of eight von Frey hair filaments (bending force ranging from 0.00g to 2g). A response was defined as a lifting, shaking, or licking of the hind paw or head, depending on the region tested. The first filament tested was 0.4g. In the absence of a response, a heavier filament (up) was tried, and in the presence of a response, a lighter filament (down) was tested. This pattern was followed for a maximum of four filaments following the first response. The peripheral threshold was measured on 9 consecutive days and the cephalic threshold was measured on days 1, 5, and 9. The effect of treatment was evaluated using a two way ANOVA including the effects of treatment, day and interaction and testing included the post hoc Holm-Sidak adjustment.

RNA Extraction and Sequencing
Mice were anesthetized with pentobarbital (somnosol), euthanized, and intracardially perfused with ice-cold saline approximately 24 h after the last injection of NTG or saline. Brains were extracted and the TG and NAc were rapidly dissected, snap-frozen, and stored at −80 • C. Total RNA was obtained from each collected nervous system region of an individual mouse following manufacturer's instructions. Steps included tissue homogenization with TRIzol (Invitrogen, Carlsbad, CA) and ceramic beads (MO BIO, Carlsbad, CA), and RNA isolation using the RNA-kit (Omega Biotek, Norcross, GA).
All 20 RNA samples representing 2 treatments (NTG and CON) and 2 nervous system regions had RNA integrity number values above 7.5 and were individually analyzed. Paired-end reads 100 nt in length were sequenced using the HiSeq 4000 (Illumina, San Diego, CA) platform. The RNA-seq datasets for this study are available in the National Center for Biotechnology Information Gene Expression Omnibus (GEO) database (identifier GSE110194). The average Phred quality score of the reads was assessed using FastQC (Andrews, 2010). The nucleotide quality score was >30 across all read positions. The read sequences were deemed of high quality and were not trimmed.

RNA Quantification and Differential Expression Analysis
The paired-end reads of each of the samples were individually mapped to the C57Bl/6J mouse genome [version GRCm38, downloaded on October 2016 from NCBI (Pruitt et al., 2007)] using Kallisto (Bray et al., 2016) with default settings. Gene expression transcript counts were imported into R (version 3.2) for analysis using tximport (release 3.5) (Soneson et al., 2015) and analyzed with edgeR using default settings (Robinson et al., 2010). Gene expression transcript counts were analyzed using a generalized linear model including the main factors of treatment group (NTG or CON) and nervous system region (TG or NAc) and the interaction between treatment and region. All genes with 5 or more reads per treatment-region group were analyzed to ensure gene expression in all groups tested for differential expression. The Benjamini-Hochberg false discovery rate (FDR) was used to adjust the differential expression P-value for multiple testing (Benjamini and Hochberg, 1995).

Functional Enrichment and Network Inference
Enrichment of Gene Ontology (GO) biological processes (BPs) and molecular functions (MFs) and KEGG pathways among the differentially expressed genes were evaluated using complementary methods (Reiner et al., 2003;Caetano-Anollés et al., 2015, 2016Nixon et al., 2015;Gonzalez-Pena et al., 2016a,b). Differentially expressed genes were analyzed with the Database for Annotation, Visualization and Integrated Discovery (DAVID; Version 6.8) using the Direct GO terms available in this repository (Huang et al., 2009). The Mus musculus genome was used as background for testing. Enrichment of each category was assessed using the Expression Analysis Systematic Explorer (EASE) score that was computed based on a one-tailed jackknifed Fisher hypergeometric exact test. The clustering of functional categories facilitated the interpretation of enriched terms. The statistical significance of each cluster of categories was assessed using an enrichment score computed as the geometric mean of the -log 10 EASE scores of categories within each cluster (Serão et al., 2013;Caetano-Anollés et al., 2015, 2016Gonzalez-Pena et al., 2016a).
Additional insights into the functional categories impacted by the treatment and nervous system region were gained using the Gene Set Enrichment Analysis (GSEA) approach implemented in the software package GSEA-P (version 2.0) (Subramanian et al., 2007). GSEA-P provides an enrichment score of functional categories in the Mus musculus Molecular Signature Database (MSigDB) that is calculated using maximum deviation of cumulative sum based on the gene-phenotype correlation (Subramanian et al., 2007).
Genes exhibiting a significant interaction between treatment and nervous system region, significant treatment main effect, or significant region main effect were queried as potential targets of enriched transcription factors. Significantly differentially expressed genes within each group were searched against the database of target genes ranked according to the associated transcription factor motifs using the iRegulon (Verfaillie et al., 2014) plugin in the Cytoscape environment (Shannon et al., 2003). A transcription factor normalized enrichment score was computed for each group where a normalized enrichment score > 3.0 corresponds to an approximate FDR between 3 and 9% (Verfaillie et al., 2014). Figure 1 depicts the profile of peripheral (left graph A) and cephalic (right graph B) mechanical threshold across the days when mice received via IP injection of NTG (black markers) or saline (white markers). The values depicted are the means and standard error of the mean (SEM) and the P-values include a Holm-Sidak adjustment for multiple comparison. The interaction between treatment and day, the main effect of treatment and the main effect of day were significant at P-value < 0.01. The contrasts between treatments within day were also significant at P-value < 0.001 (denoted with " * * * ") or P-value < 0.01 (denoted with " * * ").

Summary of RNA-Seq Measurements
Approximately 3 billion reads were generated across all 20 samples and approximately 70 million paired-end sequence readings were obtained per sample. No statistical difference in the number of mapped reads was detected between sample groups. The average percentage of the reads mapped to the mouse transcriptome was approximately 82% (±4%). The RNAseq reads were mapped to transcripts and after filtering for low counts (<5 reads) 22,071 genes were analyzed for the effects of the interaction between NTG treatment and nervous system region, and the main effects of NTG treatment and region studied.

Nitroglycerin Treatment-by-Nervous System Region Interaction Effect on Gene Expression
Twenty-five genes exhibited significant (FDR-adjusted P-value < 0.14 or P-value ≤ 2.2E-04) differential expression for the interaction between treatment and region effect ( Table 1). An expanded list of 109 genes including the previous 25 genes and others exhibiting interaction effect is provided in Table S1. The majority of these genes (96 genes) were differentially expressed between the treatment groups in the TG whereas the remaining 13 genes were significantly differentially expressed between the treatments in the NAc (Table S1A).
The second pattern has consistent relative expression between treatments albeit the significance level differed between the investigated regions (20 genes). Among the genes exhibiting consistent relative expression between treatments across nervous system regions albeit at different magnitude levels were: POU domain, class 3, transcription factor 3 (Pou3f3); erb-b2 receptor tyrosine kinase 4 (Erbb4); and solute carrier family 1 (glial high affinity glutamate transporter), member 2 (Slc1a2).
The final pattern was characterized by significant differential expression in one region with no differential expression in another region (3 genes). Genes exhibiting this profile included: natriuretic peptide type C (Nppc) prohormone, and two long noncoding RNA genes (Gm15738 and Gm39717).

Functional Analysis of the Genes Exhibiting Nitroglycerin Treatment-by-Nervous System Region Interaction Effects
Enrichment analysis revealed two clusters of enriched GO categories associated with positive regulation of transcription (enrichment score = 1.8 and 1.5, respectively) that were identified among the genes exhibiting significant treatment-bynervous system region interaction effects ( Table 2). Table S2 includes these clusters and additional enriched GO categories and KEGG pathways among genes exhibiting treatment-by-region interaction effect.
Complementary GSEA used information on differential expression between NTG and CON treatments within the studied regions. Tables 3, 4 present the most enriched (P-value < 0.005; normalized enrichment score > |1.7|) categories, selected to minimize redundancies, in the TG and NAc respectively, and the corresponding Tables S3, S4 present extended and complete lists of categories.
Within the TG (Table 3)

Main Effect of Nitroglycerin Treatment on Gene Expression
Differential gene expression in response to NTG treatment across studied nervous system regions (main effect of NTG) was investigated. The 11 genes exhibiting the highest NTG effect are listed in Table 5. An expanded list of 165 genes differentially expressed between the NTG and CON groups including the previous genes is provided in Table S5. The majority of these differentially expressed genes (86%) were over-expressed in NTG-treated relative to CON mice. Among the genes with the highest over-expression in NTG-treated compared to CON mice were: Aldehyde dehydrogenase 1 family member A1 (Aldh1a1); cytotoxic T lymphocyte-associated protein 2 beta (Ctla2b); solute carrier family 7 cationic amino acid transporter, y+ (Slc7a2); circadian rhythm period 3 gene (Per3); inositol 1,4,5-trisphosphate receptor type 2 (Itpr2); and non-coding RNA sequences including A930013F10Rik, FIGURE 1 | Peripheral (A) and cephalic (B) mechanical threshold indicator of hyperalgesia in the chronic nitroglycerin (NTG) or control (CON) mouse groups across day of treatment. Mechanical thresholds were determined before each injection using von Frey filaments. Two-way analysis of variance including the effects of treatment, day, and interaction; A. Hindpaw measurement (n = 5/group); treatment, day, interaction P-value < 0.01; *** within day Holm-Sidak testing P-value < 0.001. (B) Periorbital region caudal to the eyes measurement (n = 6/group); treatment, day, interaction P-value < 0.05; ** within day Holm-Sidak testing P-value < 0.01, ***P-value < 0.001. Results are represented as mean ± SEM.

Functional Analysis of the Genes Exhibiting Nitroglycerin Treatment Effects
The most significant DAVID clusters of functional categories enriched among the genes differentially expressed across studied nervous system regions of NTG-treated and CON mice are presented in Table 6. The extended list of DAVID clusters of enriched categories among the differentially expressed genes in NTG-treated relative to CON mice among the overand under-expressed genes in NTG-treated relative to CON mice is presented in the Table S6. Informative clusters of enriched categories include the terms: glutamatergic synapse and dopaminergic synapse (KEGG mmu04724 and mmu04728), the rhythmic process (GO:0048511), protein phosphorylation (GO:0018108), ATP binding (GO:0005524), and kinase activity (GO:0004672), and ion transport processes (GO:0006811).
Results from the GSEA functional analysis within genes over-and under-expressed in NTG-treated relative to CON mice offered complementary insights into the impact of chronic NTG treatment. Table 7 lists the most significantly enriched (P-value < 0.0005; normalized enrichment score > |1.8|) and relevant functional categories, selected to minimize redundancies and the corresponding Table S7 presents the extended and complete lists of categories uncovered by GSEA. Nucleoside binding was enriched among the genes under-expressed in NTGtreated mice across both nervous system regions ( Table 7). Conversely, anion transmembrane transport (GO:0098656) was significantly enriched among the genes over-expressed in NTGtreated relative to CON mice. The GSEA analysis was able to narrow down this transport to L-alpha-amino acid (GO:0015807) and L-glutamic acid is an important source of GABA.

Main Effect of Nervous System Region on Gene Expression
The focus of this study is understanding the disruption in transcriptome by the hyperalgesia-evoking NTG treatment in both or either TG and NA. The comparison of gene expression between nervous system regions is provided as molecular confirmation of the regions studied and for completeness. Differences in gene expression between regions were, per se, expected.
Supporting the different roles of both investigated regions, 361 genes were differentially expressed at P-value < 5.0 E-10 and log fold change > |4| ( Table S8) and 37% of these genes were over-expressed in the TG relative to NAc. Table 8 lists the most differentially expressed genes between both nervous system regions. An 8-fold over-expression of the migraine susceptibility gene potassium two pore domain channel subfamily K member 18 (Kcnk18) in TG relative to NA (P-value < 5.0 E-10) was detected. Table 9 presents the most significantly enriched DAVID clusters (enrichment score > |1.7|) of selected and relevant enriched categories, and the corresponding Table S9 presents the extended list of these categories. The most enriched categories within the top clusters include: ion transport channel activity (GO:0005216; enrichment score = 6); KEGG pathway neuroactive ligand-receptor interaction (mmu04080; enrichment score = 5.3); and KEGG pathway endocannabinoid signaling (mmu04723; enrichment score = 3.1).

Regulatory Networks
Regulatory network analysis aided in the identification of transcription factors that target many of the genes dysregulated by NTG treatment. Table 10 lists the enriched transcription factors corresponding to three complementary lists of differentially expressed genes (P-value < 0.005) between the NTG and CON groups. These lists represented differential expression within TG (338 genes), within NAc (103 genes), and across both regions (165 genes). Enriched transcription factors in the TG list included: CCAAT/enhancer binding protein alpha (Cebpa; normalized enrichment score = 5.2) and small muscle protein, X-linked (Smpx; normalized enrichment score = 4.4). Enriched transcription factors in the NAc list included: direct IAP-binding protein with low PI (Diablo; normalized 1 | Gene profiles exhibiting interaction between nitroglycerin treatment and nervous system region effect on expression.

Overall
Comparisons between interaction treatment-region levels b Gene symbol P-value FDR P-value a NTG(NAc)-NTG(TG)-CON(TG)-NTG(TG)-NTG(TG)-CON(TG)- enrichment score = 4.81) and yin and yang 1 (Yy1; normalized enrichment score = 4.34). Enriched transcription factors across regions included: E2F transcription factor 1 (E2f1; normalized enrichment score = 5.1), BarH-like homeobox 1 (Barx1; normalized enrichment score = 4.5), transcription factor GATA binding protein 3 (Gata3), and NK2 homeobox 2 (Nkx2-2). Overall, our study of regulatory networks impacted by NTG treatment confirmed our finding (Table 2) related to the important role of transcriptional regulation on the chronic NTG model studied. Figures 2, 3 depict the relationship between the transcription factors (octagons) most enriched (normalized enrichment score > 3.5 and < 40 genes to facilitate visualization) among genes (ovals) differentially expressed between NTG-treated and CON mice within the TG and NAc, respectively. Figure 4 depicts the relationship between the transcriptional factors most enriched (normalized enrichment score >3.5 and <40 genes) among the genes exhibiting a NTG treatment effect. The thresholds were used to facilitate the visualization of relationships between transcription factors and between transcription factors and their target genes. Another notable finding is that the transcription factors in the regulatory network disrupted by NTG in the TG (Figure 2) are less connected by common target genes than the regulatory network disrupted by NTG in the NAc (Figure 3) and across nervous system regions (Figure 4).

DISCUSSION
The findings from the pain threshold experiment confirmed that the chronic intermittent NTG treatment used in this study triggered mechanical hyperalgesia in mouse population studied. The NO donor NTG has been used extensively to understand the nociceptive system and pain processing (Cury et al., 2011). Nitroglycerin reliably triggers headache in normal subjects, and migraine without aura in migraine susceptible patients (Iversen et al., 1989;Christiansen et al., 1999;Afridi et al., 2005); and NTG-evoked migraine is a commonly used experimental paradigm in humans (Olesen, 2008(Olesen, , 2010. Nitroglycerin-evoked hyperalgesia in rodents has been developed as a model for sensory hypersensitivity associated with migraine (Bates et al., 2010;Markovics et al., 2012). Acute NTG treatment was previously shown to produce thermal and mechanical allodynia in mice that was reversed by the anti-migraine therapies sumatriptan (Bates et al., 2010;Markovics et al., 2012), and a calcitoningene-related peptide receptor antagonist (Capuano et al., 2014). In addition, in a transgenic mouse model of familial migraine, animals expressing a human migraine gene (casein kinase 1 delta) showed an even greater sensitivity to NTG-evoked hyperalgesia (Brennan et al., 2013). Further, NTG has also been shown to produce light-aversive behavior (Markovics et al., 2012), and increased meningeal blood flow in mice (Greco et al., 2011;Markovics et al., 2012); two hallmark characteristics of migraine.  Expanding upon the acute NTG model, we have developed a model of chronic migraine-associated pain. Figure 1, depicts the enhanced peripheral and cephalic hypersensitivity observed in the NTG mouse group. Chronic intermittent NTG produces a long-lasting and severe basal hyperalgesia, which is sensitive to migraine preventive treatment. This model has been extensively characterized pharmacologically and we have used it to identify novel therapeutic targets for migraine (Pradhan et al., 2014b;Tipton et al., 2016;Moye and Pradhan, 2017), as well as to identify mechanisms that can lead to migraine chronification (Ben Aissa et al., 2017).
Both, the NAc and TG have been associated with pain perception (Akerman et al., 2011;Schwartz et al., 2014). The present study of these regions in mice displaying hyperalgesia evoked by NTG treatment enabled to investigate molecular players that are either distinctly or similarly impacted by this treatment. Our results offer novel insights into the nervous system region-dependent and region-independent molecular mechanisms associated with NTG-evoked chronic hyperalgesia. These findings can help in the identification of solutions to ameliorated chronic pain such as migraine.

Genes and Processes Affected by Nitroglycerin Treatment in a Region-Specific Manner
Among the genes exhibiting significant treatment-by-nervous system region interaction effects, the enrichment of positive regulation of transcription processes ( Table 2, Table S2) is consistent with the established activation of the transcription factor NF-kB in migraine studies (Reuter et al., 2002;Greco et al., 2005). Likewise, the enrichment of T cell-mediated cytotoxicity in the TG (Tables 3, 4) was also reported in the TG of rats exposed to acute NTG treatment (Pedersen et al., 2016). The genes annotated to this functional category include Lrrc8a, Efnb, and Siglech listed in Table 1, Table S1.  Several of the genes exhibiting region-specific NTG effect have been associated with migraine in humans. The activity of the volume regulated anion channel protein coded by Lrrc8 has been associated with spreading through the brain focal depolarization of neurons and glial cells that occur in humans and animal models during migraines (Mongin, 2016). Susceptibility for migraine loci on Wscd1 (Table S1) were identified on a large meta-analysis in humans (Gormley et al., 2016). TIMP Metallopeptidase Inhibitor 4 (Timp4), an inhibitor of the matrix metalloproteinases, has been linked to the pathophysiology of migraine (Bernecker et al., 2011).
The unexpected enrichment of the categories sensory perception of light stimuli and retina process among the genes over-expressed in NTG-treated compared to CON mice was elucidated by the evaluation of the genes supporting this finding (Tables 3, 4). Genes supporting these enrichments included calcium channel, voltagedependent, alpha 2/delta subunit 4 (Cacna2d4), RAR-related orphan receptor beta (Rorb) and phosphodiesterase 6G, cGMP-specific, rod, gamma (Pde6g). These genes are all associated with pain and migraine (Schleithoff et al., 1999;D'Souza et al., 2008;Descalzi et al., 2017), in addition to 8 | Top 20 most differentially expressed genes (P-value < 1.0E-10) between the trigeminal ganglia (TG) and the nucleus accumbens (NAc).

Gene symbol
Gene name TG-NAc a

Gfra3
Glial cell line derived neurotrophic factor family receptor alpha 3 9.21

Trappc3l
Trafficking protein particle complex 3 like 8.83 Calca Calcitonin/calcitonin-related polypeptide, alpha 8.78 light stimuli, thus explaining the enrichment of the latter categories. The significant treatment-by-region effect on Slc32a1 (Table 1) and ErbB4 (Table S1) support the enrichment of GABA-related pathways. Slc32a1 is used as a marker for GABAergic neurons that modulate transmission of pain-related signal (Du et al., 2017). ErbB4 is mainly expressed in the GABAergic interneurons and is a receptor of neuregulin which is associated with neuropathic pain (Yau et al., 2003). Our findings offer insights into the molecular mechanism that are the target of effective migraine therapies such as anticonvulsants, topiramate, and valproic acid that disrupt the GABA pathway (Johannessen and Johannessen, 2003;Calabresi et al., 2007;Mulleners and Chronicle, 2008). In addition to the GABA pathway, Slc32a1 is mapped to the morphine addiction, nicotine addiction, and endocannabinoid signaling pathways. Likewise, Penk-derived opioid neuropeptides are known to influence chronic tension-type headaches (Langemark et al., 1995;Packard and Ham, 1997). These findings uncover intersectionality between the molecular mechanisms disrupted by NTG-treatment that evoke hyperalgesia and by drugs of abuse.
An interesting finding is that NTG treatment was associated both with lower expression of genes involved in innate immune response in the NAc and with higher expression of genes involved in adaptive immune response in the TG (Table 4). This association is in agreement with reports of the role of the trigeminovascular system, neuropeptides, and inflammatory cytokines in the pathophysiology of migraine (Farajdokht et al., 2017).
The glutamate receptor signaling pathway was enriched among the genes over-expressed in the NAc of NTG-treated mice. This result is in agreement with reports of association between migraine and increased glutamate transmission van Dongen et al., 2017). This enrichment is supported by the significant treatment-by-region interaction on Slc1a2 (Table S1), a gene that has nucleotide variants associated with migraine symptoms (García-Martín et al., 2014). The expression of Opalin, a gene that codes for a component of myelin, exhibited treatment-by-region interaction effect (Table 1). This profile corroborates the enrichment of myelin formation among genes implicated in migraine pain signaling (Eising et al., 2016).
The regulation of reactive oxygen species process was enriched among the genes under-expressed in the NAc of NTG-treated relative to CON mice ( Table 4). The under-expression of these genes including the oxygen sensor Egln2 (Table S1) may elevate the level of reactive oxygen species that have been associated with migraine (Quaegebeur et al., 2016). An alternative mechanism that may trigger hyperalgesia in NTG-treated mice may be associated with the hormone Nppc that exhibited a significant treatment-by-region interaction effect (Table S1). Nppc can produce natriuretic peptides and has been associated with hemiplegic migraine (Marchenkova et al., 2016a,b).

Genes and Processes Affected by Nitroglycerin Treatment in a Region-Independent Manner
The enrichment of dopaminergic synapses pathway among genes differentially expressed between NTG and CON mice consistently across regions ( Table 6) agrees with indications that disrupted levels of the neurotransmitters such as dopamine and noradrenaline in the corresponding synapse clefts of the pain matrix could evoke in the trigeminal system the release of calcitonin gene-related peptide, a peptide involved in the transmission of migraine pain (D'Andrea and Leon, 2010). This process can be followed by build-up of inflammatory mediators and sensitization of the TG leading to migraine (D'Andrea and Leon, 2010). Dopamine transporters has been associated with migraine and accompanying drug abuse (Cevoli et al., 2006), further reinforcing the link between NTG and drug abusetriggered hyperalgesia and migraine. Aldh1a1 and Itpr2, both associated with the dopaminergic pathway, were over-expressed between NTG and CON mice ( Table 5). Over-expression of IItpr2 has been associated with induced visceral inflammatory hypersensitivity (Qian et al., 2013) and with neuropathic pain in rats (Maratou et al., 2009). Also over-expressed in NTG, the dysregulation of Ctla2b (Table 5) supports the proposal to use cathepsin inhibitors as antihyperalgesics (Irie et al., 2008). Our findings of under-expressed serotonin receptor and phosphorylation and kinase enrichment in NTG relative to CON mice ( Table 6) can be related to reports that serotonin depletion in rats increases nociception-evoked and kinase enabled phosphorylation in the TG (Maneepak et al., 2009). Enrichment of ion transport and phosphorylation processes among genes differentially expressed between NTG-treated and CON mice have been reported in studies of hemiplegic migraine (Ophoff et al., 1996;De Fusco et al., 2003;Dichgans et al., 2005;Schack et al., 2012;Pietrobon, 2013).
The enrichment of the nucleoside binding function among the genes under-expressed in NTG-treated relative to CON mice ( Table 7) supports the use of modulators of G protein-coupled receptors that have nucleoside binding capabilities to treat pain and migraine (Venkatakrishnan et al., 2013). The enrichment of the amino acid transport process ( Table 7) is supported by the dysregulation of Slc7a2 ( Table 5). The arginine transport facilitated by Slc7a2 can be associated with migraine because arginine is a chemical precursor to NO. Likewise, Aldh1a1 was dysregulated by NTG and aldehyde dehydrogenases catalyze the conversion of NTG into NO (Sydow et al., 2004). Nitric oxide acts as a vasodilator and triggers headaches whereas inhibition of NO synthase is an effective treatment to relieve migraine attacks (Olesen, 2008). The 8-fold over-expression of the migraine susceptibility gene Kcnk18 in TG relative to NA (P-value < 5.0 E-10) is consistent with reports of over-expression of this gene in the neural-enriched human trigeminal samples (LaPaglia et al., 2018).
The enrichment of rhythmic processes among genes dysregulated by NTG treatment ( Table 6) is supported by the over-expression of Per3 in NTG relative to CON mice ( Table 5). Per3 has been associated with cluster headache (Ofte et al., 2013) and dysregulation of circadian rhythm has been associated with migraine (Solomon, 1991;Hering and Kuritzky, 1992;Pringsheim, 2002). Non-coding RNA sequences including A930013F10Rik and 1700024F13Rik (Table 5) and 4930572G02Rik and E030003E18Rik (Table S5) were overexpressed in NTG-treated relative to CON mice. Non-coding RNAs that regulate glutamate transporters have been proposed for migraine therapy (Gasparini et al., 2016).

Regulatory Networks
Regulatory network analysis of gene expression disruptions triggered by NTG offers insights into the role of transcription factors and their target genes on chronic hyperalgesia (Table 10). The transcription factors Cebpa and Smpx were enriched among the genes dysregulated in the TG of NTG-treated mice (Table 10). Cebpa has been associated with attenuation of pain hypersensitivity in the spinal cord of mice (Jiang et al., 2017). In mice, Smpx is dysregulated by disruption of neurolysin, a molecule implicated in two processes associated with hyperalgesia: pain control and blood pressure regulation (Cavalcanti et al., 2014). The lower level of interconnection among transcription factors and targets genes observed in the TG network (Figure 2) suggests that many therapeutic targets may be needed to effectively revert the network disruption triggered by NTG and associated with chronic hyperalgesia.
A notable finding in the TG regulatory network disrupted by NTG (Figure 2) is that all four enriched transcription factors (Snrnp70, Smad1, Pax6, and Cebpa) are related through other transcription factors. For example, Olig1 can have an inhibitory effect on GABAergic interneurons that attenuate neuropathic pain in rodents (Silbereis et al., 2014). Also Nkx2-2 is expressed in cells that give rise to serotonergic neurons and the modulation of these neurons has been implicated in the sensation of pain (Pattyn et al., 2004).
The transcription factors Diablo, Six6, Rfx2, and Yy1 were enriched among the genes dysregulated by NTG in the NAc (Table 10, Figure 3). Our results support reports that Diablo is elevated in rats treated with alkaloid compounds that have FIGURE 3 | Relationship between the transcription factors (octagons) most enriched (normalized enrichment score >3.5 and < 40 genes to facilitate visualization) among genes (ovals) differentially expressed between nitroglycerin-treated and control mice within the nucleus accumbens.
been linked to migraine events associated with ischemic stroke (Alsharafi et al., 2015). Yy1 was enriched among inflammatory pathway genes dysregulated in a transgenic mouse model of migraine (Eising et al., 2017). Six6 is associated with a hereditary disorder of the eye that presents with migraine symptoms (Burdon et al., 2015) and decreased expression of Rfx2 elicited increased sensitivity in a model of neuropathic pain (Wheeler et al., 2013). Six6 is a transcription factor for H2A histone family member Z (H2afz) that has been associated with the effects of abused substances (Vadasz et al., 2007;McBride et al., 2009) and migraine is one of these symptoms (Granella et al., 1987).
Disruption of NF-κb levels has been associated with migraine (Karatas et al., 2013) and acute nitroglycerin treatment in rats (Greco et al., 2005). In our study, transcription factor targets of NF-κb including Sox9 and Ahctf1 were enriched or differentially expressed whereas only NF-κb p100 subunit (Nfkb2) was over-expressed in NTG-treated relative to CON mice. This observation suggests that NF-κb may initiate pain whereas downstream molecular targets may maintain the longlasting molecular dysregulation responsible for chronic pain sensation.

CONCLUSIONS
Comparison of the NTG and CON groups confirmed that chronic NTG treatment substantially decreased the mechanical pain threshold in the treated mice. The innovative transcriptome comparison between two central nervous regions offered novel insights into region-dependent and -independent differential gene expression and networks that can be associated with hyperalgesia. The genes that exhibited opposite differential expression in response to NTG across regions included: Opalin, Efnb3, Slc32a1, Penk, Rorb, Lrrc8a, Cacna1b, H2-Eb2, Gm1673, and Egln2. The genes that exhibited consistent NTG expression pattern, albeit of different magnitudes, across regions included FIGURE 4 | Relationship between the transcription factors (octagons) most enriched (normalized enrichment score >3.5 and <40 genes to facilitate visualization) among genes (ovals) differentially expressed between nitroglycerin-treated and control mice.
The detection of differential abundance in multiple genes reported to be distinctively expressed within the regions studied supports the capacity of the experimental design used to detect molecular differences. Further validation of the differential expression patterns using additional quantitative technologies was prevented by insufficient sample from the tested mice and regions. Partially addressing this situation, our discussion focused on differentially expressed genes that were both supported by enriched functional categories and have been previously associated with allodynia and migraine phenotypes with preference to studies in rodents, NTG treatment and similar regions.
Hyperalgesia and migraine phenotypes, and therefore the underlying molecular mechanisms, are not static but rather tend to be condition-dependent. Migraine affects more females than males and aged than young subjects. Moreover, migraine, depression, and stress tend to be present as comorbidies. The present study focused on comparing the transcriptome of two CNS regions that can participate on these phenotypes in adult male mice. Follow-up studies including the comparison of phenotypes and transcriptome patterns associated with hyperalgesia in females and males, young and aged, and across behavior groups are expected to enhance the understanding of molecular players associated with hyperalgesia and migraine. Likewise, additional studies including migraine interventional therapies (sumatriptan, topiramate, and propranolol) and longitudinal studies including NTG treatment suspension and resumption will offer insights into potential multi-factorial interactions that affect the gene expression profiles associated with hyperalgesia and migraine.
The detected enrichment of positive regulation of transcription processes among the genes exhibiting significant treatment-by-nervous system region interaction effects is consistent with the established activation of the transcription factor NF-kB in migraine studies. Similarly, the enrichment of T cell-mediated cytotoxicity was also reported in the TG of rats exposed to acute NTG treatment. Several of the genes exhibiting region-specific NTG effect have been associated with migraine in humans including Lrrc8 and Wscd1.
The significant enrichment of the GABA-related pathway among genes presenting a treatment-by-region effect offer insights into the molecular mechanism that are the target of effective migraine therapies such as anticonvulsants, topiramate, and valproic acid that disrupt the GABA pathway. The enrichment of the GABA pathway and the differential expression of genes associated with addiction between treatments suggest intersectionality between the molecular mechanisms disrupted by NTG-treatment that evoke hyperalgesia and by drugs of abuse.
A thought-provoking result is simultaneous association of NTG treatment with under-expression of genes involved in innate immune response in the NAc and with over-expression of genes involved in adaptive immune response in the TG. This result confirms the participation of the trigeminovascular system, neuropeptides and inflammatory cytokines in the molecular processes of migraine. The known association between migraine and increased glutamate transmission glutamate receptor signaling pathway was confirmed in our study among the genes over-expressed in the NAc of NTG-treated mice. The enrichment of the process of regulation of reactive oxygen species process among the genes under-expressed in the NAc of NTG-treated relative to CON mice is consistent with reports of correspondence between elevated level of reactive oxygen species and migraine Regulatory network reconstruction highlighted the widespread role of several transcription factors (including Snrnp70, Smad1, Pax6, Cebpa, and Smpx) among the NTGdisrupted target genes. Overall, this work provides promising leads to help identify critical mechanisms associated with NTG-evoked hyperalgesia with therapeutic applications toward the modulation and management of chronic pain and migraine.