Evidence for a non-canonical JAK/STAT signaling pathway in the synthesis of the brain’s major ion channels and neurotransmitter receptors

Background Brain-derived neurotrophic factor (BDNF) is a major signaling molecule that the brain uses to control a vast network of intracellular cascades fundamental to properties of learning and memory, and cognition. While much is known about BDNF signaling in the healthy nervous system where it controls the mitogen activated protein kinase (MAPK) and cyclic-AMP pathways, less is known about its role in multiple brain disorders where it contributes to the dysregulated neuroplasticity seen in epilepsy and traumatic brain injury (TBI). We previously found that neurons respond to prolonged BDNF exposure (both in vivo (in models of epilepsy and TBI) and in vitro (in BDNF treated primary neuronal cultures)) by activating the Janus Kinase/Signal Transducer and Activator of Transcription (JAK/STAT) signaling pathway. This pathway is best known for its association with inflammatory cytokines in non-neuronal cells. Results Here, using deep RNA-sequencing of neurons exposed to BDNF in the presence and absence of well characterized JAK/STAT inhibitors, and without non-neuronal cells, we determine the BDNF transcriptome that is specifically regulated by agents that inhibit JAK/STAT signaling. Surprisingly, the BDNF-induced JAK/STAT transcriptome contains ion channels and neurotransmitter receptors coming from all the major classes expressed in the brain, along with key modulators of synaptic plasticity, neurogenesis, and axonal remodeling. Analysis of this dataset has revealed a unique non-canonical mechanism of JAK/STATs in neurons as differential gene expression mediated by STAT3 is not solely dependent upon phosphorylation at residue 705 and may involve a BDNF-induced interaction of STAT3 with Heterochromatin Protein 1 alpha (HP1α). Conclusions These findings suggest that the neuronal BDNF-induced JAK/STAT pathway involves more than STAT3 phosphorylation at 705, providing the first evidence for a non-canonical mechanism that may involve HP1α. Our analysis reveals that JAK/STAT signaling regulates many of the genes associated with epilepsy syndromes where BDNF levels are markedly elevated. Uncovering the mechanism of this novel form of BDNF signaling in the brain may provide a new direction for epilepsy therapeutics and open a window into the complex mechanisms of STAT3 transcriptional regulation in neurological disease. Electronic supplementary material The online version of this article (10.1186/s12864-019-6033-2) contains supplementary material, which is available to authorized users.


