|Home | About | Journals | Submit | Contact Us | Français|
Gene regulatory networks underlie the long-term changes in cell specification, growth of synaptic connections, and adaptation that occur throughout neonatal and postnatal life. Here we show that the transcriptional response in neurons is exquisitely sensitive to the temporal nature of action potential firing patterns. Neurons were electrically stimulated with the same number of action potentials, but with different inter-burst intervals. We found that these subtle alterations in the timing of action potential firing differentially regulates hundreds of genes, across many functional categories, through the activation or repression of distinct transcriptional networks. Our results demonstrate that the transcriptional response in neurons to environmental stimuli, coded in the pattern of action potential firing, can be very sensitive to the temporal nature of action potential delivery rather than the intensity of stimulation or the total number of action potentials delivered. These data identify temporal kinetics of action potential firing as critical components regulating intracellular signalling pathways and gene expression in neurons to extracellular cues during early development and throughout life.
Adaptation in the nervous system in response to external stimuli requires synthesis of new gene products in order to elicit long lasting changes in processes such as development, response to injury, learning, and memory1. Information in the environment is coded in the pattern of action-potential firing, therefore gene transcription must be regulated by the pattern of neuronal firing. Such transcriptional regulation is unique to neurons. The mechanisms that could be responsible for regulating gene expression by the temporal pattern of neuronal firing are largely unknown. Whether such regulation is limited to a few genes having special properties or extends to large numbers of genes by virtue of the dynamic response of intracellular signalling networks and cellular systems influencing transcription is unknown. How responsive the transcriptome is to very subtle alterations in patterns of action potentials remains an open, important question. There is very little data on the genomic response to the temporal features of action potential firing, where timing of action potential bursting is variable and the total delivery of action potentials over a given time period constant. Pharmacological manipulations, e.g., kainate2,3, pilocarpine4, NMDA5,6, KCl7, bicuculline8 and BDNF7 have been used to induce membrane depolarization followed by transcriptome analysis. These studies have detailed differential mRNA abundance, transcription factor regulation, regulatory sequences, chromatin occupancy of gene regulatory sites, and gene networks activated or repressed in neurons upon stimulation. However, a pharmacological approach induces stimulation that lacks temporal control of action potential firing patterns, which is critical during neural development and plasticity6,8,9. In contrast, altering the pattern of the same number of action potentials by using electrical stimulation produces a more naturalistic stimulus required to address the question of whether the abundance of RNAs and transcriptional networks are differentially sensitive to neuronal firing patterns.
We hypothesized that the temporal dynamics of action potentials would regulate the expression of hundreds of genes across broad functional categories10. In order to test the hypothesis that global mRNA abundance in neurons is dependent upon the specific pattern of action potential firing, mouse DRG neurons were electrically stimulated in cell culture with two distinct patterns of activation. These two patterns delivered the same total number of action potentials, but they were grouped into two different durations of firing at 10Hz with different inter-burst intervals. These firing patterns are within the normal range of action potential firing in DRG neurons. Action potentials were delivered at a frequency of 10Hz in 1.8second bursts, repeated at 1min intervals, or stimulated in 9second bursts, repeated at 5min intervals (termed 18/1 and 90/5 respectively). Stimulation was delivered for 2 or 5hr and the abundance of 39,430 mRNA transcripts were compared with unstimulated neurons. DRG neurons were chosen for this study because they have no dendrites and do not form synapses (autapses) with other DRG neurons in vivo or in monoculture. Therefore action potential firing pattern could be controlled by electrical stimulation without complications of neural circuit activity mediated by synaptic connections and neurotransmitters.
Membrane depolarization and action-potential dependent Ca2+ influx by voltage-gated channels serve as the primary mechanism coupling neural activity and control of transcription, transport, and stability of mRNAs in neurons11,12. (In synaptic networks, NMDA receptor and other signalling pathways are also important in excitation-transcription coupling13,14, and local submembranous Ca2+ levels can be appreciable in the absence of measurable changes in cytoplasmic Ca2+). We previously reported the magnitude and temporal kinetics of intracellular Ca2+ responses evoked in DRG neurons by both of these stimulus patterns (18/1 and 90/5) and also found that levels of the immediate early gene (IEG) c-fos did not correlate with the magnitude of the evoked intracellular Ca2+ transient, but instead correlated inversely with the interval of time between bursts of action potentials15. Even a series of single action potentials can drive expression of c-fos if they are delivered at an appropriate interval (10seconds)16. This suggested that gene expression is differentially regulated by differences in temporal dynamics of activation and inactivation of discrete/dissimilar sets of molecules within intracellular signalling networks that ultimately determine gene transcription or degradation of mRNA transcripts. In support of this hypothesis, MAPK phosphorylation, which is an important element in the signalling cascade linking membrane depolarization to gene transcription through intracellular Ca2+ fluxes, and phosphorylation of the transcription factor cAMP-response element binding protein (CREB), showed markedly different responses to these two patterns of stimulation15. MAPK is activated strongly by the 18/1 pattern of action potential firing, but not by the 90/5 pattern, whereas the transcription factor CREB is activated by both temporal patterns15. Thus expression level of specific genes that are influenced differentially by these two molecules should be selectively regulated by these two patterns of impulse firing. This hypothesis is tested here by genome-wide transcriptome analysis of all mRNAs. By extension, many other genes may be differentially regulated by temporal responses of different components in the intracellular signalling network beyond those regulated by CREB or MAPK.
Our results demonstrate that the abundance of hundreds of mRNAs in neurons across many functional gene categories can be modulated by the temporal pattern of neuronal firing. We identified functional gene categories and regulatory sequences in DNA that control activation and repression of these genes differentially in response to these two patterns of action potential firing. These genes include not only IEGs known to be responsive to action potential firing, but also categories of genes involved in neurite outgrowth, nervous system development, plasticity, and other neuronal functions.
We found that the majority of transcripts in DRG neurons were not regulated significantly by either pattern of electrical activity after 2 or 5hrs of stimulation (86.1% of transcripts, 24, 213 out of 28, 131 unique IDs). We identified hundreds of transcripts that were regulated in a stimulus pattern-specific manner (Figs 1 and and2A).2A). Significantly more transcripts were regulated after 5hr stimulation than 2hr of stimulation with the 18/1 firing pattern. In contrast, the 90/5 stimulation pattern regulated more genes at 2hr than after 5hr (Fig. 2B). Some transcripts were regulated similarly by both patterns of stimulation (Fig. 2B and C). The majority of transcripts that were up-regulated by the 90/5 stimulus pattern at 2hr were subsequently down-regulated by 5hr of stimulation. Of the 9 transcripts that were found to be up-regulated at both 2hr and 5hr stimulation at 18/1 and the 90/5 stimulus patterns, several are “core” activity-regulated genes that function as transcriptional regulators (Gadd45b, Npas4, Egr4, and Nr4a1, Nur77, and Ngfib), GTP-binding protein regulators (Gem and Rgs2) and growth factor signalling (Bdnf, Soc3). In contrast, none of the 22 transcripts downregulated in this stimulus/response category have been previously shown to be either activity-regulated or to have an established role in nervous system development. These very different gene expression profiles following electrical activity suggest that distinct gene regulatory networks are recruited by distinct action potential firing patterns (Fig. 2C). We therefore performed pathway analysis of gene-regulatory networks, signalling cascades, and gene ontologies that are recruited by distinct patterns of activity over time.
We identified distinct gene networks that were activated by the two stimulus patterns after 5hr of stimulation (Fig. 2, Tables 1 and and2).2). These included common disease and functional networks after the 18/1 (Table 1) and 90/5 (Table 2) stimulation patterns; e.g., neurological disease, cellular assembly and organization, and small molecule biochemistry. We have previously found that 18/1 and 90/5 stimulus patterns recruit Ca2+ and MAPK signalling pathways differently15,16; this was confirmed by analysis of canonical pathways enriched for both stimulation patterns. Pathways important in axonal growth and growth-cone signalling were oppositely regulated by the 18/1 and 90/5 stimulus patterns (Tables 3 and and4),4), e.g., the genes Cdk5-, Rac-, and Ngf-signalling networks. Pathways and gene networks important in axonal growth and development were oppositely regulated by the two distinct patterns of electrical stimulation (Fig. 3). Nerve growth factor (NGF) and cell adhesion (integrin) signalling pathways functioning through Rac-Rho GTP-binding proteins showed many differentially regulated transcripts (13 focus molecules up-regulated by the 18/1 stimulus and 25 down-regulated by the 90/5 stimulus). This suggests that genes in the pathways that are differentially activated by these two stimulus patterns may share transcription factor binding sites that are preferentially responsive to each of the two patterns of action potentials. We therefore analyzed conserved upstream regulatory elements of activity-regulated genes to identify over- and under-represented TFBS in the genes regulated by each stimulus pattern.
In search for TFBS enriched in genes regulated by specific action potential firing patterns, we analyzed transcriptome data by two techniques: an enhancer identification method through distant regulatory elements of co-regulated genes (DIRE)17 and by gene set enrichment analysis (GSEA)18. The sensitivity of DIRE was validated by analyzing transcriptome data where either CREB-regulated genes (526 genes), NF-κB regulated genes (403 genes), or MAPK-regulated genes (210 genes) were enriched (Fig. S2) (See methods for description of analysis and datasets for CREB, NF-κB, and MAPK). We confirmed the enrichment of several transcription-factor binding sites (TFBS) for CREB, ATF, VJUN, ATF3, and TAXCREB and to a lesser extent, NF-κB TFBS (NF-κB p65) in differentially regulated transcripts. As we expected, NF-κB target genes were enriched for NF-κB specific binding sites (p50, p65, and NF-κB) and Ets219. In contrast, MAPK acts through multiple signalling pathways and may recruit many transcription factors. We identified a more diverse set of TFBS, prolactin receptor20, p5321,22, and TATA (despite being present in 25% of upstream regions), through analysis of down-regulated transcripts by targeted deletion of Mapk in embryonic DRG neurons23 consistent with converging action of several upstream signalling networks.
We applied DIRE analysis to identify upstream regulatory elements enriched in up- and down-regulated transcripts by 18/1 or 90/5 electrical stimulation patterns and performed hierarchical clustering by TFBS site and condition of the abundance and importance factor of 311 TFBS (Fig. 4). Of all the identified TFBS, 53 were present in one condition and 39 had variable abundance across all 8 conditions (Table S3). The most abundant TFBS present in all activity-regulated transcripts were MAF, Nkx2.5, STAT, and SRF. In order to discover sites that may be recruited in a pattern- and time-dependent manner, we applied a background set of non-regulated transcripts to scale TF abundance24. The identified sites may correspond to core, enhancer, or repressor elements. Putative enhancer sites would be expected to be over-represented in up-regulated transcripts; conversely repressor sites would be enriched in down-regulated transcripts. We discovered several top-scoring TFBS (top 10 by importance factor) through analysis of co-regulated genes (Fig. 5). Thyrotroph embryonic factor (Tef), a bZIP-family TF, forms heterodimers with other transcription factors in order to act as either an enhancer or repressor25 TEF sites were enriched in up- and down-regulated transcripts by 18/1 stimulatin (5hr) and 90/5 stimulation (2hr). STAT TF sites (JAK/STAT signalling, cytokines, regeneration) were selectively over-represented in transcripts up-regulated by 18/1 stimulation at 2hr. NF-κB and NF-κB p50 enhancer binding sites were the most over-represented sites in transcripts up-regulated by 90/5 stimulation at 5hr. Conversely, CCAAT-displacement protein CUT (Cdpcr1)26,27 and GATA128,29,30, transcriptional repressors of developmental and synaptic genes, were most over-represented in down-regulated transcripts by the 90/5 stimulation pattern. NF-κB sites were significantly enriched when compared across transcripts (promoter and ECR) or TFBS (NF-κB, NF-κB p50, NF-κB p65) (Table 3).
Recognition sites for serum response factor (SRF) and CREB were significantly enriched by 18/1 stimulation at 2hr (FDR<0.05) (Table S2) when transcriptome data was independently analyzed by gene set enrichment analysis (GSEA). SRF and CREB sites were only observed after 5hr stimulation by the 90/5 pattern. As we observed with DIRE, NF-κB sites (NF-κB, NF-κB p65), were only significantly enriched by the 90/5 stimulus at 5hr. Two TFBS were over-represented across all conditions, FREAC3 (forkhead transcription factor site, FOXC1) and MYOGNF1 (myogenin/nuclear factor 1). These findings were consistent with recruitment of distinct transcriptional networks by patterned action potential firing activity.
To demonstrate specific activation of the NF-kappa transcriptional network by 90/5 electrical stimulation, we identified transcripts through DIRE and GSEA that were enriched by having multiple NF-κB sites. We performed RT-PCR on 3 genes that had enriched NF-κB TFBS: BCL6 Co-Repressor (BCOR), proto-oncogene serine/threonine-protein kinase (PIM1), and NF κB Inhibitor Epsilon (NFKBIE). All 3 transcripts were significantly increased by 5hr of the 90/5 stimulation pattern and were not changed by the 18/1 stimulus pattern (Table S1). We further analyzed the protein expression levels of PIM1 by immunoblotting, which was consistent with the mRNA data. The proto-oncogene serine/threonine-protein kinase, Pim1, controls NF-κB activation through stabilization of RelA-p65 by phosphorylation31. Pim1 was expressed in DRG neurons (Fig. 6) and pre-treatment with the general proteasome and NF-κB inhibitor MG132 (10μM) reduced both basal and activity-induced expression of Pim1.
These studies demonstrate that a critical factor in the regulation of gene expression in neurons is the temporal nature of action potential firing. We have identified large groups of genes whose mRNAs, across many functional categories, can be differentially regulated by both the temporal dynamics of action potential firing and total length of stimulus time. Thus in neurons, the transcriptome in general, and not simply a specialized set of genes, is highly regulated by the state of neural impulse activity. The implications of the broad categories of genes being regulated by the pattern of neural impulse firing is that the effect of other signalling molecules, growth factors, etc., on neurons will vary greatly depending on the state of activity of the neuron. For example, as discussed below, factors influencing gene expression that stimulate neurite outgrowth will have different responses depending on the state of impulse activity in the neuron through the broad activity-dependent changes in the transcriptome. Likewise, by regulating neural activity, neurite outgrowth and other neuronal responses could be regulated independently of signalling molecules or drug treatments. It may be possible to exploit activity-dependent regulation of the transcriptome therapeutically, for example in promoting recovery from neural injury.
We have classified these genes by the occurrence of TFBS in the promoter regions of genes specifically activated by action potential pattern and/or time of stimulation. These analyses reveal that many genes containing binding sites for well-known activity-regulated transcription factors, such as CREB and SRF are regulated without regard to the two different action potential firing patterns (18/1 and 90/5) after 2 and 5hr of stimulation. However, we have found that a subset of TFBS, such as NF-κB, are overrepresented in genes whose abundance was regulated differently by the two action potential firing patterns tested. This strongly suggests that the activation and repression of gene networks and subsequent control of gene expression in neurons, is subject to temporal control by the timing of intracellular Ca2+ transients, rather than the overall concentration of Ca2+ within the cell.
We have shown that two different patterns of electrical stimulation that produced the same number of action potentials regulate distinct gene regulatory networks. Surprisingly, very few transcripts were regulated in a time- and pattern-independent manner (Fig. 2C inset). 8 of the 9 up-regulated transcripts have been previously described as IEGs, such as the transcription factors Egr4 and Npas4, whereas Bdnf is one of the most studied IEG growth factors32. The TRP-M channel, Trpm6 (transient receptor potential cation channel subfamily M member 6), which was up regulated in all stimulus conditions, has an important role in magnesium homeostasis33, but a defined function in sensory neurons has not been shown. Historically, less interest has been placed on identifying genes or networks that are down-regulated by neural activity. The results of the present experiments revealed thousands of genes that are down-regulated in a pattern- and time-dependent manner. However, only 22 transcripts (Fig. 2C) were strongly down-regulated, independent of stimulus pattern. These genes have no shared network or function in common. Many of these genes have a role in transcriptional and translational regulation. However, by TFBS analysis, several predicted TFBS sites associated with repressor function were enriched in these genes; TAL1, Nkx2.5, ZNF21934 and NERF-ELF2 can function as repressors35,36. This is suggestive of transcriptional pathways that require derepression in order to facilitate the expression of activity-regulated genes, a less studied part of the neuronal activity-driven transcriptome. Further work is required to elucidate a role for the repressors which can perhaps lead to a common pathway or pathways of genes required to be downregulated during neuronal plasticity.
Analysis of gene networks using Ingenuity Pathway Analysis (IPA) demonstrates that the temporal dynamics of action potential firing can regulate pathways involved in growth of neurons. Figure 3 shows that the 2 patterns of stimulation used in this study have opposite effects on the abundance of many mRNAs involved in neurite outgrowth through the Rac-GTPase signalling pathway. The 18/1 stimulus pattern appears to be growth promoting whilst the 90/5 stimulus pattern downregulates many of the same genes and pathway components. Many of these same genes have been shown to regulate neurite outgrowth in cortical neurons37. It has been shown previously that growth cone migration and neurite outgrowth are regulated by patterned activity in DRG neurons38.
Our findings support the hypothesis that transcriptional control of gene expression by the temporal pattern of action potential firing occurs through recruitment of distinct networks (co-regulated sites/heterodimer binding sites) of transcription factors and binding sites. Furthermore, these results show that genes with predicted NF-κB TFBS are over-represented in genes regulated by the 90/5 pattern of stimulation (Fig. 4, Table 5, Table S2). Two independent analysis methods (DIRE and GSEA) identified sets of sites in genes that were regulated by both stimulation patterns (18/1 and 90/5) and time (2hr and 5hr). We have classified many distinct transcription factor pathways and their occurrence in genes regulated by patterned action potential firing. Binding sites for SRF were one of the most abundant of identified TFBS by DIRE (14±1%) and enriched following 2hr stimulation with the 18/1 pattern by GSEA (4 predicted TFBS, Table S2). More importantly, several CREB sites (ATF, CREBATF, and TAXCREB) were specifically enriched in the 18/1 pattern at 2hr (Fig. 4 and Table S2). These finding suggest that well characterized activity-dependent transcriptional networks, e.g. SRF alone are not sensitive to patterns of action potentials and require interactions through complexes of factors, e.g. CREB.
We identified thyrotroph embryonic factor (Tef) sites that were over-represented in up- and down-regulated transcripts following 2hr stimulation in the 90/5 pattern. Tef has previously been shown to control circadian rhythm-mediated gene expression. Tef is expressed in DRG neurons and itself downregulated by the 90/5 stimulus pattern at 5hr. TEF, like other TFs, forms heterodimers with other factors to function as enhancer or repressor complexes39,40,41. GATA1 TFBS were over-represented in down-regulated transcripts (Figs 4 and and5),5), may fine-tune TEF to form pattern-specific repressor complexes. However, transcript levels of several GATA family members is at very low levels (GATA 1–4), GATA5 and GATA6, and GATA-d1, 2a (p66α), 2b (p66β) are expressed at moderate to high levels. Despite the variable expression of GATA isoforms in unstimulated DRG neurons, transcript levels of GATA isoforms were globally repressed by patterned activity, e.g. GATA5 by 90/5 at 2hr (−2.3 fold, p<0.005) and 5hr (−1.6 fold, p<0.001). This regulation may occur through conserved upstream regions of GATA genes. Prolactin receptor (PR) TFBS were enriched in ERK-MAPK regulated transcripts (Fig. S2); PR sites were modestly enriched in up-regulated transcripts at 2hr and 5hr, irrespective of pattern.
Our findings predict that the NF-κB transcriptional network and other cis-acting elements, are sensitive to the pattern of action potential delivery. Given the importance of NF-κB in Ca2+-dependent signalling and transcriptional regulation in many cells types including neurons42 it is perhaps surprising that there is little information regarding genomic and network regulation involving the NF-κB family of transcription factors in the nervous system. Recently, Ca2+ spike duration has been shown to regulate NF-κB transcriptional activity in epithelial cells43. Electrical stimulation of the hind paw induces translocation and phosphorylation of NF-κB in small to medium DRG neurons44, supporting these findings. Furthermore, NF-κB activity affects neurite outgrowth in cultured DRG neurons45. Therefore, the NF-κB family of transcription factors could potentially play a key role in stimulus driven Ca2+ -dependent transcriptional regulation in sensory neurons.
Many other mechanisms can influence gene expression in neurons, such as regulation through epigenetic mechanisms, control of mRNA and protein stability and localization of mRNAs and proteins. However, we have demonstrated that action potential firing patterns can differentially recruit transcriptional networks in neurons leading to regulation of hundreds of transcripts.
In summary these results identify the temporal kinetics of action potentials, such as the duration of bursts of action potentials and the interval of time between repeated bursts, as key drivers in regulating gene expression in neurons through the recruitment of specific transcriptional networks to regulate hundreds of genes across all functional gene categories. The data presented here provides a large dataset on mRNA transcripts in neurons that are regulated in abundance by these different patterns of action potential firing. These data should be a valuable resource for many other investigations. The findings greatly increase our understanding of adaptive global gene regulation in response to environmental stimuli and demonstrate that gene expression in neurons is very sensitive to timing of action potential bursting, and that the activation of transcription factor networks is dependent upon the timing of action potential firing independent of the frequency of firing or the total number of action potentials delivered.
All experiments were conducted in accordance with animal study protocols approved by the NICHD Animal Care and Use Committee. For DRG cell culture, multi-compartment Campenot chambers made of Telfon, were affixed by vacuum grease to collagen-coated 35mm cell culture dishes as previously described46,47. Neurons were dissociated from spinal cords of 13.5-day old mouse embryos and plated at a density of 0.5×106 cells per ml into each side compartment in Eagle’s MEM containing 5% horse serum and 100ng/ml nerve growth factor. Non-neuronal cell proliferation was inhibited by treatment with 13μg/ml fluoro-2-deoxyurindine (Sigma, St. Louis, MO) by one day following plating for 5 days. Cultures were subsequently used for experiments 3–4 weeks after which axons extend from the side compartments into the central compartment and could be electrically stimulated.
Prior to an experiment, DRG cultures were visually inspected for acceptable outgrowth under the Teflon barrier and medium was exchanged for medium lacking NGF and sera, stimulating lids were replaced and cultures were left overnight. The following day, DRG cultures were stimulated through platinum electrodes in contact with medium on opposite sides of the barrier. Stimulation parameters and responses to stimulation have been reported previously48. For this study, either 18 actions potentials every min (18/1) or 90 action potentials every 5min (90/5) was delivered by a 6V, 0.2ms biphasic pulse for either 2hr or 5hr. Control (unstimulated) experiments were performed by handling cultures in an identical manner (from the same dissection) as experimental (stimulated) cultures, including exchanging medium and fitting stimulation electrodes, however these culture dishes were not attached to the stimulation apparatus. Unstimulated controls were prepared on the same day as stimulated cultures and RNA was prepared at the same time in both cases.
Total protein was isolated from DRG cultures using M-PER mammalian protein isolation reagent (Thermo Fisher Scientific, Waltham, MA) supplemented with HALTTM protease and phosphatase inhibitors cocktail (Thermo Fisher Scientific, Waltham, MA). Total protein was resolved by SDS-PAGE on 4–12% NuPAGE Bis-Tris gels (Life Technologies, Carlsbad, CA), transferred to PVDF membranes (Immobilon-P, Millipore, Bedford, MA) and blocked in TBS (10Mm Tris-Cl, pH 7.5, 0.9% NaCl) containing 0.1% (v/v) Triton X-100 (TBS-T) and 5% (w/v) nonfat dry milk or 5% (w/v) bovine serum albumin for 2hr at room temperature or overnight at 4°C. Membranes were probed with appropriate antibodies in TBS-T and 5% BSA overnight at 4°C or 4hr at room temperature with shaking. Primary antibodies were visualized with HRP-conjugated secondary antibodies (Amersham Pharmacia Biotech, Piscataway, NJ) at 1:4000 dilution and enhanced chemiluminescence. Primary antibodies used for immunoblotting: GAPDH (Encor, Gainesville, FL) used at 1:2000 dilution and PIM1 (Cell Signaling, Danvers, MA) used at 1:1000 dilution.
Total RNA (2μg) was extracted from unstimulated and stimulated DRG cultures after 5hr using TRIzol (Invitrogen, Carlsbad, CA) and reverse transcribed using Superscript II (Invitrogen, Carlsbad, CA) as previously described49. Semi quantitative PCR was carried out using a Roche Lightcycler using the FastStart DNA Master SYBR Green PCR reaction mix (Roche Diagnostics, Indianapolis, IN) according to manufacturer’s protocol. Data were analysed by the 2(−ΔC (T)) method with respect to unstimulated cultures.
The following primer sequences were used:
GAPDH FP: AATGCATCCTGCACCACCAAC
GAPDH RP: TGGATGCAGGGATGATGTTCTG
NPAS4 FP: TGGCATATTGTCATCAAGACGTG
NPAS4 RP: TTGAAATGGAAACTGGGATCG
EGR1 FP: AAGCACAGGAGGGAAGAGATG
EGR1 RP: TGAAGTCAAAGGGAACAGGAC
NFKBIE FP: AGAGTGAGGAGGAAGAATATG
NFKBIE RP: TTCTCCTTCTTCGTGCATGAC
PIM1 FP: CAACTCATTCCAGACTCCAGG
PIM1 RP: TGACTAAAGTTCCATCGACAGG
BCOR FP: ATGGGTGTATTGTGTGGTGAGG
BCOR RP: ACATACCACTGGATCACGTCC.
Total RNA was harvested immediately after 2 or 5hr stimulation and processed as previously described49. Each control and experimental condition was the result of 5 pooled cultures. All control and experimental manipulations were performed over a short time to minimize the effects of tissue culture variation on gene expression levels. 825ng of differently (Cy3/Cy5) labeled RNAs from two biological replicas were co-hybridized 17h over night at 65°C with each array of 4×44k Agilent 60mer G2519F mouse chips (“multiple yellow strategy”50). Chips were scanned with an Agilent G2539A dual laser scanner at 5μm pixel size/20-bit and raw data extracted with Agilent Feature Extraction software vs 11.1.1. All spots affected by local corruption or with foreground fluorescence less than twice the background in any of the 20 samples were disregarded. Data were normalized by an iterative method alternating intra- and inter-array normalization until the overall maximum error of estimate became less than 5%. Supplementary Figure 1 shows the Pearson correlation coefficients between the gene expression of all 6 sample pairs that can be formed in each condition. The high correlations prove the robustness of the data.
A gene was considered as significantly regulated in neurons subjected to α (=18/1 for 2hr, 18/1 for 5hr, 90/5 for 2hr, 90/5 for 5hr) pattern of stimulation with respect to non-stimulated (X) neurons if:
The composite criterion (1) requires that the absolute expression ratio exceeds the pooled coefficient of variability for the stimulated (pattern α) and non-stimulated neurons. The fold change cut-off was computed separately for each transcript and incorporates both biological variability and technical noise. This condition replaces the arbitrary 1.5 fold change cut-off that may be too stringent for very stably expressed genes (close values in all four replicates) or too relaxed for highly unstably expressed genes50. The p-values of the heteroscedastic t-test were computed with a Bonferroni type correction applied to the groups of spots probing the same transcript51.
A gene was considered as turned on by a particular stimulation pattern if it was not quantifiable (i.e. foreground fluorescence less than twice the background) in any of the unstimulated neuronal samples but quantifiable in all stimulated samples with that pattern. Turned off genes satisfied the opposite criterion i.e. expressed in all not stimulated samples but not expressed in any of the stimulated ones. Raw and normalized gene expression data were deposited into GEO and are publicly available under accession number GSE84872.
Hierarchical clustering and frequency distribution of Agilent microarray expression data was performed using Cluster 3 (http://bonsai.hgc.jp/~mdehoon/software/cluster/)52 and JTree View (http://jtreeview.sourceforge.net/) on Log2 normalized data. Expression data was filtered according to direction (up- and down-regulated), pattern (18/1 and 90/5), and time (2hr and 5hr) where p<0.05 and fold-change >|±1.4|. Gene ontology (GO), network, and pathway analyses was performed using pathway analysis software (Ingenuity Systems, Redwood City, CA). Four way Venn diagrams were constructed using the Java applet Venny, http://bioinfogp.cnb.csic.es/tools/venny/index.html.
Prediction of TFBS was performed using either distant regulatory elements of co-regulated genes (DIRE) or gene set enrichment analysis (GSEA). Sensitivity of DIRE in predicting TFBS enrichment was validated by analysis of gene sets corresponding to (a) VP16 constitutively active CREB (526 transcripts)9, (b) NF-κB regulated transcripts (403 transcripts) (http://www.bu.edu/nf-kb/gene-resources/target-genes/), and (c) sensory neuron specific ERK1/2 conditional knockout (210)23. DIRE accurately predicted the presence of ATF/CREB, NF-κB, and downstream MAPK TFBS (Fig. S2). For DIRE analysis of DRG transcriptome data, a background set of genes was identified as not regulated (fold-change<|1.1|) by either stimulation pattern (18/1 and 90/5) and time (2hr and 5hr). Gene symbols for co-regulated genes and background sets were uploaded to DIRE where the target elements corresponded to the Top 3 evolutionary conserved regions (ECR) and the promoter ECRs; results were reported as % abundance and importance factor. Hierarchical clustering of TFBS was performed as described for Agilent microarray data.
GSEA of gene expression data was performed with BROAD institute stand-alone application against the Molecular Signature Database (MSigDB). Microarray data (4 replicates per condition) was minimally processed (quantile normalized data of 30,564 entries collapsed to 8,546 entries and expression levels±1 S.D. of background). Expression values in each condition were analyzed as control vs. treatment (pattern and time) against a MSigDB comprising 836 gene sets (221 miRs and 615 TFs). MSigDB are pre-computed set of genes with promoter regions±2kb of the TSS containing TFBS of interest. For identified signatures with a false discovery rate (FDR) of 0.25 or less, normalized enrichment score (NES) and FDR are given; in all cases local p-values<0.05.
How to cite this article: Lee, P. R. et al. Gene networks activated by specific patterns of action potentials in dorsal root ganglia neurons. Sci. Rep. 7, 43765; doi: 10.1038/srep43765 (2017).
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We would like to thank Peter Wadeson and William Huffman for DRG neurons in multicompartment chambers, Daniel Abebe for help with experimental animals and Tanya Barrett for editing the manuscript. P.R.L., J.E.C. and R.D.F. are supported by the Intramural Program of the Eunice Kennedy Shriver National Institute of Child Health & Human Development (NICHD) of the National Institutes of Health (NIH). D.A.I. and S.I. were supported by 5R01HL092001.
The authors declare no competing financial interests.
Author Contributions Conception and design of experiments R.D.F., P.R.L. and J.E.C. Performed experiments P.R.L., D.A.I., S.I. Data analysis and bioinformatics J.E.C., P.R.L., D.A.I. Overall supervision of the project R.D.F. All authors contributed to writing and editing of manuscript. All authors have contributed to and approved the manuscript.