Temporal dynamics of neurogenomic plasticity in response to social interactions in male threespined sticklebacks

Animals exhibit dramatic immediate behavioral plasticity in response to social interactions, and brief social interactions can shape the future social landscape. However, the molecular mechanisms contributing to behavioral plasticity are unclear. Here, we show that the genome dynamically responds to social interactions with multiple waves of transcription associated with distinct molecular functions in the brain of male threespined sticklebacks, a species famous for its behavioral repertoire and evolution. Some biological functions (e.g., hormone activity) peaked soon after a brief territorial challenge and then declined, while others (e.g., immune response) peaked hours afterwards. We identify transcription factors that are predicted to coordinate waves of transcription associated with different components of behavioral plasticity. Next, using H3K27Ac as a marker of chromatin accessibility, we show that a brief territorial intrusion was sufficient to cause rapid and dramatic changes in the epigenome. Finally, we integrate the time course brain gene expression data with a transcriptional regulatory network, and link gene expression to changes in chromatin accessibility. This study reveals rapid and dramatic epigenomic plasticity in response to a brief, highly consequential social interaction.

However, we know little about the temporal and spatial dynamics of neurogenomic plasticity in response to social interactions. It is likely that there are waves of transcription associated with perceiving social information, responding to social information, maintaining a behavioral response, recovering from the social interaction and modifying future behavior [14,15]. Static experiments that measure gene expression at a single time point can only catch a glimpse of what is probably a very dynamic and coordinated process.
Studies in development have linked changes in chromatin accessibility with the time course of changes in gene expression and the activity of transcription factors operating within gene regulatory networks [16,17]. This tactic has also proven to be successful for examining acute, short-term responses of other types, for example, in response to pathogens [18,19]. However, whether the same principles apply to behavioral stimuli, and social interactions in particular, is unknown.
Here, we test the hypothesis that a brief social interaction, albeit one with strong implications for fitness, is sufficient to induce transcriptomic and epigenomic responses that change over time. We test this hypothesis in threespined sticklebacks (Gasterosteus aculeatus), a species for which successful territorial defense is critical for Darwinian fitness. Sticklebacks are small fish whose behavioral repertoire has attracted attention since the early ethologists [20]. Freshwater sticklebacks must quickly establish territories because they have a short window of opportunity to breed in the spring, and die at the end of the breeding season. Male sticklebacks typically occur in neighborhoods and function in a dynamic social environment where they vigorously defend individual nesting territories against intrusions by rival males and predators. The territory is the hub of family life, where the father constructs a nest, attracts females to mate and where he rears the offspring without any help from the mother. If a male fails to defend a territory, he will not obtain a mate and he will not produce offspring therefore effective defense of that territory is necessary for reproductive success. Like other territorial animals, male sticklebacks exhibit experience-dependent changes in behavior following a territorial intrusion, as they learn the boundaries of their territory and how to detect and repel intruders [21,22].
We provide evidence that the genome and the epigenome are highly responsive to social interactions during territory defense. We characterize transcriptomic and epigenomic plasticity in response to social interactions by measuring changes in gene expression at three points in time following a brief territorial intrusion using RNA-Seq. We compare expression in two parts of the brain containing nodes in the social decision-making network [23]: diencephalon and telencephalon. The diencephalon includes the hypothalamus-a key integrator of social information with the neuroendocrine system-while the telencephalon is a part of the forebrain, and includes the teleost homolog of the hippocampus. Using these data, we construct a transcriptional regulatory network that links temporal changes in gene expression in different parts of the brain to the activity of transcription factors operating within a gene regulatory network. Finally, we measure changes in chromatin accessibility in response to a social interaction on a genome-wide scale, using acetylated lysine 27 on histone H3 (H3K27Ac) as a marker of accessible chromatin, and link changes in chromatin accessibility to changes in gene expression. We show that many of the same principles that characterize transcriptomic and epigenomic changes unfolding over development [16,17] also apply to the brain's response to brief, but potent, social behavior.