Background
In the injured brain, it has been suggested that neurogenesis [1], mossy fiber sprouting [2], and hippocampal cell death [3] contribute to aberrant neuronal circuit reorganization. One manifestation of this reorganization is thought to be altered long-term potentiation (LTP)-induced synaptic plasticity [4], which is associated with altered levels of particular neurotransmitter receptors and ion channels [5,6]. As the functional properties of neurons change, a state of unchecked overexcitation may occur and is accompanied by decreased GABA-mediated inhibition [7]. Many believe that a major mechanism for this altered brain plasticity is the increased synthesis and release of Brain-derived Neurotrophic Factor (BDNF), a major brain signaling molecule that plays an important role throughout life. BDNF is a neurotrophin ubiquitously expressed in the whole brain with a major role in activity-dependent alterations in neuronal morphology and synaptogenesis. Protein and mRNA levels of BDNF are markedly increased after seizures in brain regions involved in epileptogenesis [8,9] and mutations in BDNF can reduce the development of spontaneous seizures in animal models [10]. Moreover, a conditional deletion of the high-affinity BDNF receptor TrkB in a subset of neurons, is sufficient to completely eliminate all behavioral evidence of the progression of epilepsy in the kindling model [11] and it has been shown that enhanced TrkB signaling can exacerbate epileptogenesis [12].
In addition to the TrkB-mediated BDNF response, the p75 Neurotrophin Receptor (p75NTR, also referred to as NGFR) is a low-affinity receptor for the cleaved mature form of BDNF while a high-affinity receptor for the uncleaved form, proBDNF. p75NTR levels are elevated after Status Epilepticus (SE) in animals and increased activation by proBDNF, like mature BDNF, increases susceptibility to seizures [13,14]. Interestingly, our labs have shown that increased synthesis of proBDNF, and not mature BDNF, is the first response of the injured brain in the pilocarpine (PILO) model of epilepsy prior to SE-induced increases of p75NTR [15]. This finding suggests that high levels of proBDNF may exert their effects through TrkB rather than p75NTR.
Studies from our laboratories demonstrate that in addition to the signaling pathways described above, increased levels of BDNF activate the Janus kinase/signal transducer and activator of transcription (JAK/STAT) pathway both in vivo, in the rat PILO model of epilepsy, and in vitro, in BDNF-treated primary cultured neurons [23]. The JAK/STAT pathway is a signaling cascade that has a prominent role in immune function and cancer development. The canonical pathway involves the binding of JAK to its target (such as a cytokine or hormone receptor), subsequent JAK phosphorylation that stimulates the recruitment and activation of STAT proteins (phospho-STAT (p-STAT)) and their movement as dimers into the nucleus where they function as transcriptional activators. Although JAK/STATs regulate transcription through this canonical pathway, more recently, non-canonical JAK/STAT signaling has been discovered where JAKs can act independently of STATs to regulate transcription in the nucleus, and STATs can function independently of their phosphorylation state to maintain genome stability via interactions with heterochromatin protein 1 (HP1) [24]. Increasing evidence also suggests that HP1 may be present at euchromatin as well as heterochromatin and involved in gene transcription [25,26]. While the role, and mechanism, of JAK/STAT signaling in the brain is still not completely understood, it is clear that it is an important part of how the brain regulates its synaptic connections [27].
Research from our laboratory has demonstrated that activation of the JAK/STAT pathway results in the expression of the inducible cAMP early repressor (ICER) via STAT3-mediated gene regulation. Moreover, we show that ICER represses the gene (Gabra1/GABRA1) coding for the α1 subunit of the type A γ-aminobutyric acid (GABA) receptor (GABAR), the major inhibitory neurotransmitter receptor in the brain, and that such repression reduces the number of α1 containing GABARs in neurons [28]. Most importantly, reduced α1 subunit gene expression also occurs in epilepsy patients [5] as well as in other disorders of the nervous system [29][30][31]. In addition, use of the JAK/STAT inhibitor, WP1066, in the PILO rat model of epilepsy reduces the number of spontaneous seizures after the latent period [32]. We now hypothesize that in addition to controlling the gene regulation of Gabra1/GABRA1, the BDNF-induced JAK/ STAT pathway controls the expression of diverse gene products that are involved in synaptic plasticity as well as in the epileptogenic process that follows brain injury. We also hypothesize that attenuating this pathway may provide a gateway to new treatments for refractory epilepsy and other brain disorders that display dysregulated neuroplasticity.
To test these hypotheses, and to gain a comprehensive understanding of the genome impacted by BDNFinduced JAK/STAT activation, we exposed primary cultured neurons to BDNF with and without JAK/STAT inhibitors and performed deep RNA-sequencing (RNAseq) in order to determine the gene set within the BDNF response that could be ascribed to JAK/STAT activation (as defined by WP1066 and Ruxolitinib (Ruxo), a potent inhibitor of JAK1/2 activation [33]). We now report evidence that suggests BDNF acts through the JAK/STAT pathway to execute changes in: the expression of ion channels, neurotransmitter and GPCR receptors, and the expression of key modulators of synaptic plasticity, neurogenesis, and axonal remodeling. Many of these functions were previously ascribed to BDNF, however, the role of the JAK/STAT pathway in their manifestation has never been previously explored. Our studies also reveal the presence of many new BDNF target genes providing a broader landscape for BDNF mediated gene regulation of the healthy and impaired nervous system.

