Fish and chips: functional genomics of social plasticity in an African cichlid fish

Brain and behavior are sculpted by a dynamic interplay between genotype and environment. Changes in social status mediate changes in a wide range of molecular, neuroanatomical, endocrine and behavioral characteristics(Oliveira et al., 2002; Fernald, 2004; Sapolsky, 2005) at different time scales (Hofmann, 2003). Nervous systems respond by modulating neural properties and neuroendocrine responses. Such rapid, short-term changes will often lead to sustained,longer-term changes in the brain and behavior via differential gene expression, leading to subsequent structural and physiological changes. However, few studies have attempted to determine the molecular underpinnings of socially regulated phenotypes within a comprehensive genomic framework,especially in vertebrates. In the present study, we apply a molecular systems approach in a highly social fish species to ask how social context regulates the molecular and physiological substrates that, in turn, sculpt subsequent behavior.

The African cichlid fish, Astatotilapia burtoni (Taxonomic Authority) (formerly Haplochromis burtoni) has become an important model system to study the mechanisms underlying socially mediated behavioral change. In this species, 20–30% of males are dominant (D), slow growing,brightly coloured and actively defend territories for mating. The remaining subordinate (S) males mimic females by schooling and displaying cryptic coloration, while experiencing faster growth(Fernald, 1977; Hofmann et al., 1999a). Subordinate males show little aggression and territoriality and, importantly,have regressed gonads and are thus not reproductive(Fernald and Hirata, 1977a; Fernald and Hirata, 1977b; Francis et al., 1993). These behavioral and physiological characteristics are plastic and influenced by the immediate social environment, such that an individual male switches between the D and S phenotypes several times during its life depending upon its relative ability to obtain and maintain access to a territory through encounters with other males (Hofmann et al., 1999a). Environmental conditions, availability of territorial shelters, relative body size and physiological condition influence the probability of acquiring and maintaining a territory. The phenotypic switch occurs over a timescale of minutes to days to weeks in both the field and the laboratory (White et al.,2002; Hofmann,2003; Burmeister et al.,2005).

In the laboratory, A. burtoni has been the focus of hormonal and molecular studies related to a broad range of phenotypic traits that are affected by the transition between the two male phenotypes(Fig. 1) (for reviews, see Hofmann and Fernald, 2000; Fernald, 2002; Hofmann, 2003; Fernald, 2004). Variation in several different components of the endocrine system reflects the complexity underlying the corresponding phenotypic switch between these dramatically different male phenotypes. Neuroendocrine pathways regulating androgen production (Parikh et al.,2006), growth (Hofmann et al.,1999a; Hofmann et al.,1999b) and stress response(Fox et al., 1997) change in a complex fashion as animals undergo phenotypic change(Fig. 1). Parikh et al.(Parikh et al., 2006)suggested that the higher levels of testosterone (T) and 11–ketotestosterone (11–KT) that they measured in D males might promote aggressive behavior, as has been shown using androgen manipulation in other fish species. Production and release of neuropeptides and neuromodulators such as gonadotropin-releasing hormone (GnRH1)(White et al., 2002) and somatostatin (Hofmann and Fernald,2000; Trainor and Hofmann,2006; Trainor and Hofmann,2007) are higher in dominant males. According to White et al.,GnRH1 mRNA levels and gonadosomatic index (GSI) (a measure of gonadal development) are positively correlated(White et al., 2002). The higher level of GnRH1 in D males coincides with territory acquisition and mating opportunity; as such, it is similar to GnRH changes in seasonally reproducing species (Amano et al.,1995; Dawson et al.,2001; Nelson,2005; Hofmann,2006). Fluctuations at the molecular level also co-vary with the observed behavioral switch of male phenotype. In the case of GnRH(Parhar et al., 2005; Au et al., 2006) and somatostatin (Trainor and Hofmann,2006), specific receptor sub-types are regulated according to social status. Furthermore, the membrane properties of GnRH1-expressing neurons, reduced excitability in S males correspond with the differences in peptide release between D and S males(Greenwood and Fernald,2004).

Manipulation of the social environment allows experimental control of the phenotypic switch (Francis et al.,1993; White et al.,2002; Burmeister et al.,2005). The ease of experimental manipulation is paired with a wealth of ecological and evolutionary information available for haplochromine cichlids in general (for a review, see Kocher, 2004; Salzburger et al., 2005). For example, there is thought to be a trade-off between reproduction and survival,such that territorial males, the sole reproducers, are also subject to higher mortality rates through predation, probably due to their conspicuous coloration (Fernald and Hirata,1977a; Fernald and Hirata,1977b; Maan et al.,2008). Furthermore, A. burtoni is currently undergoing whole-genome sequencing, along with three other African cichlid species (see www.broad.mit.edu/models/tilapia/),and thus offers an unrivaled laboratory-based model system for the genomic analysis of complex and ecologically relevant phenotypes.

cDNA microarrays for transcript profiling, have become a powerful tool when applied to species of behavioral, ecological or evolutionary interest [e.g. alternative life histories (Aubin-Horth et al., 2005), cooperative breeding(Aubin-Horth et al., 2007),social behavior (Grozinger and Robinson,2002; Whitfield et al.,2003; Robinson et al.,2005), response to heat stress(Buckley et al., 2006),response to environmental estrogens(Martyniuk et al., 2006) and physiology of drug addiction (Rhodes and Crabbe, 2005)]. In the present study, we employ a microarray platform that contains many known candidate genes as well as ∼4000 brain-derived cichlid cDNAs (Renn et al.,2004). Such a combined candidate gene and genomic strategy allows hypothesis-driven and discovery-based experiments on a single platform. The obviously complex nature of behavioral traits – such as this socially regulated, reversible switch between D and S phenotypes – requires a discovery-based approach in order to identify the many genes involved.