Spatiotemporal dynamics of the transcriptomic response to territorial challenge
Within both brain regions, we identified genes whose expression was influenced by a territorial challenge at three time points: 30, 60 and 120 minutes (S1 Table). The greatest transcriptional response to a territorial challenge occurred 60 minutes after the challenge (Fig 1A). Generally, gene expression was down-regulated 30 and 60 minutes after the social challenge but was up-regulated at the 120 minute time point in diencephalon ( Fig 1A).
The transcriptomic response to a territorial challenge changed rapidly over time in both brain regions (Fig 1B and 1C); in fact, there was little overlap between the differentially expressed genes detected at each time point. Functional enrichment analysis revealed that early responding genes (30 minutes) were related to hormones and post-translational modifications (PTMs), whereas a strong signature of genes related to metabolism dominated at the 60 minute time point (consistent with [24]). By 120 minutes, differentially expressed functions shifted toward transcription, immune response and homeostasis (Fig 1G, S2 Table).
Not surprisingly, we detected strong differences in gene expression between brain regions. However, there were some genes that were differentially expressed in both brain regions, and these genes showed a remarkably concordant quantitative pattern of expression across brain regions at 30 and 60 minutes (Fig 1D and 1E). Specifically, genes that were strongly upregulated in diencephalon in response to a territorial challenge were also strongly upregulated in telencephalon (correlation > 0.9). The pattern at 120 minutes was different, with a subset of genes (n = 18, hereafter referred to as 'discordant genes', S3 Table) showing the opposite pattern of regulation in the two brain regions. These 18 genes, which were upregulated in diencephalon and downregulated in telencephalon after the territorial challenge (Fig 1F), are primarily related to visual perception and include retinal genes (e.g. rom1b, rom 1a, opsins), and down (red)regulated genes at 30, 60 and 120 minutes after a territorial challenge in diencephalon and telencephalon. Overlap between differentially expressed genes across time points in diencephalon (B) and telencephalon (C). Correlation between expression in diencephalon (X axis) and telencephalon (Y axis) at 30 min (D), 60 min (E) and 120 min (F) after a territorial challenge. The numbers in the Venn diagram indicate the number of differentially expressed genes in each brain region and the overlap between them at a given time. Scatterplots show the expression pattern of the genes that were shared between brain regions at circadian genes (e.g. crx, opsins) and phosphodiesterases, which have been repeatedly duplicated in the stickleback genome and acquired new functions [25].