BDNF signaling and the identification of epilepsyassociated gene networks
Primary cultured neurons were treated with BDNF (0.7 nM, a physiologically relevant level) or vehicle control (water) to determine the genome response to prolonged receptor activation (4 h) (Fig. 1). Using RNA-sequencing for open discovery, we identified a total of 2869 differentially expressed genes (DEGs) whose levels change (increase or decrease) in response to BDNF signaling ( Fig. 2a) (CTRL vs V + B gene list, Additional file 5). Inspection of the gene set using Ingenuity Pathway Analysis (IPA) [34] revealed that only 83 of these genes have been previously associated with BDNF (as determined using the Ingenuity Knowledge Base, a curated database of all historical scientific information in IPA databases), suggesting that the majority of genes identified in this study are novel in their association with the BDNF signaling pathway. A subset of 194 potential BDNF regulated genes are associated with "seizure disorders" (BDNF and epilepsy, Additional file 5), an umbrella A C B Fig. 1 Schematic of experimental design and analysis. a Timeline representation of the drug treatment protocol, including the 1 h drug or DMSO (vehicle, V) pretreatment (initiated at − 1 h), administration of aqueous BDNF (B) or water (vehicle) treatment (Time = 0 h) and time of cell collection and RNA extraction (4 h). b 6-well plate showing each one of the 5 treatment groups and the abbreviations by which they will be referred to in the text: DMSO+Water (CTRL), DMSO+BDNF (V + B), 100 nM Ruxolitinib+BDNF (RX1 + B), 10 μM Ruxolitinib+BDNF (RX2 + B), 10 μM WP1066 + BDNF (WP + B). c Diagram represents the approach taken to identify differential gene expression (DEG) in response to JAK/STAT pathway inhibition. BDNF through its receptors activates multiple signaling pathways (represented by blue arrows) that impact transcription. Differential expression analysis comparing V + B vs. CTRL treatment groups (green box) will reveal the total set of genes that are regulated by BDNF-induced intracellular signaling pathways. To identify BDNF DEGs that are specific to JAK/STAT signaling, comparisons are made between the groups pretreated with JAK inhibitors (RX1 + B, RX2 + B, WP + B) and the group pretreated with vehicle (V + B) (purple box); DEGs from this comparison (where JAK/STAT inhibition reverses BDNF stimulation or inhibition) are thought to be associated with the JAK/STAT pathway category in IPA that includes all epilepsy classes as represented in the curated knowledge base. Moreover, seizure disorder is the primary most significant neurological disorder identified by enrichment analysis of the gene pool (p = 3.29E-27).
In addition to the 194 epilepsy-associated genes, pathway analysis suggests that increased levels of BDNF alters the expression of essential neurobiological gene sets that are required for brain development and that are dysregulated during epileptogenesis in animal models and human patients. For instance, the differentially expressed set of transcripts from BDNF-treated primary neurons is highly enriched for genes essential to axonal guidance (p = 1.09E-11); further, there is a striking enrichment that is specific to all major classes of ion channels and neurotransmitter/neuropeptide receptors, with 133 of 599 total genes related to one of these categories regulated by prolonged BDNF exposure (4 h) (p = 2.66E-11) (Fig. 2b). Specifically, the differentially expressed transcripts for receptors include those activated by GABA, glutamate, acetylcholine, dopamine, opioids, serotonin, galanin, and neuropeptide Y, covering both metabotropic and ionotropic receptor subtypes. Most classes of ion channel genes are also represented, including those for potassium, sodium, and calcium, where mutations have been shown to underlie specific genetic epilepsies [35][36][37]. There is also a specific decrease (1.73 fold) in the expression of Gabra1 that is a member of the GABAR gene cluster which we have previously demonstrated to be regulated by BDNF using qRT-PCR in extracts of BDNF-treated primary neurons and hippocampal tissue extracted from the in vivo PILO model of temporal lobe epilepsy [23]. As expected, the neuronal RNA-seq dataset revealed that prolonged exposure to BDNF increases BDNF transcript levels, a finding that is well established in the literature [38]. Given our previous identification of a novel JAK/STAT pathway in neurons that is regulated by BDNF, we asked whether any members of the JAK/STAT pathway might be regulated at a transcriptional level in response to BDNF receptor activation. Our dataset indicates significant enrichment for transcripts linked to JAK/STAT signaling (p = 6.03E-5) and in particular, an enrichment of Jak1 (1.6-fold enrichment) and Jak2 (1.8-fold enrichment).

