Weisner PA, Chen CY, Sun Y, Yoo J, Kao WC, Zhang H, Baltz ET, Troy JM, Stubbs L
Abstract
AUTS2 was originally discovered as the gene disrupted by a translocation in human twins with Autism spectrum disorder, intellectual disability, and epilepsy. Since that initial finding, AUTS2-linked mutations and variants have been associated with a very broad array of neuropsychiatric disorders, suggesting that AUTS2 is required for fundamental steps of neurodevelopment. However, genotype-phenotype correlations in this region are complicated, because most mutations could also involve neighboring genes. Of particular interest is the nearest downstream neighbor of AUTS2,GALNT17, which encodes a brain-expressed N-acetylgalactosaminyltransferase of unknown brain function. Here we describe a mouse (Mus musculus) mutation, T(5G2;8A1)GSO (abbreviated 16Gso), a reciprocal translocation that breaks between Auts2 and Galnt17 and dysregulates both genes. Despite this complex regulatory effect, 16Gso homozygotes model certain human AUTS2-linked phenotypes very well. In addition to abnormalities in growth, craniofacial structure, learning and memory, and behavior, 16Gso homozygotes display distinct pathologies of the cerebellum and hippocampus that are similar to those associated with autism and other types of AUTS2-linked neurological disease. Analyzing mutant cerebellar and hippocampal transcriptomes to explain this pathology, we identified disturbances in pathways related to neuron and synapse maturation, neurotransmitter signaling, and cellular stress, suggesting possible cellular mechanisms. These pathways, coupled with the translocation’s selective effects on Auts2 isoforms and coordinated dysregulation of Galnt17, suggest novel hypotheses regarding the etiology of the human “AUTS2 syndrome” and the wide array of neurodevelopmental disorders linked to variance in this genomic region.
Read more at G3
![Published: July 13, 2017https://doi.org/10.1371/journal.pgen.1006840 Article Authors Metrics Comments Related Content Abstract Author summary Introduction Results Discussion Methods Supporting information Acknowledgments Author Contributions References Reader Comments (0) Media Coverage (1) Figures Abstract 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. Author summary Social interactions provoke changes in the brain and behavior but their underlying molecular mechanisms remain obscure. Male sticklebacks are small fish whose fitness depends on their ability to defend a territory. Here, by measuring the time course of gene expression in response to a territorial challenge in two brain regions, we show that a single brief territorial intrusion provoked waves of gene expression that persisted for hours afterwards, with waves of transcription associated with distinct biological processes. Moreover, a single territorial challenge caused dramatic changes to the epigenome. Changes in chromatin accessibility corresponded to changes in gene expression, and to the activity of transcription factors operating within gene regulatory networks. This study reveals rapid and dramatic epigenomic plasticity in response to a brief, highly consequential social interaction. These results suggest 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. Figures Fig 4Table 1Fig 1Fig 2Fig 3Fig 4Table 1Fig 1Fig 2Fig 3 Citation: Bukhari SA, Saul MC, Seward CH, Zhang H, Bensky M, James N, et al. (2017) Temporal dynamics of neurogenomic plasticity in response to social interactions in male threespined sticklebacks. PLoS Genet13(7): e1006840. https://doi.org/10.1371/journal.pgen.1006840 Editor: David Clayton, Queen Mary University of London, UNITED KINGDOM Received: December 23, 2016; Accepted: May 27, 2017; Published: July 13, 2017 Copyright: © 2017 Bukhari et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: The RNA-Seq data and ChIP-seq data are accessible on GEO (https://www.ncbi.nlm.nih.gov/geo/) with this ID: GSE96673 Funding: This work was supported by grant #SFLife 291812 from the Simons Foundation to LS and Gene E. Robinson (https://www.simonsfoundation.org/), by National Science Foundation grant # 1121980 to AMB and Katie McGhee (www.nsf.gov) and the National Institutes of Health grant # 2R01GM082937 to AMB (www.NIH.gov). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist. Introduction Animals exhibit remarkable behavioral plasticity. Social interactions in particular can provoke moment-to-moment changes in behavior. These changes are coordinated at the neural level, but social interactions also elicit transcriptional changes within the brains of behaving animals [1]. For example, genome-wide transcription studies show that roughly ~10% of the genome responds to a mating opportunity [2–7], predation risk [8–10], or a territorial challenge [11–13]. 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. Results 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). thumbnail Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. Brain region-specific changes in gene expression in response to a territorial challenge over time. (A) Numbers of up- (blue) 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 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. https://doi.org/10.1371/journal.pgen.1006840.g001 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), 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 DEGx (FDR < 0.1, S4 Table). We next used hierarchical clustering of the DEGx 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. thumbnail Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Hierarchical clustering of genes whose expression profiles changed over time in response to a territorial challenge (DEGx) and their functional enrichments. Hierarchical clustering grouped together genes with similar expression profiles over time. 13 clusters were identified in diencephalon (D1-D13, A). 12 clusters were identified in telencephalon (T1-T12, B). Each line represents the expression pattern of a different gene, where positive fold change indicates upregulation and negative fold change indicates downregulation in response to a territorial challenge. Clusters of genes with similar expression profiles (columns) had different GO molecular functions associated with them (rows); C) diencephalon; D) telencephalon. Some clusters did not have significant functional enrichment. https://doi.org/10.1371/journal.pgen.1006840.g002 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 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 DEGx 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). thumbnail Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Network of interacting transcription factors (TFs) in the transcriptional regulatory network highlighting enrichments of TFs in clusters of DEGx. 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. https://doi.org/10.1371/journal.pgen.1006840.g003 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 DEGx 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 (DEGx, Fig 4C). The DEGx associated with nearby DAPs (hereafter referred to as DAPDEGx) 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). thumbnail Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Connecting gene expression and chromatin accessibility in diencephalon. (A) Fold change of differentially accessible peaks at 30 minutes and 120 minutes; blue indicates up in experimental, red indicates down in experimental. (B) Functional enrichment (molecular function) of genes associated with differentially accessible peaks at 30 minutes and 120 minutes. Blue indicates up in challenged, red indicates down in challenged. (C) Overlap of genes whose expression profile changed over time in response to a social interaction (DEGx) with genes associated with differentially accessible peaks at 30 minutes and 120 minutes. The overlap between DEGx and accessibility at 120 minutes is statistically significant (P 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 makeBigWig.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 DEGx 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 genome-scale 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 (PANTHER 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](https://i0.wp.com/stubbs.igb.illinois.edu/wp-content/uploads/2017/11/journal.pgen_.1006840.g003.png?resize=750%2C499&ssl=1)