Waves of transcription in response to a territorial challenge
In order to find genes that changed in a coordinated fashion in response to a territorial challenge, we first analyzed the gene expression data by testing for main effects and interactions between them. We built separate generalized linear models for each brain region and were particularly interested in genes whose time course of expression was influenced by the territorial intrusion (time x treatment interaction term). There were 758 and 739 such genes in diencephalon and telencephalon, respectively, hereafter referred to as DEG x (FDR < 0.1, S4 Table).
We next used hierarchical clustering of the DEG x to determine whether there were clusters of genes that changed in concert together. We identified 12 and 13 clusters in diencephalon and telencephalon, respectively (Fig 2A and 2B, S4 Table). Each cluster had a particular expression profile over time in response to a territorial challenge. For example, cluster D1 comprised a set of genes that were downregulated at 30 mins, upregulated at 60 mins and downregulated again at 120 mins. On the other hand, cluster D2 comprised a set of genes that were upregulated at 30 mins, downregulated at 60 mins and then strongly upregulated at 120 mins.
Functional analysis revealed that genes with similar time-course profiles also tended to have similar functions. For instance, cluster D1 comprised hormonal genes such as tshb, prl, cga, lhb and gh1, and a nuclear receptor transcription factor nr5a1b, which binds to a prl (encoding prolactin) promoter [26]. In contrast, cluster T2 included transcription factors such as pax7 (both a and b paralogs), irx (2a, 3a and 5a), tfap2b, shox, sp5l, DMBX1 and pou4f2 and two hormonal genes known to be very important to social behavior (avp and oxt). The clustering of hormonal genes with transcription factors suggests a complex interplay between hormones with transcription factors in response to a territorial challenge. Clusters D2 and T3 included the genes that exhibit the discordant pattern of expression across brain regions at 120 minutes ( Fig 1F).
Functional enrichment analysis confirmed that each cluster of genes is associated with its own unique set of functions ( Fig 2C and 2D, S5 Table). GO terms enriched in cluster D10, for example, were not shared with any other cluster, while cluster D4 also associated with its own unique set of GO terms. These results are consistent with the hypothesis that there are waves of transcription associated with different biological functions following a social interaction. Some biological functions peak early then subside (e.g., cluster T3), while others peak at 60 minutes (e.g., D1, D8, T5, T8); still others peak hours (120 mins, e.g., D2, D4, T1, T4) following a social interaction.
Transcription factors within a gene regulatory network coordinate waves of transcription Next, to identify genes that regulate transcriptional changes in the stickleback brain in response to a territorial challenge, we reconstructed a transcriptional regulatory network (TRN) model using the ASTRIX approach [27]. We used the gene expression data to identify regulatory interactions between transcription factors and their predicted target genes (see methods). ASTRIX infers a genome-scale TRN model capable of making quantitative a time point. Note the cluster of genes in the lower right corner of 1f, hereafter referred to as 'discordant genes', which were differentially expressed in both brain regions at 120 minutes but in opposite directions: they were upregulated in diencephalon and downregulated in telencephalon. (G) Functional enrichment of DEGs by time point (columns) and by brain region (rows), shown as revigo-like MDS graphs. Blue indicates enrichment of up-regulated genes, red indicates enrichment of down-regulated genes. Groups of terms with similar functions are highlighted. predictions about the expression levels of genes given the expression values of the transcription factors. The full TRN is in S1 Fig. We then integrated the DEG x from the hierarchical clustering analysis with the TRN in order to find transcription factors that are predicted to regulate the clusters. This integration proved to be insightful because it connected dynamic gene expression to interacting transcription factors within a gene regulatory network. For example, the transcription factors dlx4a, grhl3 and si:ch211-157c3.4 were predicted to regulate cluster D9 (enriched for energy metabolism and immune response) and were connected to each other in the network. This analysis therefore allows us to identify transcription factors within a gene regulatory network that we hypothesize are regulating clusters of genes that change in a coordinated fashion in response to a territorial challenge (Fig 3).
Indeed, closer examination of the dynamics of expression of transcription factors and their targets revealed that many of the transcription factors in the TRN showed expression patterns consistent with the cluster they were predicted to regulate. For example, the expression pattern over time of irf8 and cebpb was very similar to the expression pattern of their predicted targets (D4, D5, D6, D9: no change, down, up).
The TRN offered a number of insights into the spatiotemporal dynamics of gene expression in response to a territorial challenge. For example, the TRN can help explain striking patterns in the gene expression results, such as the discordant genes that were upregulated in diencephalon and downregulated in telencephalon 120 minutes after a social challenge ( Fig 1F). The discordant genes were in clusters D2 and T3, which were predicted to be regulated by the set of connected transcription factors otx5, vsx1 and crx. Crx, implicated with circadian rhythm in addition to visual functions [28], is noteworthy because its expression profile was consistent with the expression pattern of the discordant genes: crx was upregulated in diencephalon and downregulated in telencephalon at 120 minutes. A full list of the TFs in the TRN that were enriched in the DEG x is in S6 Table. Linking changes in gene expression to changes in chromatin accessibility While chromatin is suspected to change relatively slowly compared to mRNA in adult tissues, changes in chromatin accessibility can be an important driver of changes in gene expression [29]. However, little is known about the impact of short-term behavioral interactions on the chromatin landscape.
To test the hypothesis that a brief social interaction has consequences for the epigenome, we used chromatin immunoprecipitation on histone H3 subunits with acetylated lysine 27 (H3k27Ac ChIP-Seq), a marker of accessible chromatin, to assess changes in genome-wide chromatin accessibility at two time points (30 and 120 minutes) following a territorial challenge in diencephalon. These experiments revealed tens of thousands of H3K27Ac peaks in each sample tested with robust p values and enrichment (S7 Table; see Methods for details). We distinguish between areas of the genome that were accessible in controls ('baseline accessible peaks') from areas of the genome whose accessibility changed in response to a territorial challenge (i.e. differed between control and experimental males, 'differentially accessible peaks', DAPs).
Most of the genes were accessible at baseline (S2 Fig). There were 23656 and 18797 baseline accessible peaks (!4-fold change in peak difference between sample and input, and p < 10 −4 ) associated with 12630 and 11723 genes (within 20kb) at 30 minutes and 120 minutes, respectively (S8 Table).
However, there were a large number of genes whose accessibility was affected by a territorial challenge, particularly 120 minutes following the challenge (Fig 4A and 4B, S8 Table). There were 2868 differentially accessible peaks (DAPs) that were associated with 1975 genes (within 20kb, DAPs; 2-fold and p < 10 −4 ). Representative DAPs are shown in Fig 4D, 4E and 4F. Many (n = 97) of the peaks that were differentially accessible at 120 minutes were near genes whose expression profile changed over time in response to a territorial challenge (DEG x, Fig 4C). The DEG x associated with nearby DAPs (hereafter referred to as DAPDEG x ) were not a random set, but were enriched for specific functional categories: functions related to stimulus response, cell signaling and development were highly enriched in this gene set (S9 Table). Each node represents a TF. Slices of pie correspond to different clusters in diencephalon or telencephalon; the key to the clusters is in the lower left corner. A full orange slice represents a diencephalon cluster. A purple half slice represents a telencephalon cluster, a purple and orange slice represents clusters in both brain regions. For example, cebpb is predicted to regulate D4, D5, D6, T3, T4, T9. The territorial challenge had dramatic consequences for chromatin in terms of peak size, with fold enrichment or depletion in specific peaks after challenge as high as 30-40x ( Fig 4A). There was a general trend toward more accessibility at the 120 min following a territorial challenge, consistent with the general pattern of up-regulation of DEGs in diencephalon at the 120 min time point (Fig 1A).
Computational analysis suggests a small set of peaks (associated with 24 genes) that were inaccessible in the control group but which became accessible in response to a territorial challenge at 120 minutes, and are possible pioneer factors (S9 Table). A representative sample is in Fig 4F, which shows a differential peak within 5kb upstream of C4B. The peaks that were inaccessible at baseline but became accessible in response to a territorial challenge included several that are near genes associated with the immune response, e.g. irg1, lcp1, ccr8.1, pstpip1b, PRF1, C4B, zc3h12a. Table 1 illustrates how changes in chromatin are linked to the activity of transcription factors in the TRN and the expression of their targets over time. All of the transcription factors in Table 1 are in the TRN and the genes encoding these TFs were all associated with DAPs that either became accessible, or became more accessible, in response to a territorial challenge at 120 minutes. All but one of these TFs (NFATC3) regulate clusters of genes that are upregulated at 120 minutes (e.g. clusters D9, D12, D6, D12, D5). Several transcription factors (NFATC3, Overlap of genes whose expression profile changed over time in response to a social interaction (DEG x ) with genes associated with differentially accessible peaks at 30 minutes and 120 minutes. The overlap between DEG x and accessibility at 120 minutes is statistically significant (P<0.0001). (D-F) Examples of differentially accessible peaks around DEG x . Separate tracks are shown for H3K27Ac peaks in control 30 min, experimental 30 min, control 120 min, experimental 120 min, and H3k4Me3, which marks the location of the promoter. (D) Pparg (a TF in D9 and also present TRN) was more accessible at 120 minutes and was also up-regulated at 120 minutes. (E) P2ry12 (cluster D9) is purinergic receptor involved in synaptic plasticity [65] that was more accessible in controls at 30 minutes then become more accessible in experimental animals at 120 minutes. P2ry12 is known to stimulate microglia migration toward neuronal damage [66]. (F) C4B (cluster D9) was not accessible at baseline but became accessible at 120 mins. irf8, pparg and cebpb) were themselves differentially expressed over time, and their expression tracks the expression of their targets. The overall pattern of chromatin becoming more accessible at 120 minutes suggests that there are transcriptomic consequences of a brief territorial challenge that persist for more than two hours afterwards.