Neuropharmacological dissection of BDNF-induced JAK/ STAT signaling reveals a non-canonical pathway in neurons
Primary cultured neurons were pre-treated with agents known to inhibit the JAK/STAT pathway ( Fig. 1) by interfering with JAK kinase activity, either directly by inhibiting JAK phosphorylation or stimulating JAK degradation. Two different antagonists were employed in these studies: the small molecule WP1066 that has been shown to inhibit JAK2/STAT3 signaling and degrade total JAK2 [39] and Ruxolitinib (Ruxo), currently in use as a cancer therapeutic, which acts on JAK1 and JAK2 by forming two hydrogen bonds with the hinge region of the JAK protein (the segment that connects the N-lobe to the C-lobe of the kinase domain) and reduces phosphorylation of the JAK2 protein at Tyr1007/1008 [33]. Based on our earlier discovery that siRNA knockdown of STAT3 reverses BDNF-induced downregulation in the levels of Gabra1 subunit mRNAs [23], as does WP1066, we used RNA-seq to identify the full set of genes responsive to JAK inhibition. This strategy enables us to unmask those genes in the BDNF transcriptome that may be JAK/STAT dependent.
Comparing the transcriptomes BDNF (V + B) or 100 nM Ruxo + BDNF (RX1 + B), we were surprised to find that there was no difference between these two transcriptomes, as statistically assessed using DSEQ2 (0 DEGs) (FDR = 0.05, Wald test). To make sure that Ruxo treatment at 100 nM was active in the primary cultures of our studies, Western blot was performed using the same drug stock as used in cultures that generated the RNA-seq libraries. In this control experiment, Ruxo (100 nM) potently inhibited STAT3 phosphorylation at A C B Fig. 3 RX1 and RX2 reduce levels of phospho-STAT3 but only RX2 blocks ICER induction. a Representative Western Blot analysis of whole-cell protein extracts from primary neurons 9 DIV pretreated with100nM Ruxolitinib (RX1), 10uM Ruxo (RX2), or DMSO for 1 h before the addition of BDNF. Cells were collected 30 min after BDNF administration and probed with anti pSTAT3, STAT3 and β-Actin. b Quantitation of signals from densitometry is displayed as mean percent change (±SEM) relative to the DMSO+Water control group (**p < 0.01). Ruxolitinib significantly reduces the levels of P-STAT3 with no change in total STAT3. c Graphical representation of results from real-time PCR analysis using specific Taqman probe and primers. RNA was extracted from cells collected 4 h after BDNF administration. Transcript levels are shown as the mean values (±SEM) of the ratio relative to the DMSO+Water control group for Icer n = 6. Note, that while RX1 blocks pSTAT3, it does not block the effect of BDNF treatment on ICER induction, while WP and RX2 do 30 min (80% + − 3.07 reduction of BDNF signal) and at 4 h (73% reduction of BDNF signal, Fig. 3). Therefore, despite reducing STAT3 phosphorylation at Tyr 705, our data taken together demonstrate that 100 nM Ruxo does not alter BDNF-induced gene regulation, including changes in ICER expression (Fig. 3c). This finding suggests that the BDNF-induced JAK/STAT pathway is not STAT3 phosphorylation-dependent even though knockdown of STAT3 with siRNAs prevents ICER induction and Gabra1 downregulation [23]. Note that in our current studies BDNF does not significantly increase STAT3 phosphorylation but Ruxo still inhibits the basal state of its activation. These observations favor the hypothesis that neurons use a non-canonical mechanism of JAK/STAT signaling which may be relevant to epilepsy, and potentially to learning and memory where BDNF is a critical signaling molecule, but require further experimentation.
Using pharmacology to reduce the dimensionality of the dataset and identify potential BDNF-induced JAK/STAT target genes While we discovered that co-treatment of neurons with BDNF (V + B) vs. RX1 + B produced similar transcriptomes, we found that a marked increase in the concentration of Ruxo (10 μM) plus BDNF (RX2 + B) attenuated the expression of multiple genes that were altered in response to BDNF, and most importantly, those that we previously identified in the PILO model of epilepsy. We next compared the RX2 DEG as determined by DESEQ2 (FDR = 0.05, Wald test) to those generated from cultures treated with 10 μM WP1066 plus BDNF (WP + B) and generated a subset of genes whose BDNFinduced expression was reversed back to vehicle control by both drug treatments. Using these JAK/STAT inhibitors as probes for unmasking potential genes that are regulated by the BDNF-induced JAK/STAT pathway, we identified 1559 of the 2659 BDNF target genes that may be most associated with JAK/STAT signaling in neurons. Interrogation of IPA networks specific to this shared gene set revealed that 131 genes are epilepsy-related and seizure disorder is still the top most significant neurological disease (P = 7.31E-22). Considering that 131 of the 194 BDNF-regulated epilepsy-linked genes are reversed by JAK/STAT inhibitors (68%), it suggests that JAK/STAT signaling may be a large component of BDNF-related epileptogenesis. This can be best appreciated by looking at the heat map displayed in Fig. 2c which depicts the expression levels of all 194 epilepsylinked genes that are changed by BDNF treatment. The results clearly show that WP + B (10 μM) and Ruxo (RX2 + B, 10 μM) block BDNF-induced changes in neuronal gene expression for epilepsy-linked genes and both appear to have a similar gene expression signature to that of DMSO+Water control. Note that Ruxo at its lower concentration (RX1 + B, 100 nM) looks similar to BDNF alone.

Gene ontology and pathway analysis
Gene Ontology (GO) term analysis of the WP + B/RX2 + B dataset, performed using EnrichR [40,41], reveals that the 1559 overlapping genes are significantly associated with molecular functions related to receptor and channel activity (found in 8/10 top molecular functions) (Fig. 4a,  top panel). All of the top 10 biological process GO terms include the words "synaptic transmission" (Fig. 4a, middle panel). These consistent findings strongly suggest that both WP and Ruxo have the capacity to regulate synaptic transmission at the gene expression level in the brain. When compared with the total list of receptor/ion channel related genes that were regulated by BDNF (as depicted in Fig. 2b, top panel), 107/133 genes (80%) were reversed by treatment with JAK/STAT inhibitors (Fig. 2b, bottom panel). This is best seen by the color reversal in the receptors for all but a few that were not significantly affected by pre-treatment with the inhibitor prior to BDNF activation (displayed in white).
In addition, neurotransmitter receptors, axonal guidance, GPCR and c-AMP signaling were found to be amongst the top KEGG pathways associated with the shared gene set of WP + B and RX2 + B (Fig. 4a, bottom  panel). Members of the JAK/STAT signaling pathway were also identified as significant gene targets (P = 1.044E-3) suggesting that BDNF may also participate in the feedback regulation of the pathway. Figure 4b is a schematic representation showing many of the component molecules that make up the JAK/STAT pathway, with expression of all the colored genes being affected by BDNF treatments and all the ones with a purple border reversed by JAK/STAT inhibitor treatments WP + B and RX2 + B.
In addition to the significant relationship to epilepsy identified by IPA (see Additional file 1: Table S1), the overlapping WP/RX2 gene set had a significant link to neuroinflammation despite the neuronal cultures being largely devoid of any glial cells, 41 genes of the BDNFinduced set of 56 neuroinflammation-linked genes were reversed by WP + B or RX2 + B treatment (73% of the gene set) (Fig. 4c).

