Social Regulation of Gene Expression in Threespine Sticklebacks

Identifying genes that are differentially expressed in response to social interactions is informative for understanding the molecular basis of social behavior. To address this question, we described changes in gene expression as a result of differences in the extent of social interactions. We housed threespine stickleback (Gasterosteus aculeatus) females in either group conditions or individually for one week, then measured levels of gene expression in three brain regions using RNA-sequencing. We found that numerous genes in the hindbrain/cerebellum had altered expression in response to group or individual housing. However, relatively few genes were differentially expressed in either the diencephalon or telencephalon. The list of genes upregulated in fish from social groups included many genes related to neural development and cell adhesion as well as genes with functions in sensory signaling, stress, and social and reproductive behavior. The list of genes expressed at higher levels in individually-housed fish included several genes previously identified as regulated by social interactions in other animals. The identified genes are interesting targets for future research on the molecular mechanisms of normal social interactions.


Introduction
Social interactions with conspecifics are found across all animal taxa, and the fundamental processes that govern social behavior are highly conserved. Among vertebrates, the core brain circuitry and key neuropeptides and neuromodulators that mediate social behavior are shared ( [1,2]; but see [3]). Furthermore, recent work has shown that gene networks that regulate social behavior are even conserved across invertebrates and mammals [4].
To identify genes and molecular pathways involved in social behavior, previous studies have examined animals with different social experiences to determine which genes show changes in expression [5][6][7]. These studies have either examined the expression of candidate genes or have employed expression arrays or transcriptome sequencing to more globally sample gene expression changes [5][6][7][8]. Global expression studies in vertebrates have identified numerous genes that are socially regulated, highlighting genes not previously associated with social behavior [4,[8][9][10][11][12][13][14][15]. These studies have been informative for dissecting the molecular mechanisms of sociality [4,5].
Here we sought to identify the genes that play a role in normal interactions among fish in a social group. We used threespine sticklebacks (Gasterosteus aculeatus), which are a longstanding model for studies of social behavior and have a wealth of genomic resources available, which facilitates transcriptomic analyses [16,17]. Marine sticklebacks are highly social, and are typically found in social groups [17,18]. We modulated the extent of social interactions of individual fish by housing fish either in social groups or individually for a one-week period. This manipulation should permit detection of a state change that is not due to the process of isolation (i.e. not within several hours), but also avoids the detrimental effects of long-term isolation on increasing stress and anxiety [19]. We then used RNA-sequencing (RNA-seq) to compare gene expression in brains of group-or individually-housed fish.

Fish and sample collection
Fish were from a lab-reared population of Japanese Pacific Ocean marine fish originally derived from the Bekanbeushi River in Japan. Fish were reared in 110-L tanks in 3.5 ppt seawater (Instant Ocean, United Pet Group, Blacksburg, VA) at 16 C, and under 16 h light / 8 h dark lighting conditions. Fish were fed Artemia nauplii and mysis shrimp. All fish were treated in accordance with the guidelines of the Institutional Animal Care and Use Committee of the Fred Hutchinson Cancer Research Center (FHCRC), protocol number 1575.
For social housing manipulation, fish from a single community tank were caught and transferred to four new 38-L tanks. Fish were either housed individually (n = 2 tanks) or in groups of eight mixed sex fish (n = 2 tanks). After one week of individual or group housing, we removed a single fish from each tank for analysis such that we had two individually-housed and two group-housed fish. We replicated this experiment with a second tank of fish so that we had a total of four biological replicates for both individually-and group-housed fish, from two original home tanks. Gonads were visually inspected to identify sex and maturity. Only prereproductive females were included in the experiment. Fish were euthanized with MS-222 and their brains were removed into RNA-later (Life Technologies, Carlsbad, CA) and stored at -20 C. Brains of individual fish were then dissected into three portions: 1) the telencephalon, 2) the diencephalon, pituitary, and rostral midbrain, and 3) the caudal midbrain, hindbrain, and cerebellum. We will refer to these portions as telencephalon, diencephalon, and hindbrain/cerebellum for simplicity. Tissue was homogenized using a pellet pestle (Kimble-Chase, Vineland, NJ) and total RNA was isolated using Trizol (Life Technologies, Carlsbad, CA). We performed the dissection and RNA isolation in separate batches on two different days, such that fish from one experimental replicate (i.e. home tank of origin) were processed on the same day.