Discussion
By integrating different types of transcriptomic and epigenomic data (RNA-Seq, H3K27Ac ChIP-Seq) with rigorous computational analyses, we show heretofore underappreciated consequences of social interactions for the brain transcriptome and epigenome. We detected dramatic changes in gene expression over time in response to a brief territorial challenge: clusters of genes enriched for particular biological functions changed in a coordinated fashion, over a period extending for hours afterwards. Using a TRN and generalized linear model, we linked changes in gene expression to the activity of transcription factors, which we propose to be factors that regulate them. Moreover, we demonstrate that a brief social interaction was sufficient to cause changes in the accessibility of chromatin elements throughout the genome, including many linked to DEGs. While conventional wisdom is that chromatin changes relatively slowly in adult tissues, there is some precedent for our findings of rapid response in adult brain; for example, epigenetic responses to strong stimuli such as cocaine can happen quickly, e.g. within an hour [30]. Indeed, there is growing evidence from the learning and memory literature that chromatin can be very dynamic [31][32][33][34], and changes in chromatin accessibility in response to a social challenge have been reported in other species [35]. The magnitude of epigenomic plasticity we observed in response to a territorial challenge is also noteworthy. Hundreds of genes were closely linked to differentially accessible chromatin peaks, and for many of these we found dramatic differences in the degree of accessibility, measured by enrichment for H3K27Ac, following a social interaction (Fig 4B). Indeed, a territorial challenge was sufficient to cause some genes that were not clearly associated with accessible chromatin prior to a territorial challenge to become accessible afterwards (S10 Table). We hypothesize that a territorial intrusion provoked dramatic responses at the transcriptomic and epigenomic level in male sticklebacks because successful territory defense is so consequential in this species, with strong implications for fitness. Table 1. Integrating TFs with DEG x and chromatin accessibility. These TFs are in the TRN and are enriched in the DAPDEG x with accessibility indicated. Some of the TFs (in bold) were differentially expressed and in a cluster. The general expression pattern of their cluster is indicated. A complete set of enrichments is in S10 Changes in gene expression over time were similarly dramatic, consistent with the hypothesis that there are waves of transcription associated with different temporal aspects of behavioral plasticity-some genes are involved in the initial reaction to a conspecific, others with responding to social information and still others involved in recovery and preparing for the future. The early hormonal response parallels the time course of the neuroendocrine response to aggression, which involves both the hypothalamic-pituitary-adrenal (interrenal in fishes) axis and the hypothalamic-pituitary-gonadal axis in vertebrates [36], including in sticklebacks [37]. Interestingly, prolactin-a hormone associated with maternal care-was also recruited in response to a territorial challenge. This result is consistent with the hypothesis that aggression and parental care share common neuroendocrine and neurogenomic substrates [36]. The relatively large number of upregulated DEG and increased chromatin accessibility at the 120 minute time point implies that much of the neurogenomic response to a brief territorial intrusion is related to recovery and preparing for the future, rather than producing the immediate behavioral response [38]. The increased accessibility and expression of genes related to immunity and learning at 120 minutes is consistent with this idea. For example, GO terms related to learning (calmodulin binding and calcium ion binding, involved in the activation of CamK) were enriched in clusters of genes that peak at 120 mins (D2 and D3). Similarly, the expression of CAMKK1 (important for long term memory [39], for example, changed over time in response to a territorial challenge and was upregulated at 120 minutes. Finally, actin binding, important for late long term potentiation and long term memory [40], was implicated in the differentially accessible genes (Fig 4C). Indeed, there is growing appreciation for the relationship between immunity and learning, especially spatial learning [41,42]. Increased chromatin accessibility at 120 minutes is also consistent with the idea that the transcriptomic response to social interactions might be even faster in the future, i.e. priming.
The involvement of learning and memory-related genes makes sense in light of the biology of territorial animals [43]. During an intrusion, territory holders gain information about the spatial boundaries of their territory, the competitive ability of their neighbors and their own resource holding potential, and use that information to guide future behavior. Indeed, territorial animals improve their ability to detect and repel intruders with experience [44][45][46][47] and the behavioral literature is rife with examples of experience-dependent changes following a territorial challenge such as the prior residency advantage [48], the winner effect [49] and the dear enemy phenomenon [50]. Social interactions during territory defense are especially likely to influence spatial learning. For example, fishes actively patrol sites where they've had previous encounters with intruders [51,52]. We speculate that a brief territorial challenge triggers the expression of learning-related genes and that changes in chromatin are associated with the formation of memories of where the social interaction occurred.
A growing number of studies are implicating metabolic genes with aggression [24,[53][54][55]. Consistent with this, the expression and accessibility of peroxisome proliferator activated receptors gamma (pparg), which participates in the regulation of lipid metabolism and glucose homeostasis, changed over time in response to a territorial challenge. Cebpb is another transcription factor enriched in the DAPDEG x which directly binds at the pparg promoter and can regulate its expression [56]. Other studies have shown that pparg is expressed in the hypothalamus and is important for CNS energy balance [57,58]. For instance, pparg agonists, which are insulin-sensitization drugs such as thiazolidinedione (TZD), are widely prescribed to diabetes mellitus 2 patients [57]. Pparg and its targets are downregulated at 60 min and then up-regulated at 120 minutes, possibly reflecting changes in energy balance and homeostasis following an aggressive interaction.
From an ethological perspective it is staggering to consider these results in light of the richness of social animals' lives. Animals that live in social groups are constantly engaged in social interactions. Indeed, rates of territorial intrusions in natural populations of sticklebacks have been reported to be as high as 76 intrusions per hour [59]. Moreover, territory holders interact not only with competitors but also with predators, potential mates and offspring. How animals in natural populations behave during these interactions influences their current and future social environment as well as their fitness. Our results prompt the intriguing hypothesis that meaningful social interactions (even brief ones) can provoke waves of transcription and changes to the epigenome which lead to changes in neural functioning, and those changes are a mechanism by which animals update their assessment of their social world.

Animals
Adult males were collected from Putah Creek, a freshwater population, in spring 2013 and maintained in the lab on a 16:8 (L:D) photoperiod and at 18˚C in separate 9-liter tanks. Males were provided with nesting material including algae, sand and gravel and were visually isolated from neighbors. All males were in the 'territorial' phase of the nesting cycle, i.e. defending a territory. Sneaking is rare in this population.

Territorial challenge
We employed a method to simulate a territorial challenge initially developed by van Iersel [60] and used in previous studies [13,24]. Males were randomly assigned to either the experimental or control group. Males in the experimental group were presented with a smaller, unrelated male intruder confined to a flask. Males in the control groups were presented with an empty flask. At the same time as a confined intruder was introduced to an experimental male's tank, an empty flask was introduced into a paired control male's tank. After 5 min the flask was removed, and after a predetermined period (see below) males were quickly netted and sacrificed by decapitation within seconds following an IACUC approved protocol (#15077) of the University of Illinois at Urbana-Champaign.

RNA sequencing
Tissue preparation. Males for RNA Sequencing were collected 30, 60 or 120 minutes after the flask was introduced, with n = 10 males per time point. Heads were flash frozen in liquid nitrogen and the telencephalon and diencephalon were carefully dissected and placed individually in Eppendorf tubes containing 500 μL of TRIzol Reagent (Life Technologies). Total RNA was isolated immediately using TRIzol Reagent according to the manufacturer's recommendation and subsequently purified on columns with the RNeasy kit (QIAGEN). RNA was eluted in a total volume of 30 μL in RNase-free water. Samples were treated with DNase (QIA-GEN) to remove genomic DNA during the extraction procedure. RNA quantity was assessed using a Nanodrop spectrophotometer (Thermo Scientific), and RNA quality was assessed using the Agilent Bioanalyzer 2100 (RIN 7.5-10). RNA was immediately stored at −80˚C until used in sequencing library preparation.
Library preparation. Poly-A RNA was enriched from 1-2 μg of total RNA by using Dynabeads Oligo(dT)25 (Life Technologies), following the manufacturer's protocol. Two rounds of poly(A) enrichment were performed with a final elution in 14μL of water. The poly-Aenriched RNA was used to prepare RNA-Seq libraries, using the NEXTflex Directional RNAseq Kit (dUTP based) with Illumina compatible adaptors (Bioo Scientific). Manufacturer's instructions were followed and 13-15 cycles of PCR amplification were performed depending on the starting input of total RNA. Libraries were quantified on a Qubit fluorometer, using the dsDNA High Sensitivity Assay Kit (Life Technologies), and library size was assessed on a Bioanalyzer High Sensitivity DNA chip (Agilent). Libraries were pooled and diluted to a final concentration of 10 nM. Final library pools were quantified using real-time PCR, using the Illumina compatible kit and standards (KAPA) by the W. M. Keck Center for Comparative and Functional Genomics at the Roy J. Carver Biotechnology Center (University of Illinois). Single-end sequencing was performed on an Illumina HiSeq 2500 instrument by the W. M. Keck Center for Comparative and Functional Genomics at the Roy J. Carver Biotechnology Center (University of Illinois). The samples were sequenced on 20 lanes.

ChIP sequencing
Tissue preparation. Diencephalons from a new set of males were collected for ChIP-seq at 30 or 120 minutes after the flask was introduced. Prior to nuclei isolation, brain tissue was pooled into groups of 5 and kept at 0˚C in PBS with Protease Inhibitor Cocktail (PIC, Roche 04693132001). Tissue was homogenized by motor pestle and then fixed in PBS+PIC with 1% formaldehyde for 10 minutes. The fixing reaction was stopped with addition of Glycine to a final concentration of 0.125M. Fixed cells were washed 2x with PBS+PIC to remove formaldehyde. Washed cells were lysed to nuclei with L1 lysis solution-50 mM Tris-HCl (pH 8.0), 2 mM EDTA, 0.1% v/v NP-40, 10% v/v glycerol, and protease inhibitor cocktail (PIC)-for 30 minutes on ice. Cell debris was washed away with PBS and PIC. Nuclei were then pelleted and frozen on dry ice. Prior to pelleting, a small (2 μL) sample of nuclei was taken, stained with Trypan, and checked for quality and quantification via hemacytometer. Nuclei were sonicated at high power for 7 x 7 minute cycles (30 s on, 30 s off) in a Diagenode Biorupter Sonicator (Diagenode). Remaining cellular debris was pelleted by centrifugation for 10 minutes at 13,000 x g.
Fragmented chromatin was processed in duplicate for histone H3K27Ac ChIP with Diagenode iDeal ChIP kits, according to manufacturer's specifications with minor adjustments. Six million nuclei were used for each IP. 25 μl of each IP was reserved for input samples. Technical replicate inputs were pooled to 50 μl. 2 g of H3K27Ac antibody (Abcam ab4729) was used for each IP. An additional wash in TE buffer was performed after the initial four IP washes.
Library preparation. After ChIP, IP DNA was quantified by Qubit with a dsDNA High Sensitivity quantification kit (Invitrogen). Libraries were prepared using KAPA LTP library kits, with protocol as written, using Bioo index adapters. Libraries were size selected using AmpureXP beads, with protocol as written, selecting for DNA between 200-500bp in size. Library quality was checked by a Qubit flourometer and Bioanalyzer. Samples were sequenced with an Illumina HiSeq 2500 sequencer.

ChIP-Seq informatics
Libraries from each technical replicate and the input control were sequenced with average depth of 7.6 M reads with average quality score > 35. Technical replicates were pooled with average sequence depth of 16M reads; numbers and QC scores for reads and peaks are summarized for each sample and pools in S11 Table. Sequence data were mapped with Bowtie2 [61] to the Gasterosteus aculeatus reference genome (the repeat masked reference genome, Ensembl release 75), using default settings, yielding 3.87-7.96 uniquely mapped reads (averaging approximately 1X whole genome coverage of the stickleback genome for each replicate). Mapped sequence data were analyzed for peaks using HOMER (Hypergeometric Optimization of Motif EnRichment) v4.7 [62]. Samples were converted into tag directories, and QC was performed using read mapping and GC bias statistics. Histone peaks were then called from the Tag Directories with default factor settings, except local filtering was disabled (-L 0) and input filtering was set at three-fold over background (-F 3), to increase the sensitivity of the peak calling and identify individual subunits of multi-histone peaks, identifying tens of thousands peaks for each sample with average tag counts ranging from 42.7-58.1. Replicates were assessed for correlation, displaying >80% correlation in these filtered peaks across the two samples, which were then pooled for final peak identification. Peaks were highly associated with annotated gene promoters with average distance to transcription start sites (TSS) ranging from 75-328.4 bp (S11 Table), as expected for H3K27Ac [62]; these data confirmed the robustness of the ChIP-Seq data. After peak calling, peak files were annotated to the stickleback genome using HOMER's annotation script to assign peaks to genes, and associate peaks with differential expression data. BigWiggle pileup files were generated using HOMER's makeBig-Wig.pl script with default settings.

Defining differentially expressed genes (DEGs)
HTSeq read counts were generated for genes using stickleback genome annotation. Any reads that fell in multiple genes were excluded from the analysis. We included genes with at least one counts per million (cpm) in at least two samples. Count data were TMM (trimmed mean of M-values) normalized in R using edgeR. To assess differential expression a nested interaction model (~time+treatment:time) was fitted separately for diencephalon and telencephalon in edgeR (see edgeR manual section 3.3.2). A tagwise dispersion estimate was used after computing common and trended dispersions. Finally, to call differential expression between treatment groups, a 'glm' approach was used. We FDR-adjusted the p-values from all contrasts at once. A FDR cutoff < 0.1 was used to call for differentially expressed genes.

Hierarchical clustering analysis
An agglomerative clustering was done separately on DEG x from each brain region. A hierarchical dendogram was generated using hclust function in R (R version 3.2.2), whereas "ward. D" objective criterion was used to merge the pair of cluster at each step. Trees were cut at height 25 to obtain clusters. Each cluster's fold change values at each time point were plotted as profile plots using ggplot2 in R.

Defining differentially accessible peaks (DAPs)
H3K27Ac peaks and their differences between experimental and control groups were calculated at 30 minutes and 120 minutes in each brain region using HOMER's getDifferentialPeaks functionality. For each time point and brain region, two sets of results were calculated: one treating the experimental group as background and the other treating the control group as background. An H3K27ac peak was termed to be "differentially accessible" if it had a fold change of larger than 2 in either set of results, and if it had a p-value less than 10 −4 . Differential peak sets were then annotated using a custom R script to search for all transcripts with transcript start or end sites within 20 kb on all Ensembl-annotated splice variants built using biomaRt.

Associating DAPs with genes
A chromatin domain was defined for each gene in the Ensembl build (v.1.75) of the stickleback genome. First, for each transcript corresponding to a gene, a window was defined that began 20 kb upstream of the transcription start site and ended 20 kb downstream of the transcription end site. A 20 kb window was chosen based on the estimated intergene interval in the stickleback genome. Next, this window was truncated so that it did not intersect with any transcript of any other gene. The union of these windows for all transcripts of a gene constituted that gene's domain. All peaks that had any overlap with the domain of a gene were considered as potential regulators of that gene's expression.

Transcriptional regulatory network (TRN) analysis
ASTRIX uses gene expression data to identify regulatory interactions between transcription factors and their target genes. A previous study validated ASTRIX-generated TF-target associations using data from ModENCODE, REDfly and DROID databases [27]. The predicted targets of TFs were defined as those genes that share very high mutual information (P < 10 −6 ) with a TF, and can be predicted quantitatively with high accuracy (Root Mean Square Deviation (RMSD) < 0.33 i.e prediction error less than 1/3rd of each gene expression profile's standard deviation. The list of putative TFs in the stickleback genome was obtained from the Animal Transcription Factor Database. Given TFs and targets sets, ASTRIX infers a genomescale TRN model capable of making quantitative predictions about the expression levels of genes given the expression values of the transcription factors. The ASTRIX algorithm was previously used to infer a TRN model for the honeybee brain that showed remarkably high accuracy in predicting behavior-specific gene expression changes. ASTRIX identified transcription factors that are central actors in regulating aggression, maturation and foraging behaviors in the honey bee brain [27]. Transcription factors that are predicted to regulate a cluster (from the hierarchical clustering analysis) were determined according to whether they had a significant number of targets in a cluster as assessed by a Bonferroni FDR-corrected hypergeometric test. TFs with at least 3 targets were used (S12 Table) and a FDR cutoff of < 0.05 was used to call for significant associations.

Functional analysis
We derived GO assignments, using protein family annotations from the database PANTHER [63]. Stickleback protein sequences were blasted against all genomes in the database (PAN-THER 9.0 *85 genomes). This procedure assigns proteins to PANTHER families on the basis of structural information as well as phylogenetic information. Genes were then annotated using GO information derived from the *82 sequenced genomes in the PANTHER database.
GO analysis were performed in R using TopGo v.2.16.0 and Fisher's exact test. A p-value cut off of <0.01 was used to select for significantly enriched functional terms wherever possible. For visualization we found dissimilarity among GO terms using zebrafish as closest organism and "Wang" algorithm in GOSemSim package [64]. We then plotted their similarity using the non-metric isoMDS function in MASS. We used the individual terms and the genes inside each term to manually annotate names for clusters appearing in MDS plots. This study has been submitted to GEO. The RNA-Seq data and ChIP-Seq data are accessible with this GEO ID: GSE96673.
Supporting information S1 Fig. TRN highlighting top hubs in the network and region-specific DEGs. The TRN contains 352 TFs, which regulate 1155 genes through a total of 3683 interactions. The top 20 TFs ("hubs") with the highest number of targets (over 30 each) are highlighted in pink. Target genes that are differentially expressed in Diencephalon or Telencephalon (CFDR < 0.1) are shown as blue or green nodes respectively. Genes in the TRN that were differentially expressed in both the regions are shown in orange.   Table. Differentially expressed genes whose expression profile changed over time in response to a territorial challenge (DEG x ), and their associated cluster with separate tabs for each brain region. (XLSX) S5 Table. Functional enrichment of DEG x (biological processes) with separate tabs for each brain region. GO terms were further summarized according to their biological significance in the 'comments' column. (XLSX) S6  Table. Functional enrichment of DAPs with separate columns for DAPDEG x and genes near peaks that were not accessible at baseline but became accessible in response to a territorial challenge (!AccDAPDEG x ). There are separate tabs for biological processes and molecular function. (XLSX) S10 Table. Integrating TFs with DEG x and chromatin accessibility. These TFs are in the TRN and are enriched in the DAPDEG x with accessibility indicated. Some of the TFs (in bold) were differentially expressed and in a cluster. The general expression pattern of their cluster is indicated. A subset of this table is in Table 1. (DOCX) S11 Table. H3K27Ac ChIP-Seq samples read mapping and quality control information. (XLSX) S12 Table. Transcription factors and their targets as predicted by ASTRIX generated TRN. (XLSX) S1 Data. RNA-Seq quality control. This directory contains six subdirectories. The subdirectory "correlation" contains correlation heatmaps among all samples in diencephalon and telencephalon respectively. The subdirectory "mapping_stat" contains read mapping information on genomic features for each sample. The subdirectory "MDS" shows three dimensional MDS plots of the samples. The subdirectory "ReadDuplication" contains read duplication distributions for each sample. The subdirectory "ReadQuality" contains reads quality information for each sample plotted as both boxplots and heatmaps. The subdirectory "RPKMSaturation" contains information about read depth saturation for each sample as assessed by RPKM resamplings. All transcripts were divided into four quantiles based on their expression and a relative difference of observed and real RPKM values are plotted for each sample. (ZIP)