Functional relationships between genes as described in the literature
Pathway and GO analysis identified the following epileptogenic processes, in addition to inflammation, as being enriched in the WP + B/RX2 + B overlapping data set: synaptic plasticity, receptors and ion channels, epilepsy, transcription factors, proliferation, and neurogenesis. In order to visualize the functional connectivity of the genes involved in epileptogenesis, we created a functional network of the most significant genes as determined by a p-value of < 1 × 10^10 and a fold-change > 2 (Additional file 3: Figure S1). Genes that met these conditions were separated into the lists of the epileptogenic processes, as described above. In order to relate these genes to one another and their epileptogenic functions, we created a network that portrays their function, expression fold change, direction, and functional protein relationships in Cytoscape [42] (Fig. 5a). Using GENEMania, a Cytoscape tool containing a database of experimental data showing protein-protein interactions [43], we determined potential functional relationships. Each node in the network represents a gene while the edges represents a protein interaction. Node color is determined by the epileptogenic category. Due to the interplay between these functions, many genes within one category overlap with another category. For instance, Egr3 plays a role in 5 of the 6 categories: synaptic plasticity, epilepsy, transcription factors, proliferation, and neurogenesis. The color in this case is determined by whichever category comes first in this list: synaptic plasticity, receptors and ion channels, epilepsy, transcription factors, proliferation, and neurogenesis, where the order is determined by the number of genes in that category from smallest to largest. This was done to avoid having the majority of the nodes colored for the largest category (neurogenesis). For Egr3, it is yellow representing synaptic plasticity, however, its location in the network is near the center because the proximity of the node to other categories represents its link to those secondary categories. In addition, the Venn diagram in 5B also shows the degree of overlap for each category. For instance, there are 108 total genes related to receptors and ion channels. In the list, 61 of those did not belong to any other category, 11 were also linked to neurogenesis, 5 to receptors, neurogenesis and proliferation,14 to receptors and epilepsy, 6 to receptors, epilepsy and Differentially expressed genes exclusive to BDNF/RX2 contain a rich set of epilepsy-associated targets not present in the WP1066 exclusive list To further dissect the regulatory contribution of the two JAK/STAT inhibitors, we examined the list of DEGs (as determined by DESEQ2, FDR = 0.05, Wald test) that differed between WP + B and RX2 + B treatments. Four hundred seventeen DEGs are reversed by RX2 + B treatment that are absent with WP + B and 610 DEGs are reversed by WP + B that are not by RX2 + B (Fig. 2a, bottom panel). Interestingly, the genes whose expression is exclusively reversed by RX2 + B are significantly related to neurodegeneration, and brain-related movement disorders (i.e., Huntington's disease) as well as to epilepsy. On the other hand, those exclusively regulated by WP + B were not highly associated with epilepsy or other brain disorders. Additional file 4: Figure S2 contains the EnrichR generated KEGG pathway analysis for this gene set (Additional file 4: Figure S2A). Nine of the 10 top related canonical pathways differ from the KEGG list of overlapping WP/RX2 DEGs (Fig. 4b). Top genes involved in the KEGG pathways listed are also shown in a clustergram (Additional file 4: Figure S2B). The list of genes exclusively regulated by RX2 + B that are linked to epilepsy, Huntington's disease, and neurodegeneration (E) are also presented (Additional file 4: Figure S2C-E).

qRT-PCR gene validation of datasets
Many of the genes identified in these RNA-seq studies were originally reported by us in the literature (ICER, Gabra1, Egr3), using a variety of molecular approaches [23,32,44]. To further validate the new datasets, we generated an additional set of cultures, originating from the embryos of multiple pregnant rats, and chose the following genes for qRT-PCR because of their interest in the field of epilepsy: Dopamine Receptor D5 (Drd5), Galinin Receptor 1 (Galr1), γ2 subunit of the GABAR (Gabrg2), Glutamate metabotropic receptor 1 (Grm1), and their relationship to the JAK/STAT pathway: Inducible cAMP early repressor (Icer), myeloid leukemia cell differentiation protein (Mcl1), Suppressor Of Cytokine Signaling 3 (Socs3), Cyclin D1 (Ccnd1), and C-X-C Motif Chemokine Receptor 4 (Cxcr4) (Fig. 6). In most cases, findings from qRT-PCR were consistent with those of RNA-seq. Quantitated results from select validations are shown (Fig. 6a-d).
A results table for all candidates where validation was performed is included in Fig. 6e. All significance was determined using One-Way ANOVA followed by Tukey's test for multiple comparisons.

HP1α association with STAT3 in neurons
As a first step to uncover the non-canonical mechanism underlying BDNF-induced JAK/STAT signaling in neurons, we asked whether BDNF treatment alters the association of STAT3 with Heterochromatin Protein 1 Alpha (HP1α), a mediator of gene silencing that is found in heterochromatin. An association of STAT3 with HP1α has been previously reported in the literature [49] and we expected to find that BDNF might decrease association thereby releasing the DNA for the subsequent binding of transcription factors and transcriptional elongation. To our surprise, however, BDNF increased association of the two proteins which was blocked by the JAK/STAT inhibitors WP and RX2 (Fig. 8).

Discussion
There is a strong body of evidence suggesting that BDNF signaling through its TrkB receptor contributes in a major way to the development of TLE, as suggested by animal models of epilepsy [11,50] and their potential relationship to human epilepsy brain pathology [51]. BDNF plays a critical role in the cell biological processes that regulate the structural plasticity of dentate granule cells and the enhanced excitation associated with changes in LTP and compromised GABA-mediated inhibition [17,46,50,52]. Although the role BDNF plays in brain disorders is better understood now than at the time of its discovery, many key aspects of its complex regulatory programs remain a mystery. We designed our study to gather a comprehensive list of genes whose expression is either upregulated or downregulated by increased exposure to BDNF, as is believed to occur in response to seizures and other brain insults such as traumatic brain injury [53]. The use of primary rat neurons treated with recombinant BDNF allows us to pinpoint the effects of BDNF without the background of additional signaling pathways that may be activated in TLE models. Our study provides a comprehensive list of the changes in the transcriptome that occur in response to a relatively short but continuous exposure to BDNF (4 h) at a concentration of 0.7 nM that Fig. 7 The WP + B/RX2 + B overlapping datasets contain previously validated STAT3 ChIP-sequencing target genes. List of all WP + B/RX2 + B genes (308) that are validated STAT3 ChIP-seq targets (as determined by Zhang, et.al, 2013) in alphabetical order. Highlighted in yellow are genes of special interest Fig. 8 Binding of HP1 to STAT3 is enhanced by BDNF signaling in neurons. a Representative Western Blot image of coimmunoprecipitation assays using whole-cell extracts of 9-DIV primary cortical neurons treated (30 min) with Water or 0.7 nM BDNF. Protein was precipitated with a specific HP1 antibody (Cell Signaling Technologies #2616) and the elution analyzed by SDS-PAGE. Membrane was immunoblotted with an anti-STAT3 antibody (Cell Signaling Technologies #4904) to detect co-association of HP1 with STAT3, (n = 4). b Quantitation of STAT3 association with HP1 from dataset in "A". c Representative Western Blot image: Inhibitors of JAK/STAT signaling (WP1066 (WP) and Ruxolitinib 2 (RX2)) reduce HP1α antibody precipitation of STAT3. Data was normalized to βactin levels which did not change upon BDNF treatment, * = p < 0.05 is relevant to endogenous levels in neurons that activate TrKB receptors and not p75NTR (unpublished data from our laboratory). A collection of the BDNF regulated genes within the list have been previously identified both in treated cultures and in vivo by our laboratory and others. They include Gabra1, the inducible cAMP early repressor (ICER) (which negatively regulates Gabra1) [23], BDNF itself [38], and multiple ion channel genes, as well as genes that code for important markers of neurogenesis and cell proliferation [54,55].
In addition to being consistent with the findings in the literature, our results reveal new BDNF-induced gene targets that are linked to epilepsy. Continuous exposure to BDNF alters the expression of 194 epilepsy-linked genes identified by IPA. The following genes are of particular interest to us because of their marked level of transcriptional alteration and their involvement in epilepsy: Dopamine receptor D5 (Drd5)-A D1-type dopamine receptor whose knockdown has been shown to be protective against seizures [56]. In our dataset, the neuronal expression of Drd5 increases nearly 15-fold in response to BDNF treatment. Galinin receptor 1 (Galr1)-Galanin is a neuropeptide in the brain that inhibits glutamate release and galanin agonists delivered to the brain can inhibit seizures [57]. Our dataset reveals a 35fold increase in Galr1 expression, demonstrating the complex medley of both protective and pathogenic changes in genome expression that are invoked by increasing levels of neuronal BDNF. Glutamate metabotropic receptor 1 (Grm1)-Grm1 has been associated with the alterations in synaptic plasticity and glutamate disturbances that result in neuronal overexcitation, a hallmark of epilepsy [48,58]. Grm1 mutations can also cause spinocerebellar ataxia type 44, a genetic disease that can present with seizures. In our dataset, BDNF treatment led to a 2.45-fold increase in the expression of Grm1. γ2 subunit of the GABAR (Gabrg2)-the Gabrg2 subunit of GABAR is important for benzodiazepine action and mutations lead to epilepsy [59][60][61][62]. In this study, RNA-seq revealed a 2.38-fold decrease in the expression of Gabrg2 in response to BDNF.
Our new results suggest that sustained levels of BDNF at endogenous concentrations (0.7 nM) can dysregulate a major segment of the epilepsy-associated transcriptome. Recent discovery of the role of BDNF in cancer suggests that understanding how BDNF orchestrates its diverse effects on genome dysregulation will also provide a new window on cancer chemotherapies [63][64][65]. Surprisingly, beyond its activation of downstream targets like the cAMP regulatory element binding (CREB) protein, whose relationship to BDNF has been rigorously described [66,67], we still know little about the role of other BDNF controlled transcriptional regulators and their intracellular signaling partners that contribute to disease.
Using the TLE model of epilepsy, we previously showed that during the latent period after pilocarpine-induced SE, there is a marked upregulation in activated CREB and an induction in the expression of the inducible cAMP early repressor (ICER) [23]. We also reported that activated CREB interacts with ICER to downregulate the levels of Gabra1 transcripts leading to a decrease in α1-containing GABARs in hippocampal neurons. In our studies, we discovered that BDNF was linked to the JAK/STAT pathway as 1) STAT3 knockdown blocked ICER induction in response to BDNF, 2) STAT3 was found on the ICER promoter in response to BDNF treatment (as measured by chromatin immunoprecipitation (ChIP), 3) JAK/STAT inhibitors (pyridone 6 and WP1066) blocked ICER induction and Gabra1 downregulation, and 4) in vivo application of WP1066 at the time of pilocarpine-induced SE reduced the frequency of subsequent spontaneous seizures in these animals [32].
Based on these findings, we asked whether the BDNF transcriptome, as identified in our prior studies, could be separated into those where the BDNF response was blocked by exposure to JAK/STAT inhibitors (Fig. 1) and those that were regulated by other pathways. We chose two inhibitors, WP1066 (a JAK2 inhibitor with diverse mechanisms of inhibition [39]), that we showed was active in primary cultures and in vivo, and Ruxolitinib (Ruxo), a high-affinity JAK1/2 inhibitor that blocks JAK activation and is currently used in the treatment of patients with high-risk myelofibrosis. To our surprise, when we compared the transcriptome of cells exposed to BDNF, and to low dose Ruxo plus BDNF (RX1,100 nM + B), we did not see a significant difference between the two datasets. As we show in Fig. 3, this was not due to any issue with the ability of the stock drug to block STAT3 phosphorylation in sister dishes. Phosphorylation-independent STAT3 activity, however, has been described in multiple contexts [49,68] suggesting that JAK/STAT signaling in neurons may act in a non-canonical fashion, at least when activated by prolonged BDNF signaling as studied in our research. A role for unphosphorylated STATs in genome stability and for nuclear JAK2 in the regulation of histone 3 Y41 phosphorylation has been described [69,70] and is the subject of our future studies.
Despite the lack of response of the transcriptome to BDNF in the presence of the high affinity concentration of Ruxo (RX1, 100 nM), which inhibits STAT3 phosphorylation, WP1066 and Ruxo at 10 μM (RX2) produced substantial changes in gene expression (Fig. 2a). By comparing the two datasets (WP1066 + B and RX2 + B) we were able to decrease the dimensionality of the BDNF transcriptome into those most relevant by their association with epilepsy. In fact, using enrichment analysis, we discovered an association with epilepsy that is more pronounced than any other neurological disease, due to the large proportion of differentially expressed epilepsy-associated genes present within the WP + B and RX2 + B subset (68%) (see Additional file 2: Table S2).
Pathway analysis also revealed that JAK/STAT inhibition had a marked effect on reversing the BDNF-induced transcriptional response for genes organized in pathways that control neurogenesis, neuroinflammation, synaptic plasticity, receptor expression, and synaptic transmission. Alterations in these functions exist in the epileptogenic brain, which is known to have altered neurogenesis, mossy-fiber sprouting, inflammation and neuronal pruning leading to increased plasticity, and altered excitation and inhibition at the synapse [2,3,71]. They also expand upon the recent findings that JAK/ STAT signaling plays an important role in learning and memory [27], suggesting that impairment in such signaling may be implicated in the cognitive comorbidities of the condition. In the future, it will be important to study these changes in the context of an epileptogenic animal model in vivo where the STAT3 neuronal response can be separated from its nonneuronal counterparts in astrocytes and microglia.
Given that multiple epilepsy-associated genes have been shown to have functional binding sites for STAT3, including BDNF itself (Fig. 7), and considering the substantial relationship between prolonged exposure to BDNF and the ability of JAK/STAT inhibitors to reverse the differential expression of epilepsy associated genes in our studies, we hypothesize that identifying the mechanism that connects BDNF signaling to this novel non-canonical role of JAK/STATs in neurons has the potential to reveal new strategies in the treatment of intractable epilepsy, or its comorbidities, and provide a new window on how neuropharmacology can be used to identify the complex cross-talk between different intracellular signaling pathways in different cell types.
As a first step towards identifying the novel mechanism(s) that drive differential gene expression in response to dysregulated BDNF, we discovered that, contrary to our expectation, BDNF increases the association of HP1α with STAT3 (see Fig. 8). Perhaps this is not so surprising as HP1α has also been associated with euchromatin and is particularly important in gene induction [25,26]. Future studies will be aimed at learning more about the functional consequences of such protein/protein interactions and their potential role in disease.