In the present study, we analyze gene transcript patterns for reproductively active D males and reproductively suppressed S males as well as for brooding females. We then compare the expression patterns for each of the three phenotypes in order to identify gene sets associated with reproduction or dominance behavior, providing insight into the molecular modularity underlying these phenotypes. Next, we annotate the array features according to gene ontology (GO), with the goal of identifying gene regulation within molecular categories free of a priori expectations or experimenter bias. Finally, we examine variation in gene expression patterns between individual animals within a social phenotype and ask whether any of these variable genes are also those that are differentially expressed between the social phenotypes.

Astatotilapia burtoni (Günter) from a wild-caught stock population were kept in aquaria under conditions mimicking their natural environment (Fernald and Hirata,1977b): pH 8.0, 28°C water temperature and 12 h:12 h light:dark cycle. Gravel substrate and terracotta shelters allowed the establishment and maintenance of territories necessary for reproduction(Fernald and Hirata, 1977a). Fish were fed every morning with cichlid pellets and flakes. All work was carried out in compliance with Animal Care and Use Guidelines.

Males were marked using colored bead combinations attached near to the dorsal fin. Nine groups of 2–3 males with 2–3 gravid females were established in half of a 100 liters aquarium. Each group was visually isolated from neighboring fish. Ten-minute behavioral observations were made approximately twice every week for five weeks by directly observing chasing,threat, display, border threat, courtship, flee, schooling and carousel (as described by Fernald, 1977). Behavioral measures were used to calculate a dominance index (DI), i.e. the sum of all aggressive behaviors observed (threat, chasing, border threat and carousel) minus the number of submissive behaviors observed (flee) and a reproductive index (RI), i.e. the total of all reproductive behaviors (court,dig, spawn). On the final day of the experiment, males of each status were taken for gene expression profiling only if they had continuously expressed either the D or S phenotype for the past 28 days. Six mouth-brooding females were also selected for expression profiling. The 18 fish used for gene expression analysis were obtained from a total of seven different groups on the basis of stable phenotype. For the purpose of hybridization design the animals were treated as equivalent within phenotype, though comparisons between individuals can be made indirectly (see below). Close inspection of the behavioral and microarray data did not reveal any covariation between animals derived from the same social group. Standard length, body mass and gonad mass were measured for each fish in order to calculate the GSI as gonad mass/body mass. Condition factor (CF) was calculated based on the residuals from the regression of body mass on standard length for before and after the experiment (r²=0.95; P=2.79×10^–6), and growth rate (GR) was calculated as the relative change in standard length over the course of the experiment. Whole brains were dissected and stored in RNAlater (Ambion,Austin, TX, USA) within 5 min of initial tank disturbance.

The array used in the present study has numerous redundant features, i.e. two or more features represent the same gene. We exploited this property for quality control purposes and to assess the sensitivity of the approach. For example, the array includes four independent features that represent the neuropeptide GnRH1 (two previously cloned cDNAs and two obtained independently from the cDNA library), known to be expressed in only ∼300 neurons and upregulated in D males (Davis and Fernald,1990; White et al.,2002). All four features were similarly significantly differentially expressed, being upregulated with dominance, and thus serve as both a biological and a technical control demonstrating the sensitivity even with whole-brain RNA.