RNA-seq
Barcoded RNA libraries from 24 samples (eight fish each with three brain regions) were generated in the FHCRC Genomics facility using Illumina's TruSeq RNA Sample Prep Kit v2 (Illumina, San Diego, CA) and a Sciclone NGS Workstation (PerkinElmer, Waltham, MA). Libraries were multiplexed, split across three lanes, and 50-bp paired-end sequences were generated on an Illumina HiSeq 2500 (Illumina, San Diego, CA). Demultiplexing was performed using Illumina's CASAVA v1.8.2 software, allowing for a single mismatch in the index read. Fastq files have been deposited to the Sequence Read Archive (Study Accession SRP056943). We used a local instance of Galaxy [20][21][22] to perform alignment and to quantify reads aligning to genes. Reads were first aligned to the stickleback genome (BroadS1 [16]) using the default parameters in tophat2 (version 2.0.9, Galaxy tool version 0.6 [23]). Next, reads that fell within predicted genes (Ensembl genes 76) were counted using htseq-count ("Count reads in features with htseq-count" Galaxy tool v1.0 [24]). In htseq-count, we used the following parameters:-q-m intersection-nonempty-s no-a 0-t exon-i gene_id. The resulting matrix was exported from Galaxy and imported into R (http://r-project.org) where we used edgeR, version 3.8.6, [25] to identify differentially expressed genes. A multidimensional scaling (MDS) plot was generated in edgeR. We also calculated the biological coefficient of variation (BCV) of samples using edgeR.
We first analyzed expression differences as a function of brain region, independent of social environment, by performing three analyses: telencephalon vs. diencephalon and hindbrain/cerebellum; diencephalon vs. telencephalon and hindbrain/cerebellum; and hindbrain/cerebellum vs. telencephalon and diencephalon. We filtered out genes that did not have at least 1 count per million reads in at least two samples. We present and discuss the top 10 differentially expressed genes for each brain region, all of which were significant at a False Discovery Rate (FDR) of P < 0.05.
To identify genes differentially expressed as a function of social environment, we next performed a General Linear Model (GLM) analysis separately for each brain region by comparing read counts in group-and individually-housed fish. We included experimental replicate (1 or 2) as a factor in the model to control for home tank of origin and RNA isolation-batch effects. We filtered out genes that did not have at least 1 count per million reads in at least two samples. Differentially expressed genes were those that had FDR of 0.05. We present and discuss the genes upregulated in group-and individually-housed fish separately, so for simplicity we report the log 2 fold change (log2FC) as positive for both comparisons.

Functional annotation and enrichment analysis
We used DAVID to perform functional annotation and enrichment analysis [26]. DAVID tests enrichment of Gene Ontology (GO) terms, as well as other annotation categories including Interpro domains, KEGG pathways, and SMART protein domains. Ensembl gene identifiers were first converted to zfin identifiers specifically for these analyses. Fold-enrichment of all significant up-or down-regulated genes was calculated over the background gene list, which included all genes expressed in the hindbrain/cerebellum. Functional annotation terms that were significantly enriched are reported, and are organized into clusters based on DAVID's functional annotation clustering.
We also tested for enrichment of glutamate receptors in genes upregulated in group-housed fish. We counted the number of glutamate receptor and GABA receptor genes in the upregulated list and the list of all genes expressed in the hindbrain/cerebellum. We then used the test of equal proportions in R to determine whether there was significant enrichment of these gene classes.

Results and Discussion
Sequencing generated an average of 42 ± 2 million total reads per sample, of which 88 ± 1% aligned to the genome. Of the aligned reads, 40 ± 2% fell within a predicted gene, thus were counted by htseq-count, and included in the analysis. Genes expressed at low levels were not included, leaving a total of 17,095 genes for the telencephalon, 17,553 for the diencephalon, and 17,081 for the hindbrain/cerebellum.

Differential expression as a function of brain region
We first compared gene expression as a function of brain region, independent of social housing condition. A multidimensional scaling plot showed clear separation of samples based on brain region (Fig 1). The top ten differentially expressed genes in each brain area based on log2FC included genes with known functions in these parts of the brain (Table 1). For example, the top ten genes enriched in the hindbrain/cerebellum were all known or predicted homeobox (Hox) transcription factors (Table 1). Hox genes are involved in hindbrain patterning during development and are expressed in the adult brain [27]. Eight of the top ten differentially expressed genes in the diencephalon encode pituitary hormones (Table 1), which was expected as this portion of the brain contained the pituitary. The other two diencephalon-enriched genes were Nr5a1a, which is expressed in the diencephalon of zebrafish [28] and Mibp, whose function in the brain has not been studied. In the telencephalon, the top ten differentially expressed genes included: 1) genes that are known to be involved in forebrain patterning and/or used as forebrain markers (Eomesa and Emx3 [29], Tbr1b [30], Scgn [31]), 2) a gene expressed in the forebrain of zebrafish (Rtn4rl2b [32]), and 3) genes with unclear function in the brain (Apod, Ctrb1). The genes identified as being highly enriched in specific brain regions may prove to be useful markers of different neuronal populations in future neuroanatomy studies in sticklebacks and other fish.