Conclusion
We have employed deep RNA-seq to determine the response of BDNF on the transcriptome of neurons and have employed a neuropharmacological approach to reduce the dimensionality of the datasets to identify BDNF responsive genes whose transcriptional polarity is reversed by two well characterized JAK/STAT inhibitors. Our results show that this subset of genes contains many with previous association to epilepsy in both animal models and human patients. The BDNF-induced JAK/ STAT gene set is highly enriched for genes involved in synaptic neurotransmission, and contains targets from all the major classes of ion channels and neurotransmitter receptors of the brain. In particular, the gene set includes those known to play a role in epileptogenesis through the regulation of synaptic plasticity, neurogenesis, transcriptional regulation, neuroinflammation and proliferation. Some of the genes are known to contain functional STAT recognition sites as verified by EN-CODE, while others remain to be interrogated. Pharmacological analysis also revealed that phosphorylation of STAT3 at Tyr705 most likely does not control BDNF-induced JAK/STAT regulation of genome expression in neurons, suggesting that the mechanism is non-canonical. This conclusion is consistent with the additional observation that BDNF regulates the association of HP1α with STAT3, suggesting that it either regulates heterochromatin structure or plays a novel role in euchromatin. Most importantly, the RNA-seq datasets highlight the importance of the JAK/STAT pathway in neurons and stimulate further discovery in this area based on its potential relevance to brain function and the need to identify new therapeutic strategies for the treatment of intractable neurological disorders.