Brains were homogenized (Tissue tearor, Biospec products, Bartlesville, OK,USA) and total RNA was extracted according to standard Trizol (Invitrogen,Carlsbad, CA, USA) and phase-lock gel (Eppendorf; Westbury, NY, USA)protocols. RNA integrity was determined on the Bioanalyzer (Agilent; Santa Clara, CA, USA) and spectrophotometer (Agilent) prior to indirect RNA labeling protocol, starting with 2 μg of total RNA according to Renn et al.(Renn et al., 2004). Briefly,each RNA sample was labeled twice, once with Cy3 and once with Cy5. After purification from unincorporated label, each sample was divided in two,combined with the appropriate samples and every individual was compared with two individuals from each of the other two phenotypes in a balanced loop design incorporating a dye-reversal (see below for how to access annotation of submitted data at NCBI's GEOdatabase for further details). Targets were hybridized to the brain-specific cDNA array from A. burtoni(Renn et al., 2004) at 65°C for 12–16 h. Arrays were scanned on an Axon 4000B arrays scanner (Genepix 4.0; Molecular Devices, Sunnyvale, CA, USA). The loop design allows for the direct comparison of samples of interest, thus offering greater statistical power with fewer replicates(Churchill, 2002).

After filtering for bad feature morphology, hybridization artifacts and low intensity (<2 s.d. above local background), raw data were imported into R software [v.1.0 R-Development team, 2004, Vienna, Austria(R Development Core Team 2006)] and normalized using the Linear Models for Microarray Data package [LIMMA v.1.6.6 (Smyth et al.,2004)]. Background-subtracted mean intensities were normalized using a within-array printtip-lowess normalization and used to calculate ratios for a Bayesian analysis of gene expression levels [BAGEL v.3.6(Townsend and Hartl, 2002)]. BAGEL takes advantage of additional information obtained from transitive comparisons of individuals in loop designs experiments(Townsend and Hartl, 2002; Churchill, 2002). Genes represented by more than one feature on the array were only counted as significant if at least one representative passed the significance threshold and the full complement survived Fisher's test of combining probabilities from multiple tests of significance [p.794 in Sokal and Rohlf(Sokal and Rohlf, 1995]. All raw and processed data are available at GEOdatabase(www.ncbi.nlm.nih.gov/projects/geo/),sample numbers GSM267785–819 of series GSE10624.

The clone templates for PCR amplification were end sequenced(Salzburger et al., 2008),resulting in 4258 expressed sequence tags (ESTs), 3670 of which were deemed to be of high quality and have been submitted to GenBank (accession number CN472211–CN468542) and are maintained at the Dana-Farber Cancer Institute as a GeneIndex(http://compbio.dfci.harvard.edu/tgi/cgi-bin/tgi/gimain.pl?gudb=a_burtoni)such that 1280 clones combine into 399 tentative contigs (TC), leaving 2381 singleton sequences for a total of 2780 unique sequences. Half (49.5%) of the sequences (1408 out of 2842) could be annotated with a `best hit' at a threshold of e<sup>–12</sup> or better according to BLAST alignment to UniProt (v.5.2). GO terms were applied to 869 unique cichlid genes by transitive annotation, meaning that the GO annotations for a cichlid gene's best hit were collected and used in the present study for further analysis. In order to avoid species bias, we collected GO terms from all genes with the same name as the best hit annotation. GO annotations also include confidence codes. The less reliable annotations derive from `inferred sequence similarity(ISS)' or `inferred electronic annotation (IEA)'(The Gene Ontology Consortium,2000), therefore we excluded those annotations when transitively applying GO terms to the ESTs represented on the A. burtonimicroarray. The resulting GO graphs [referred to as directed acyclic graphs(DAGs)] were then `slimmed' to 183 terms and a total of 4102 total annotations. For the `slimming' process, the leaf-most nodes that were selected to contain a minimum of 10 annotated cichlid sequences, and parent nodes were retained only if an additional 10 cichlid sequences were annotated at that level.

Over- and under-representation of GO terms for a regulated set of genes was determined in Cytoscape (Shannon et al.,2003) using the Biological Network Gene Ontology tool, BiNGO(Maere et al., 2005), which relies upon hypergeometric statistical significance. As GO categories are highly non-independent, the statistical treatment of these terms is still under discussion (Ge et al.,2003). Also, owing to the small number of genes for each ontology term and the relatively small number of genes that are regulated, there is less statistical power to identify significantly under-represented GO terms. For these reasons, we use GO analysis as a hypothesis-generating tool and report only uncorrected hypergeometric P-values.

Prior to clustering, features representing replicate ESTs were collapsed by combining probabilities from multiple tests of significance [p. 794 in Sokal and Rohlf (Sokal and Rohlf,1995)], and the mean expression level was determined for each set of features. A hierarchical clustering analysis was applied to the list of genes that were significantly regulated according to D, S and brooding female phenotypes. The estimated gene expression levels were used to obtain the dissimilarity matrix by applying Euclidian distance measure, which integrates the effects of amplitude of ratios as well as direction (correlation) in patterns. Clustering analysis of gene expression patterns of each individual was performed using the hclust function in R software v.2.0.1. Clustering was based on dissimilarity measures obtained using the dist functions in the stats package. The consensus tree and bootstrap confidence values for each tree node were obtained with the consensus function in the maanova package(Wu et al., 2002). The consensus tree and confidence values were calculated as the proportion of trees obtained with bootstrapped datasets that agreed with the original tree. Each bootstrapped tree was based on the Euclidian distance matrix calculated for each of 1000 permutated gene expression profile datasets obtained by resampling with replacement. Alternate clustering methods and different measures of distance are available and are similarly appropriate for gene expression analysis. Hierarchical clustering, based upon all features on the full array rather than on regulated genes only, provided a similar tree with reduced confidence values at each node (not shown). The heatmapfunction and colour options in the package gplots were used to visualize clusters of gene expression, the z-transformed expression ratios were grouped by k-means function in the stats package of R and ordered as such, while the samples were ordered according to the consensus hierarchical cluster. GO terms provide a means to address the possible functional relationship of a cluster of genes that are coordinately regulated. However, no statistically significant over- or under-representation of GO terms was seen for any of the gene groups identified by the k-means clusters according to a hypergeometric test (not shown).

Individually tagged males were placed into community tanks with females such that naturalistic dominance hierarchies were established among the males according to standard protocol (Hofmann et al., 1999a). Regular observations over the course of five weeks identified males of definitive dominant (D) and subordinate (S) phenotypes. The territory-holding males (N=6) behaved more aggressively than S males (N=6), as reflected by their significantly higher DI (D males,15.7±1.7; S males, –3.5±0.8; t-test, P<0.001) and RI (D males, 2.12±0.57; S males,0.0±0.0; t-test, P=0.004). The GSI varied continuously across the two groups, showing the expected (yet, in this study,non-significant) trend for higher GSI in dominant males (D males,0.85±0.23; S males, 0.56±0.04; t-test, P=0.25). There was no significant difference in growth rate (GR) over the course of the experiment (D males, 1.25±0.215; S males,1.39±0.285; t-test, P=0.69) or CF (residuals) at the time of sacrifice (D males, 0.01±0.053; S males,–0.01±0.061; t-test, P=0.58). Note, however,that CF differed marginally when the social communities were established at the beginning of the experiment (D males, 0.03±0.051; S males,–0.03±0.068; t-test, P=0.07), which could suggest that increased condition may be advantageous when ascending to dominance in the first place. Specific behaviors also varied significantly and predictably between the D and S males (Fig. 2): increased presence of an eye-bar in D (P<0.001);increased schooling behavior of S (P=0.0005); more fleeing behavior by S (P=0.004); and more frequent chasing behavior(P=0.0003), border threat behavior (P=0.08) and sexual behavior (P=0.003) by D. These results are representative of D and S males in A. burtoni as previously described (e.g. Francis et al., 1993; Hofmann and Fernald, 2000; White et al., 2002).

For the first investigation of our microarray data, we conducted Bayesian analysis of gene expression (Townsend and Hartl, 2002) by treating the six individuals of each phenotype as biological replicates and calculated gene expression differences at the phenotype level (D, S or brooding female). When comparing the different phenotypes, array features representing 3647 cichlid unique sequences passed quality filters in a sufficient number of hybridizations to be considered in the results. D and S phenotypes showed significant differences in gene expression for 171 genes [Bayesian posterior probability, BPP>0.99) (87 upregulated; 84 downregulated with dominance (supplementary material Tables S1 and S2)]. Permutation analysis of the data shows that, at this significance threshold of BPP>0.99, only ∼17% of these 171 genes could be expected to show significant differential expression by chance, a reasonable rate of false positives (Table 1). Therefore, almost five percent of the genes studied were differentially regulated in the brain according to male dominance phenotype. This percentage is considerably smaller than that found in typical honey bee colonies (39%),where nurse and forager phenotypes are not only distinguished by their social role but also differ in their daily behavioral routines and the environment in which they move about (Whitfield et al.,2003). These results also show phenotype-specific regulation of a considerably smaller portion of the genome compared with the 15% observed difference between alternative mating tactics in a study of Atlantic salmon males, which similarly compared sexually mature and immature males(Aubin-Horth et al., 2005).

The microarray was designed to include candidate genes previously studied in the context of social dominance and other behavioral contexts in A. burtoni (Hofmann, 2003; Fernald, 2004). The inclusion of known candidate genes allowed us to test multiple hypotheses and also offered validation of the microarray results by comparison with previous studies for some of these genes (summarized in Fig. 1). For instance,peptidergic neurons in the pre-optic area (POA) and other brain regions express several neurohormones (e.g. arginine vasotocin, GnRH1, galanin) as well as neurohormone receptors and steroid receptors, which have previously been shown in separate studies to play a role in the regulation of social behavior. Below, for the first time, we provide a combined analysis of these neuroendocrine pathways in A. burtoni.

As predicted from previous studies using ribonuclease protection assays and in situ hybridization, among the three GnRH neuropeptide genes that are expressed in the brain of fish, only GnRH1, the form expressed in the POA(White et al., 1994; White et al., 2002), showed highly significant regulation in the microarray results (BPP=0.9998), with D males having higher levels. Given the small number of cells expressing this neuropeptide (∼300) (Soma et al.,1996; Munz, 1999),confirmation of GnRH1 regulation by our microarray analysis provides an important cross-validation and confirms the sensitivity of the array even when using whole-brain RNA. As predicted from previous studies(White et al., 1994; White et al., 2002), our results also confirmed that the other two forms of GnRH, GnRH2, (midbrain) and GnRH3 (terminal nerve), are not regulated according to male social phenotype(BPP=0.366 and 0.700, respectively). None of the GnRH receptor sub-types on the array were significantly regulated, although studies have demonstrated their regulation in the pituitary relative to sexual maturity and social status (Parhar et al., 2005; Au et al., 2006).

Galanin, a neuropeptide that links metabolic activity and reproduction through regulation of GnRH release (reviewed by Kageyama et al., 2005; Tortorella et al., 2007), was marginally upregulated in D males (BPP=0.9592). There is considerable evidence from mammals that galanin reduces nociception(Wiesenfeldhallin et al.,1992), increases food intake(Schick et al., 1993) and stress reactivity (Holmes et al.,2002), plays a role in regulation of sexual behavior, and is itself regulated by GnRH and estrogen(Gabriel et al., 1993). Specifically, in fish, galanin is thought to play a role in regulation of food intake and is widely distributed in the brain, being localized to the olfactory bulb, telencephalon, hypothalamus, midbrain and posterior brain(reviewed by Volkoff et al.,2005), as well as the pituitary(Jadhao and Pinelli, 2001) and peripheral tissues (Johnsson et al.,2001). The results of the present study suggest the intriguing possibility that galanin might be upregulated in D males as a response to reduced food intake and constant challenges to their social status by other individuals. Future studies will test these novel hypotheses.

Arginine vasotocin (AVT; represented by multiple clones on the array), the non-mammalian homolog of arginine vasopressin (AVP), was among the most strongly regulated genes in this study, being upregulated in the brains of D males (BPP>0.9999). AVP/AVT has been implicated in the regulation of social behavior across vertebrates, including aggression and social affiliation(Goodson, 1998; Goodson and Adkins-Regan,1999; Winslow et al.,1993). In teleost fish, AVT is known to play a role in male mating tactics [peacock blenny (Grober et al.,2002; Carneiro et al.,2003); midshipman (Goodson and Bass, 2001)], as well as in the behavioral regulation of sex change and the associated territory acquisition [bluehead wrasse(Semsar and Godwin, 2003; Semsar and Godwin, 2004)]. AVT is also associated with territorial aggression [damselfish(Santangelo and Bass, 2006)]as well as dominant and territorial behavior in both the male and female of a breeding pair when compared with their subordinate helpers in another cichlid species, the cooperative breeding Neolamprologus pulcher(Aubin-Horth et al., 2007). Our data from A. burtoni suggest a role for AVT in regulation of dominance and are consistent with an in situ hybridization-based study (Greenwood et al.,2008). The AVT V1a receptor, which plays a fundamental role in affiliative behaviors in voles (Lim et al., 2004), was not represented on the array.

The enzyme aromatase, which converts testosterone to estrogen, is important in sex determination (Nakamura and Kobayashi, 2005), sex change in fish(Black et al., 2005; Marsh et al., 2006) and regulation of social behavior (Hallgren et al., 2006). There are two isoforms of aromatase, one localized to the brain and the other to the gonads, and both are represented on the microarray. In the present study, D males showed increased neural expression of the brain form (two features on the array, BPP=0.9914; 0.9997) but not the gonad form of aromatase (BPP=0.3192). This result suggests that the elevated testosterone levels found in D males(Parikh et al., 2006) may affect aggression, courtship or dominance through aromatization and subsequent action via estrogen receptors in the brain. In birds, aromatase activity increases during the territorial period and correlates with aggression (e.g. Soma et al.,2003; Silverin et al.,2004). Blocking brain aromatase reduces male courtship in guppies(Hallgren et al., 2006),further suggesting estrogen-mediated neuroendocrine regulation of reproductive behavior for some species. However, studies on gonadal sex change in fish(Black et al., 2005; Marsh et al., 2006) suggest the opposite relationship between brain aromatase and male aggression and thus a more complex mechanism, possibly involving differences in receptor expression, binding proteins or anatomical localization. Estrogen receptors did not show differential regulation on the array, although they may have been expected to, according to Burmeister et al.(Burmeister et al., 2007). The inability to reliably detect differences in receptor gene expression is probably due to small localized effects that are masked by whole-brain gene expression levels.

We bioinformatically annotated the ESTs obtained from the cichlid cDNA library features represented on the microarray (see Materials and methods). Several of the genes thus identified fall into categories that represent candidates likely to play a role in the social regulation of a complex phenotype (Table 2). In the present study, these genes are considered to be `novel genes' rather than`candidate genes' because the annotation process does not involve rigorous manual curation of genes a priori, which was employed for the candidate genes discussed above. In addition to many genes involved in cellular metabolism that are differentially regulated between the two social phenotypes, we found genes encoding structural proteins, cell-cycle regulators, specific transcription factors, a plethora of neuropeptides,components of the neurosecretory machinery and neurotransmitter receptors.

Genes coding for structural proteins, such as tubulin and actin, and proteins that bind scaffold elements, such as septin 7 and ELF-1a, were more highly expressed in D males reflecting the observed differences in soma size between D and S for pre-optic neurons expressing GnRH1(Francis et al., 1993) and somatostatin (Hofmann and Fernald,2000). Furthermore, genes involved in axonal growth, neuromodulin[also known to play a role in modeling of sex-specific brain regions(Simerly, 2002)] and neuroserpin (Miranda and Lomas,2006), were also upregulated in D males. Taken together, the regulation of this gene set strongly suggests increased neuronal re-wiring in D males not previously reported and possibly similar in scale to the massive remodeling of neural circuits seen in seasonal accession to territoriality and mating accompanied by increased testosterone levels in song birds(Devoogd and Nottebohm, 1981;reviewed in Arnold, 1992). It is particularly intriguing that neuroserpins may play a role in anxiety and sexual behaviors. Specifically, neuroserpin-deficient rats showed decreases in exploratory behavior together with increases in anxiety and neophobia(Madani et al., 2003). In swordtail fish, neuroserpin expression increased in the brain of females exposed to an attractive male compared with females exposed to a non-attractive male (Cummings et al.,2008). This association of neuroserpin with social behavior is intriguing in that it may enable dominant males to approach and interact with novel stimuli such as competitors and potential mates.

Similarly, several cell-cycle regulators(Table 2) were significantly regulated in D and S phenotypes, suggesting that the extent of neurogenesis and subsequent cell death may also differ between these phenotypes, a hypothesis consistent with the finding that cell proliferation in the brain is correlated with high social status in rainbow trout(Sorensen et al., 2007). While there is currently no other evidence for plasticity of this kind in neuroanatomical structures outside the POA in A. burtoni, gross neuroanatomical differences that correspond to species' typical reproductive strategies have been identified in other cichlid species(Pollen et al., 2007).

Genes encoding neuropeptides and protein hormones that have not been previously studied in this system were perhaps the most striking, although not unexpected, class of genes regulated according to social status. In addition to the neuropeptides GnRH1, AVT and galanin discussed above, we found somatotropin, prolactin and somatolactin [all members of the growth hormone(GH) family of genes], as well as proopiomelanocortin (POMC), to be upregulated in D males. Interestingly, a similar pattern of endocrine gene regulation (GH, prolactin, somatolactin, POMC) is observed in Atlantic salmon,such that the expression profile for the early maturing `sneaker' male compared with immature males matches the profile observed in the present study for reproductive D males, suggesting conserved function of these pathways(Aubin-Horth et al., 2005). We found cholecystokinin (CCK) and natriuretic peptide to be upregulated in S males.

Since the activity of pituitary somatotrophes is associated with testis maturity and is stimulated by high levels of GnRH in several fish species(reviewed by Legac et al.,1993; Yu and Peter,1991), we suggest that the increased expression of growth-related genes in D males is probably related to gonad maturation. Somatolactin, which thus far has been found only in teleost fish, is involved in both growth(Forsyth and Wallis, 2002) and color change (Fukamachi and Sugimoto,2004), two plastic traits associated with social dominance in A. burtoni. Importantly, the observed increase in expression of GH is consistent with the previous finding that circulating GH levels are higher in D males (Hofmann et al.,1999b). Additionally, the growth hormone-releasing hormone(GHRH)/GH axis facilitates territorial behavior in A. burtoni(Hofmann et al., 1999b; Trainor and Hofmann, 2006)(note that GHRH was not represented on the array). Finally, antagonists of the neuropeptide somatostatin (which inhibits GH production and release in the pituitary; this gene was not represented on the array) inhibit aggressive behavior in A. burtoni males without affecting sexual behavior(Trainor and Hofmann,2006).

In addition to neuropeptide genes, we found many genes involved in production, maturation, release and reception of neuropeptides and neurotransmitters to be differentially expressed between social phenotypes. For example, secretory granule proteins, such as a homolog of synaptophysin as well as members of the granin family of acidic proteins, were upregulated in D males. These are notably found in a wide variety of endocrine and neuro-endocrine cells (for reviews, see Gerst, 1999; Helle, 2004), and the regulation pattern found in the present study may simply be a consequence of increased neuroendocrine activity in D males.

Although few neurotransmitter receptors were present on the array, there is one very intriguing result to report: a GABA-(A) receptor was upregulated in D males whereas at least two subunits of the kainate-type glutamate receptor were upregulated in S males. Clearly, these differences were not predicted nor do we have, at this point, a clear understanding of their specific roles in regulating social dominance. However, the fact that these receptors are affected by (or affect) social phenotypes probably reflects the interconnected nature of neurotransmitter systems and suggests that regulation at all levels of neural circuits may underlie the transition from one social phenotype to another. Interestingly, an antagonistic relationship between GABA and kainate signaling has been suggested in the regulation of reproductive physiology in other systems (e.g. Sagrillo et al.,1996; Chu and Moenter,2005; Clarkson and Herbison,2006; Liu et al.,2006) (see Molecular modules underlying dominance and reproduction).

The GO annotation scheme applied to the cichlid microarray allows rigorous statistical analysis for over- and under-representation of particular molecular functions, biological processes, and cellular locations in genes that are differentially expressed in each male phenotype. Despite the biased nature of the GO terms due to their origin and application in model organisms and directed research, this tool offers a mechanism for statistical analysis of microarray results according to function(Shaw et al., 1999). These terms, unlike specific gene names, avoid experimenter bias and cross-referencing between experiments and even between species and relate experimental results between organisms and platforms. Of the 171 features regulated according to male social phenotype, GO terms could be applied to 85. Analysis at all GO levels revealed 22 categories that are statistically over-or under-represented among the genes that are regulated by social phenotype(P<0.05; compared with their representation among all genes above threshold) (Fig. 3). Using permutation analysis, we determined that only five GO terms were expected to show significant over- or under-representation by chance alone (hypergeometric test P<0.05, only one GO term at P<0.01). Importantly,this unbiased statistical approach confirms our observation discussed earlier that cytoskeleton/structural molecules as well as hormone signaling are upregulated in D males (Table 3). Furthermore, cation/potassium transport pathways appear to be important building blocks for each male phenotype. In D males, the biological processes of GTP binding, iron ion binding and motor activity were significantly enriched, whereas in S males potassium ion transport, regulation of cellular cation transport and ligand-gated ion channel function were activated. Although difficult to interpret directly, this bioinformatic approach results in a considerable data reduction, facilitates comparisons across species and platforms, and provides a framework of hypotheses for future studies on the molecular underpinnings of socially regulated brain function.

The notion that biological entities (e.g. cognitive tasks, developmental programs, neural circuits, metabolic pathways) operate as functional and discrete (i.e. largely non-overlapping) units, or modules, is not new(Fodor, 1983; Redies and Puelles, 2001; Schlosser and Wagner, 2004; Op de Beeck et al., 2008). In molecular systems biology, a module can simply be defined as a set of co-regulated genes or proteins (Segal et al., 2004). Many such modules may serve as building blocks for the assembly of more complex processes. To date, most studies in this area have primarily been concerned with molecular and cellular networks and pathways in simple unicellular systems (Hartwell et al., 1999; Wolf and Arkin,2003). However, the ultimate challenge in the biology of complex systems is the integration across many levels of biological organization, from molecules to whole organisms, in an ever-changing environment. In the following, we make an initial attempt at such an integration of genomic data with physiological and behavioral phenotypes to provide a comprehensive conceptual framework for understanding phenotypically plastic traits. Specifically, we examine male dominance phenotypes in terms of molecular modules of socially controlled traits such as aggression, territoriality,reproduction and growth.

As we also obtained neural expression profiles of brooding females, we examined variation in transcript levels in relation to sex. 569 genes were regulated according to sex (316 male upregulated, 253 female upregulated)(Table 1), a number that is far greater than that observed for dominance phenotypes within males (171 genes or∼5%). This suggests that while the switch between dominance phenotypes is multifaceted (see Fig. 1), the difference in brain gene expression profiles between males and females is even more dramatic, affecting ∼16% of the genes on the array. Interestingly,this proportion is comparable to the 15% observed difference between alternative mating tactics in males of Atlantic salmon(Aubin-Horth et al., 2005). A considerable proportion of the 171 genes associated with male social phenotypes was also regulated according to sex (supplementary material Tables S1 and S2). However, among the social status-regulated genes, there were similar proportions of female-enriched and male-enriched genes(Fig. 4). Among the 87 dominant upregulated genes, were 11 female-enriched and 20 male-enriched genes(hypergeometric test; P=0.09), and among the 84 subordinate upregulated genes there were 16 female-enriched genes and 21 male-enriched genes (hypergeometric test; P=0.13). In other words, S males do not appear to be molecularly feminized nor are D males simply `super-males'. This result makes sense in light of the reproductive state of these animals. Although female behavior is, in many ways, similar to that of S males, both brooding females and D males are reproductively active. Furthermore, just as the metabolic demands of maintaining a territory are associated with reduced growth in D males (Hofmann et al.,1999a), mouth-brooding females starve while incubating their offspring and exhibit a marked reduction in body mass(Mrowk, 1984). By contrast, S males do not reproduce, and metabolic energy is directed toward growth(Hofmann et al., 1999a). Thus,the 11 genes (including synaptophysin, neuroserpin and GABA-receptor)upregulated in both D males and brooding females may be part of a module facilitating reproduction and/or reducing growth and may not necessarily be involved in dominance behavior per se.

Eighteen of the 316 male-enriched genes were even more highly expressed in D males compared with S males. We suggest that these genes, which include several structural proteins (actin, tubulin, ELF-1) and hormones (GnRH1, AVT,somatolactin, POMC) (supplementary material Table S1) are part of a`super-male' module, which probably plays a role in the aggressive and/or sexual behaviors typical of the D male(Fig. 4). Even though we know from previous studies using other approaches that the absolute levels of GnRH1(Francis et al., 1993; White et al., 2002) and GH(Hofmann et al., 1999b) are higher in D males, it is important to note that cDNA microarrays measure only relative transcript levels. Future studies will have to determine whether these genes must be activated in order to produce a dominant phenotype or whether their expression needs to be suppressed for the subordinate phenotype to occur. Following a similar rationale, we suggest that the 21 male-enriched genes (including neuroD and kainate receptor) that were further upregulated in S males relative to D males represent a masculinizing module in some sense. In this case, reduced expression would be critical for the manifestation of behaviors linked to aspects of social dominance and reproduction. Consequently, the 21 genes in this masculinizing module could be considered to be complementary to the super-male module. Finally, as brooding females and S males all display submissive behaviors and very little aggression, it is tempting to suggest that the 16 genes (including natriuretic peptide)expressed in a similar pattern in these two phenotypes are part of a`subordination module' (Fig. 4).

The opposing pattern of regulation for receptor expression in two classic neurotransmitter systems that we observed between male social phenotypes (see above) is also maintained in relation to sex. The gene encoding a GABA-(A)receptor was upregulated in D males and in females whereas the kainate-type glutamate receptor was upregulated in males in general and particularly in S males. There is a wealth of research that ties the GABA-(A) receptor to the regulation of the hypothalamic–pituitary–gonadal axis viaintegration of steroid feedback to GnRH neurons (for a review, see Sagrillo et al., 1996). In mammals, GABA has mixed inhibitory and excitatory effects on the release of GnRH, due in part to a developmental switch from GABA-(A) receptor depolarization to hyperpolarization(Clarkson and Herbison, 2006). In fish, depending on the species, GABA has either excitatory or inhibitory effects on GnRH release (Trudeau et al.,2000). Interestingly, glutamate-controlled GABA release has been implicated in GnRH regulation (Chu and Moenter, 2005; Clarkson and Herbison, 2006). The kainate system has also been proposed to underlie observed sex differences in the mechanisms of the neural glucocorticoid/stress response. Female mice show less atrophy of hippocampal neurons in response to elevated glucocorticoid levels, possibly due to the increased expression of NMDA, AMPA and kainate glutamate receptor sub-types(Liu et al., 2006). In A. burtoni, the increased kainate receptor expression seen in S males could similarly provide a neuroprotective effect against the elevated cortisol levels seen in S males during specific social situations(Fox et al., 1997). Future pharmacological and neurohistochemical experiments will elucidate the mechanistic interactions of these neurotransmitter systems in relation to sex and social behavior.

Taken together, the molecular systems analysis in the present study supports the notion that transcript patterns may indeed be organized in a modular fashion and can be strongly associated with behavioral and/or physiological traits associated with social phenotypes or sex in either concordant or contrasting ways. Additionally, we can exploit expression variation between phenotypes for tentatively annotating gene function and predicting functional roles of these genes.

To appreciate the importance of variation in gene expression for phenotypic plasticity, we need to evaluate the expression differences between individuals. To determine the extent to which individuals of the same phenotype differ in their expression profiles, we estimated transcript levels for each individual separately. Individual profiles were then clustered for similarity according to the estimated transcript levels using the gene list that had been identified as significantly regulated between any two phenotypes(BPP>0.99) (Euclidian distance matrix based on resampling for bootstrap confidence levels) (Fig. 5).

Sex was a strong factor in the clustering of individual expression profiles, as five of the six females clustered separately from all males. While males of the same social phenotype tended to have similar brain expression profiles, two D individuals (males 4 and 5) clustered with the S males. The analysis of behavioral data of these males did not suggest any obvious difference from other D males (Fig. 2). It is unlikely (since we controlled for age and size and based on our observations) that these two D males were experiencing dominance for the first time. Also, at the behavior level, these two D males were not similar to each other. While male 5 did show the fewest chases and the most border threats (Fig. 2) of all D males, the behaviors of male 4 were not in any way different from those of the other D males. While it is possible that there is more than one molecular substrate for constructing and maintaining a dominant male phenotype, the expression profiles of these two D males are no more similar to each other than to S males.

When individual profiles were clustered according to estimated transcript levels using all genes on the array that passed filtering for every fish,rather than only those genes that showed significant regulation, similar results were obtained (data not shown). While the bootstrapped confidence values were lower, the same four D males formed a cluster, suggesting that gene regulation according to sex and social status account for the greatest variation among individuals in this study. Principal component analysis (PCA)corroborated this conclusion (supplementary material Fig. S1).

While, overall, the clustering resulted in a robust separation of the three phenotypes according to sex and social status, the individual variation,already apparent at this level, prompted further inquiry into expression variation between individuals of the same phenotype. In voles, expression patterns of oxytocin receptor for females and vasopressin receptor for males is correlated with individual variation in social and anxiety-related behaviors (Olazabal and Young,2006). Similarly, male mice show individual variation in estrogen receptor distribution that correlates with aggressive behavior(Trainor et al., 2006). Previous studies in A. burtoni and other cichlid species have also reported strong covariation patterns within a social phenotype between the expression of candidate genes and specific phenotype measures [somatostatin correlated with aggression (Trainor and Hofmann, 2006); AVT correlated with hormone titers(Aubin-Horth et al., 2007)]. We therefore asked whether significant differences in gene expression existed between individuals of the same phenotype. Furthermore, we tested whether significant differences in gene expression between individuals of the same phenotype could be found among those genes that are differentially regulated between phenotypes.

In order to investigate the degree of individual variation in gene expression, we determined the number of genes that varied significantly in expression between individuals of the same phenotype and compared it with the number that varied between individuals of different phenotypes. Although statistical power was lower because we had only four technical replicates per individual as opposed to six biological replicates in the analyses above(Clark and Townsend, 2007), we were able to measure the average number of genes significantly regulated(BPP>0.99) in each possible pairwise comparison of two individuals within and between phenotypes. For intra-phenotype variations among males, a mean(±s.e.m.) of 82.8±4.5 (2.3% of all genes analyzed) genes varied between any two D males, while a mean of 92.3±4.3 (2.6%) genes varied between any two S males. The identity of these genes was substantially different for each pairwise comparison, such that 38% of all the array features varied in at least one intra-phenotype comparison. Interestingly, for inter-phenotype variation between individual males, a mean of 132±8.4(3.6%) genes varied between any two males of differing social phenotype, which was not significantly different from the intra-phenotype variation(t-test, P=0.13). Although it is difficult to set an equivalent threshold for significant variation between individuals and between phenotypes, this high degree of individual variation is consistent with other studies that have examined genome-scale individual differences in gene expression in order to study the molecular basis of natural variation. Whitehead and Crawford found that 69% of the metabolic pathway genes showed significant variation between individuals within a population, while only 12%were significantly regulated between populations adapted to different temperatures (Whitehead and Crawford,2006). Similarly, in yeast, up to 50% of the expressed genes show significantly different levels of expression among individual strains(Brem et al., 2002). Not only do we find a similar number of genes to be regulated between individuals of the same or different phenotypes, we also find no significant difference in coefficient of variation of expression level for sets of these individually regulated genes. This result indicates that absolute gene expression levels vary between individuals of the same phenotype as much as between phenotypes.

While low variation within a phenotype for those genes that define that phenotype can be expected, an alternative hypothesis posits that those genes that define the social phenotype vary across individuals displaying that phenotype in a manner associated with variation in physiological and behavioral traits (e.g. Trainor and Hofmann, 2006; Aubin-Horth et al., 2007; Cummings et al.,2008). In support of this notion, we found a statistically significant over-representation of intra-phenotype regulated genes among those that were regulated by social status (Table 4). Approximately 72% of the 87 genes upregulated in the dominant phenotype were also significantly regulated among individuals within a phenotype (41 among D and 39 among S, 24 of which are shared). Similarly, 64%of the 84 genes that were upregulated in the S phenotype were also significantly regulated among individuals (47 among D and 52 among S, 38 of which are shared). This variation cannot be explained by technical variation in array hybridizations: for a given animal, an array feature must show consistent results across four dye-reversed hybridizations before it can be identified as regulated across individuals according to our statistical analysis. Rather, the results of the present study show that even considerable and potentially important within-phenotype variation in gene expression can give rise to reliable and readily identifiable between-phenotype differences. Future integrative studies will help determine whether the observed variation between individuals is caused by, or causes, subtle phenotypic differences or represents a dramatic, alternative molecular mechanism for constructing the same phenotype.

By including candidate genes on the microarray, we have validated the discovery-based approach that identifies the gene expression patterns for reproductively active D males and reproductively suppressed S males as well as for brooding females. While the regulation of neuroendocrine genes was predicted from previous research, the unexpected and novel discovery of opposing roles for two classic neurotransmitter systems opens exciting new avenues for future research. Furthermore, the increased neuronal remodeling activity in D males suggests that the known neuroanatomical changes observed in pre-optic GnRH and somatostatin neurons may extend to other cell types and/or additional brain regions. Using the GO framework to interrogate the data for statistical significance reinforces the importance of hormonal regulation and highlights the hitherto under-appreciated roles of cytoskeletal components in addition to neurotransmitter pathways. Interestingly, expression profiles vary among individuals within a male phenotype for roughly the same number of genes and with similar magnitudes of transcript abundance as seen between male phenotypes. Specifically, there is a surprisingly high level of variation among those genes that molecularly define the very phenotypes in the first place. Taken together, our genome-scale analysis of molecular systems in the brain has identified complex patterns of gene expression that are associated with a socially regulated switch in behavioral phenotype.

LIST OF ABBREVIATIONS

 
• AMPA

  alpha-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid

 
• AVP

  arginine vasopressin

 
• AVT

  arginine vasotocin

 
• BPP

  Bayesian posterior probability

 
• CF

  condition factor

 
• CCK

  cholecystokinin

 
• DAG

  directed acyclic graphs

 
• DI

  dominance index

 
• EST

  expressed sequence tags

 
• FDR

  false discovery rate

 
• GABA

  gamma-aminobutyric acid

 
• GH

  growth hormone

 
• GHRH

  growth hormone-releasing hormone

 
• GnRH

  gonadotrophin-releasing hormone

 
• GO

  gene ontology

 
• GR

  growth rate

 
• GSI

  gonadosomatic index

 
• IEA

  inferred electronic annotation

 
• ISS

  inferred sequence similiarity

 
• KT

  ketotestosterone

 
• NMDA

  N-methyl-D-aspartic acid

 
• POA

  pre-optic area

 
• POMC

  proopiomelanocortin

 
• RI

  reproductive index

 
• T

  testosterone

 
• TC

  tentative contigs

All authors contributed to the design, implementation, analysis and presentation of this experimental data set. We are grateful to Josiah Altschuler and Melinda Snitow for animal care, Sarah Annis for assistance with the behavioral experiments, Christian Landry for programming and photography,Amir Karger for bioinformatics advice, and Victoria Zero, Brian Dias, Lin Huffman and Kim Hoke for comments on earlier versions of this manuscript. We thank the H ofmann laboratory and the members of the Bauer Center for Genomics Research for stimulating discussions. This research was supported by an NRSA post-doctoral fellowship to S.C.P.R., a postdoctoral fellowship from FQRNT(Fonds Québécois de la Recherche sur la Nature et les Technologies) and a postdoctoral fellowship from the Natural Science and Engineering Research Council of Canada (NSERC) to N.A.H., and by National Institutes of Health grant NIGMS GM068763, the Bauer Center for Genomics Research (H.A.H.).