Differential expression as a function of social housing
We next identified genes that were differentially expressed as a result of social experience. There were numerous genes that were differentially expressed in the hindbrain/cerebellum (985 higher in group and 401 higher in isolate; all significant genes are shown in S1 File; the top 25 are shown in Tables 2 and 3). However, few genes were differentially expressed in either the diencephalon (5 higher in group) or telencephalon (1 higher in isolate). Four of the five differentially expressed genes in the diencephalon ( Table 2) were also upregulated in the hindbrain/ cerebellum of group-housed fish (S1 File; hindbrain/cerebellum values: Cyr61: log2FC = 2.4; FDR < 0.008, Tgm8: log2FC = 1.5; FDR = 0.008, Etv5a: log2FC = 0.9; FDR < 0.001, and Fam46d: log2FC = 1; FDR < 0.012). The fifth gene, novel gene ENSGACG00000012907, was not differentially expressed in the hindbrain/cerebellum (log2FC = 0.3; FDR = 0.13). Etv5a is a Multidimensional scaling plot shows leading log2 fold-change (log2FC) differences between samples. Brain regions are colored as follows: blue = telencephalon; red = diencephalon; black = hindbrain/cerebellum. Inset shows schematic of brain with the same colors representing dissected brain regions. Circles = group-housed samples; diamonds = individually-housed samples. transcription factor involved in specification of dopaminergic cells in C. elegans, and has been shown to co-localize with diencephalic dopaminergic cell populations in fish [33]. Cyr61 is expressed at the midbrain-hindbrain boundary in developing zebrafish, but its function is unknown [34]. Tgm8 shares distant homology with the transglutaminase family, which are enzymes involved in protein cross-linking [35]. Tgm8 was highly differentially expressed in all three brain regions, although it did not reach an FDR threshold of p < 0.05 in the telencephalon (higher in group; log2FC = 1.9; FDR = 0.13). Fam46d has an unknown neural function but is known to be expressed at higher levels in a mouse model of autism [36]. The single gene that was differentially expressed in the telencephalon is Proca1, whose function is unknown other than it is found in a protein complex with the cell division gene cyclin A1 (Table 3). It was interesting that many genes were differentially expressed in the hindbrain/cerebellum compared with few in either the telencephalon and diencephalon, which both contain nuclei known to be involved in the control of social behavior [37]. There are several possible explanations for this result. First, the hindbrain and cerebellum may indeed show a greater response to this alteration in social housing than the rest of the brain. Social interactions are associated with sensory stimulation, and this is reduced in individually-housed fish. The hindbrain serves as a primary sensory relay for several senses, and thus may show an increased transcriptional response to this manipulation. Alternatively, lack of detection of differentially expressed genes in the telencephalon and diencephalon could theoretically result from increased heterogeneity of these regions compared with the hindbrain/cerebellum. However, the coefficient of variation is similar across all brain regions (telencephalon BCV = 0.212; diencephalon BCV = 0.206; 2.3 0.008 S100u S100 calcium binding protein U All five significant genes from diencephalon and top 25 from hindbrain/cerebellum are shown; no genes were significantly upregulated in the telencephalon. Log2FC = log2 fold-change, FDR = false discovery rate, Symbol = gene name, NA = novel gene with no associated name. hindbrain/cerebellum BCV = 0.201), suggesting that this is not the cause in this case. Moreover, another study of stickleback gene expression differences that dissected the brain into similar portions did detect gene expression differences in all regions 30 min after social stimulation [14]. In that study, the diencephalon had the largest number of differentially expressed genes, whereas the telencephalon had the fewest. Thus, it is likely that there are differences in which brain regions respond to different stimuli. In addition, timing of stimulus exposure likely has an important impact on differential gene expression; this should be tested more thoroughly in future studies.