RNA extraction and library preparation
At the end of the 4-h treatment period, cells were collected and RNA extracted using the Qiagen RNeasy Micro Kit with DNAase treatment (Cat. 74,004). RNA was run in the Agilent Bioanalyzer to determine RNA integrity number (RIN). RIN of all samples was ≥9.0. mRNA was selected by using the NEBNext® Poly(A) mRNA Magnetic Isolation Module (E7490S) before using the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (E7420S) with Agencourt AMPure XP Clean up beads (A63881) to make RNA into strand specific cDNA libraries with multiplexing barcodes from the NEBNext® Multiplex Oligos for Illu-mina® (Index Primers Set 1) kit (E7335S).

Illumina sequencing
Sequencing was performed at the Boston University Microarray Core and University of Chicago Genome Center. Fragmentation and concentration were analyzed on the Agilent Bioanalyzer and Fragment analyzer prior to sample pooling. Runs were either sequenced at 75 bp single-end reads on the Illumina NextSeq sequencer to a depth of 100-120 million reads per sample or at 50 bp single-end reads on the Illumina HiSeq sequencer to a depth of 90-100 reads per sample. All samples were multiplexed and pooled using barcodes to separate biological replicates.

RNA-seq bioinformatics
Data analysis was performed using Strand NGS software, Version 3.2, Build 237,248.© Strand Life Sciences, Bangalore, India. Proprietary Strand NGS algorithms were used to conduct read alignment to the rat genome build rn6 RefSeq Genes (2016.05.11). Reads were filtered based on the following parameters: Max number of novel splices: 1, Min alignment score: 90, Number of gaps allowed: 5, number of matches to be output for each read: 1, exclude reads with alignment length less than: 25, ignore reads with matches more than: 5, Trim 3′ end with average base quality less than: 10. RNA quantification was performed by the Strand NGS software, which counts the total number of reads that map to each gene and exon in the genome, and is reported as raw counts. Normalization of the raw counts was achieved using the DESeq method, which accounts for difference in the total number of reads between samples and log transforms the data [72]. R version 3.1.2 was used with DEseq2 version 1.6.1 to determine differential expression FDR = 0.05, Wald test type. The list of significant genes was then used in subsequent analyses. Disease involvement data was analyzed through the use of IPA with proprietary statistical algorithms (QIAGEN Inc., https://www.qiagenbioinformatics.com /products/ ingenuitypathway-analysis). Cytoscape 3.5.1 software was used with GeneMANIA for functional network generation with [42,43]. EnrichR was used for GO-term analysis and KEGG pathway analysis [40,41].

Real-time quantitative reverse transcriptase PCR (qRT-PCR)
Primary neurons were prepared, treated, and RNA extracted, as described above for RNA-seq. The Applied Biosystems 1-step RNA to Ct (cat. 4,392,938) kit was used and protocol was followed as directed: 50 ng of RNA template, 10 μL of 2x mastermix, 0.5 μL of each 40x primer/probe kit (Ppia and the gene of interest), 0.5 μL of the enzyme mix, and nuclease free water was added to each 20 μL reaction. Each reaction was performed with 2 technical replicates in a 384-well plate. All taqman primer/probe kits were ordered from Thermo Fisher Scientific.