Genes upregulated in the hindbrain of group-housed fish
The 25 genes that were higher in the hindbrain of group-housed fish, based on fold-change, are shown in Table 2. Many of these genes were involved in developmental processes. The Hox genes and retinoic acid receptor (Rarb) are specifically involved in hindbrain development [38]. Several additional genes are otherwise implicated in neural development (Agrn [39] and Syne2a [40]) or intellectual disability (Herc2 [41] and Kcnk9 [42]), and Stard9 is involved in cell division [43]. Functional annotation and enrichment analysis echoed the finding that developmental genes are strongly enriched in the list of genes upregulated in group-housed fish ( Table 4). All of the significantly enriched functional clusters were related to development, including cell morphogenesis and neural development, cell adhesion, plexin/semaphorin signaling, and EGF signaling (Table 4). Semaphorins and EGF signaling are involved in neural development [44,45]. Increased activity of developmental processes is suggestive of more arborization and neurogenesis in group-housed fish. There is ongoing neurogenesis in the hindbrain/cerebellum of sticklebacks [46], and previous work has shown that sensory stimulation, including social housing, can alter levels of neurogenesis in other fish [47]. It is possible that the upregulated gene expression of developmental genes in the hindbrain/cerebellum of group-housed fish is related to increased sensory function due to higher levels of sensory stimulation.
Other genes in the top 25 upregulated genes included Szt2 and Sst1.1. Szt2 mutant mice have a lower seizure threshold [48]. Somatostatin (Sst1.1) has previously been implicated in decreasing growth as well as decreasing aggressive behavior in fish [49,50]. Social isolation can Table 4. Functional annotation and clustering of genes expressed at higher levels in group-housed fish.

Cluster
Term Description Fold Enrichment lead to increased aggression in fish [51]. It would be interesting to determine whether group housed sticklebacks have slower growth and reduced aggression than individually-housed fish. In addition, it would be interesting to manipulate somatostatin levels [49] and determine whether there was an impact on growth and gene expression. The list of 985 genes upregulated in the hindbrain/cerebellum as a result of group housing included many other interesting genes in addition to those presented in Table 2. We will highlight a few here, although the entire list can be found in S1 File. Many enriched genes were in neurotransmitter or neuromodulator pathways. First, several genes related to acetylcholine synthesis and signaling were higher in group-housed fish: acetylcholinesterase (Ache, ENSGACG00000000728; log2FC = 0.9; FDR = 0.009), choline o-acetyltransferase (Chat, ENSGACG00000002482; log2FC = 0.8; FDR = 0.008), and the muscarinic acetylcholine receptor, Chrm2a (ENSGACG00000019948; log2FC = 1.2; FDR = 0.04) (S1 File). Acetylcholinergic cells are found in cranial sensory and motor nuclei and throughout the reticular formation of the hindbrain [52]. Chrm2a also expressed in cranial nuclei [53]. Given these expression patterns, we speculate that increased acetylcholine signaling is related to higher levels of sensory processing due to more sensory stimulation in the group-housing environment.
Another signaling pathway gene that was differentially expressed was prostaglandin F2 receptor inhibitor (Ptgfrn; ENSGACG00000014419; log2FC = 0.8; FDR = 0.013). Prostaglandin F2α signaling increases fish reproductive physiology [57] and behavior [58]. The females in social groups were exposed to males but isolated females were not, so it may be that mixed-sex housing facilitates reproduction. Investigating levels of reproductive hormones would directly address this question.
Opiate signaling pathway genes were also regulated as a function of social status. Prepronociceptin a (Pnoca, ENSGACG00000014805; log2FC = 1.7; FDR = 0.003) and its receptor, opiate receptor-like 1 (Oprl1; ENSGACG00000010479; log2FC = 1.1; FDR = 0.02), were both expressed at higher levels in fish in social housing. Interestingly, Pnoc and Oprl1 (aka NOP) were also found to be higher in brains of mice housed in groups than in mice housed in isolation [59]. Nociceptin signaling decreases stress and anxiety in mammals [60]. It may be that social interactions in group-housed fish lead to increased nociceptin signaling, which results in reduced stress and anxiety. Alternatively, individually-housed fish might have decreased levels of nociception signaling.
Finally, 14 glutamate receptor subtypes were found in the list of significantly upregulated genes in socially housed fish (S1 File; Gria1a, Gria4b, Grik2, Grik3, Grik5, Grin2ab, Grin2b, Grin2bb, Grin2ca, Grin2db, Grip2b, Grm3, Grm5, Grm8). Because the glutamate receptor family is quite large, we tested to see whether this was a specific enrichment or was simply a result of there being a large number of glutamate receptor genes in the entire gene list. We also compared the level of enrichment of another large neurotransmitter receptor family, the GABA receptors. This analysis showed that glutamate but not GABA receptors were significantly enriched in fish housed in social groups (Χ 2 = 43, P < 0.00001; Fig 2).

Genes upregulated in the hindbrain of individually-housed fish
We next examined genes that were higher in the hindbrain of individually-housed fish ( Table 3). The list of the top 25 genes with the highest fold-change contained genes with diverse functions. For example, Slc16a1 has been implicated in neurogenesis in zebrafish [61]. Ddit4 may play a role in development through interactions with Wnt/beta catenin signaling [62]. There were several transcription factors with varied functions (Tcf24, Mxd3, Sox1b). Gatm is involved in creatine synthesis. Interestingly, Mad2l1bp, which has homology to a gene involved in cell division and the spindle checkpoint pathway, was also found to be regulated by social interactions in other populations of sticklebacks. Specifically, it was higher in males following a territorial intrusion [14]. Finally, novel gene ENSGACG00000001910 has homology to the MHC class 1 antigen family. A gene from this family was previously shown to be expressed at higher levels in brains of female than male cichlids [11].
The entire list of 401 genes upregulated in individually-housed fish included several other genes with interesting functions, and is shown in S1 File. One of these was an enzyme involved in steroid biosynthesis, hydroxysteroid (17-beta) dehydrogenase 7, which was expressed at higher levels (Hsd17b7; ENSGACG00000016134; log2FC = 0.6; FDR = 0.02). Hsd17b7 is involved in the biosynthesis of cholesterol and sex steroids, and thus may play a role in regulating steroid hormone abundance in the brain. Another gene upregulated in individually-housed Socially-Regulated Genes in Sticklebacks fish, MAD2 mitotic arrest deficient-like 1, was also shown to be higher in brains of isolated rats (Mad2l1; ENSGACG00000001594; log2FC = 0.8; FDR = 0.007) [13].
We next performed functional annotation and enrichment analysis of the list of genes upregulated in individually-housed fish. Relatively few categories were enriched, and they included genes related to RNA processing (Table 5).

Conclusions
In summary, we found that manipulating social housing impacted the expression of genes predominantly in the hindbrain/cerebellum. In group-housed fish, many of the upregulated genes were in developmental signaling pathways, and functional annotation reinforced the conclusion that there was enrichment of development-related genes in this dataset. These results suggest that fish in group-housing environments experience more neurogenesis or more axon and dendrite outgrowth. Alternatively, because many developmental genes act as repressors, it may be that upregulated expression of these genes is actually associated with decreased neurogenesis. It would be interesting to distinguish between these possibilities by directly by comparing levels of cell division and differentiation on a cellular level. Other differentially expressed genes were involved in stress/anxiety, social behavior, and possibly sensory processing. These findings suggest interesting directions for future research on the molecular control of normal social interactions in sticklebacks and other systems. In the future it could also be interesting to evaluate different timescales of experimental manipulation, for instance social isolation for an entire lifetime or across evolutionary timescales [63].
Supporting Information S1 File. List of all differentially expressed genes in the hindbrain. File contains a list of all hindbrain genes that were significantly upregulated (FDR < 0.05) in group-and individuallyhoused fish, on two separate worksheets. (XLSX)

Acknowledgments
We thank Shaugnessy McCann for fish care; Jeff Delrow, Andy Marty, Ryan Basom and the FHCRC Genomics Shared Resource for RNA-seq library preparation and sequencing; and Brian Claywell for help with local Galaxy setup.

Author Contributions
Conceived and designed the experiments: AKG CLP. Performed the experiments: AKG. Analyzed the data: AKG. Wrote the paper: AKG CLP. Socially-Regulated Genes in Sticklebacks