Comparative gene expression profiles for highly similar aggressive phenotypes in male and female cichlid fishes (Julidochromis) Free

Parental care demands an investment of parental resources, and represents a trade-off in fitness for each parent (Trivers, 1972). In biparental species, males and females share the fitness burden of parental care but, despite this, males and females frequently differ in the amount and nature of parental care they provide, showing strong biases such that sex-specific behaviors develop. For a wide range of biparental species that exhibit sex-specific behaviors, it is common to find males engaged in territorial defense while females are involved in more direct care of offspring [e.g. birds (Burger, 1981; Creelman and Storey, 1991; Fraser et al., 2002); and fish (Keenleyside and Bietz, 1981; Lavery and Reebs, 1994; Wisenden, et al., 1995; Itzkowitz et al., 2001; Itzkowitz et al., 2003; Itzkowitz et al., 2005)].

In a minority of species, the conventional bias of sex-specific behavior is reversed such that females engage in aggressive behavior while males perform the majority of direct parental care [e.g. pipefish (Berglund and Rosenqvist, 2003); and Wattled Jacanas (Emlen and Wrege, 2004)]. Reversal of sex-specific behavior in these species is often assumed to be the result of females having a higher potential reproductive rate (Eens and Pixten, 2000; Cluttonbrock and Vincent, 1991), ultimately resulting in stronger sexual selection on females. This ultimate explanation has rarely been evaluated at a genetic level (Jones et al., 2001) but is well supported by studies that show plasticity of sex-specific behavior under ecological fluctuations that reverse the direction of sexual selection (Forsgren et al., 2004; Gwynne, 1985). Even when these ecologically based proximate causes of reversed sex-specific behavior are identified, the mechanistic explanations often remain elusive. Despite investigation of hormone levels [e.g. pipefish (Mayer et al., 1993); and birds (reviewed in Eens and Pixten, 2000)] and gene expression (Voigt and Goymann, 2008), it has been difficult to determine the extent to which the mechanisms that underlie aggressive behavior in females are similar to those that underlie aggressive behavior in males. Species that show reversed sex-specific parental care behavior provide the opportunity to examine the proximate mechanisms that underlie reversed phenotypes.

The African cichlids provide an excellent model with which to address the proximate mechanistic explanations of sex-specific behavior. Their adaptive radiation has been accompanied by a vast amount of behavioral diversification, including several independent transitions in parental care strategy (Goodwin et al., 1998). A number of studies demonstrate that within the family Cichlidae biparental behaviors are divided into stable sex-specific behavioral phenotypes. The existing body of research suggests that biparental cichlids generally show sex-biased behavior, with females performing the majority of parental care and males utilizing their increased size for territory and brood defense (Keenleyside and Bietz, 1981; Itzkowitz, 1984; Kuwamura et al., 1986; Kuwamura et al., 1997). The cichlid genus Julidochromis is composed of biparental substrate brooders (Konings, 1998) that show variation in specialized parental care such that several species in this genus exhibit the conventional pattern of sex-specific behavior while Julidochromis marlieri shows reversed sex-specific behavior (Barlow and Lee, 2005). Julidochromis marlieri females are large and highly aggressive (Yamagashi and Kohda, 1996; Barlow and Lee, 2005) (K. Wood and S.C.P.R., unpublished), and J. marlieri males spend significantly more time in the nest when eggs are present (Schumer, 2009). Although estimates vary, according to recent phylogenetic analysis, Julidochromis transcriptus and J. marlieri diverged ∼3.7 million years ago (Sturmbauer et al., 2010). In our present study, the molecular mechanisms that underlie sex-specific behavior in J. marlieri males and females were directly compared with those of another Julidochromis species, J. transcriptus, which has conventional sex-specific behavior (Awata et al., 2006; Schumer, 2009). The present study focuses on investigating whether the proximate mechanisms underlying dominant and subordinate behavior are the same in J. transcriptus and J. marlieri or whether these mechanisms differ due to reversal of sex-specific behavior in J. marlieri.

Recent advances in the ability to identify the genes involved in complex phenotypes means that we are no longer constrained to mutant screens in model organisms [circadian rhythm in Drosophila (Konopka and Benzer, 1971); chemotaxis in Caenorhabditis elegans (Dusenbery, 1980); and mating in Drosophila (Yamamoto et al., 1997)] and can even move beyond the identification of natural allelic variation [foraging behavior in Drosophila (Osborne et al., 1997)] to take on a genome-wide comparative approach. The development of reliable methods to globally assay gene expression levels and associate these gene expression patterns with behavior has revolutionized our understanding of the relationship between gene expression and complex behaviors such as aggression (Edwards et al., 2006), addiction (Lewohl et al., 2000; Rhodes and Crabbe, 2005; Buitenhuis et al., 2009), maternal care (Gammie et al., 2007), social behavior (Whitfield et al., 2003) (reviewed by Robinson et al., 2008), and even mate choice (Cummings et al., 2008). Among African cichlids, microarray techniques have been used to examine the basis of a variety of behaviors including dominance and subordinate behavior in male Astatotilapia burtoni (Renn et al., 2008), various phenotypes in Neolamprologus pulcher breeding groups (Aubin-Horth et al., 2007) and different mating strategies in species from the Ectodini clade (Machado et al., 2009). Combined, those studies have identified a number of genes whose expression strongly correlates with aggressive behavior or with sex, thus demonstrating that certain sex-specific behaviors in cichlids are accompanied by differences in gene expression.

Few studies to date have compared gene expression patterns between similar phenotypes in different species. Microarrays can be used to ask these comparative questions but become less useful for comparing phenotypes between species as genetic divergence increases (Machado et al., 2009; Renn et al., 2004; Cohen et al., 2007). In closely related species, microarrays have been successfully used to investigate whether differences in phenotype are reflected by changes in phenotype-specific gene expression. Honeybee species have slightly different species-specific foraging dances; gene expression profiles of three species producing the foraging dance revealed striking similarities between species but also identified a number of species-specific differences in gene expression in genes related to motor control and synaptic activity (Sen Sarma et al., 2009). Comparative work investigating the genes underlying a similar foraging behavior in honeybees and paper wasps showed that a small set of genes was shared between these independently evolved phenotypes (Toth et al., 2010). Such studies expand the potential questions posed in microarray studies from those focused on assessing gene expression differences within species to broader questions about the evolution of phenotypic differences between species.

The main aim of the present study was to identify the proximate mechanisms for the variation in sex-specific roles observed among species of the genus Julidochromis. We asked whether when a female J. marlieri achieves the conventionally male-like dominant phenotype, she does so through a pattern of neural gene expression similar to that of a J. transcriptus male or whether this reversal of the behavioral phenotype is orchestrated through a novel pattern of brain gene expression. The former would suggest the co-option of homologous molecular mechanisms to produce a male-like dominant phenotype in these females whereas the latter would suggest convergent evolution to produce superficially similar behaviors in the opposite sexes of these two lineages.

Using cDNA microarrays we characterized differences in brain gene expression between sex within each species to identify genes that are associated with sex-specific behavior in each species. We then compared the list of female-biased genes in J. marlieri with the list of male-biased genes in J. transcriptus to ask whether female aggression in J. marlieri relies upon the same genetic mechanisms as male aggression in J. transcriptus. We also compared the list of male-biased genes from J. marlieri with the list of female-biased genes from J. transcriptus to identify genes that are associated with subordinate behaviors, regardless of sex and species. Keeping in mind that downregulation of gene expression level may be equally important for the observed phenotypes, the observed regulation of these subordinate-biased genes may also be important for the dominant phenotype. Although the majority of gene expression was both phenotype and species specific (∼10% of differentially expressed genes shared between phenotypes), we found convincing evidence for a core set of genes regulating aggressive behavior with conserved expression between species.

Fish were maintained at pH 8.0 and 650 μS salinity with a 5% daily water change, and light cycle including 11.5 h of light with a 30 min fade for dawn and dusk to mimic conditions in nature. Pairs of Julidochromis marlieri, Poll 1956, and Julidochromis transcriptus, Matthes 1959, were allowed to form naturally from mixed male–female populations in 113.6 l tanks. Two individuals were considered to have formed a ‘pair’ when they were observed over multiple days to have established a territory together. After a pair was established it was transferred to a 113.6 l observation tank that consisted of a neighboring pair of the same species separated by a transparent acrylic divider. Each pair was provided with ∼2 cm of gravel substrate for digging and an artificial nest crevice. The nests were made of two clay tiles (15×15×1 cm) with an entrance width of 8 cm, designed to mimic the caves and crevices that are used as nests in the wild (Awata et al., 2006). Julidochromis transcriptus and J. marlieri stocks had been maintained in-house for three years at the time of the study. First filial generation J. transcriptus individuals were obtained from Atlantis Tropical Fish Hatchery (Gardiner, NY, USA), and J. marlieri individuals were obtained from George Barlow, Berkeley (CA, USA), formerly at UC Berkeley (age uncertain).

All experiments were carried out according to animal protocol IACUC no. 1032007.

Eight pairs (four pairs of each species) were checked every morning between 08:00 h and 10:00 h for eggs and observed weekly using the event recorder Jwatcher (Blumstein and Daniel, 2007). Ten-minute behavioral observations were used to confirm that individuals were exhibiting the sex-specific behavior typical of their species. Males and females within a pair were observed simultaneously. The ethogram included event measures of aggressive behaviors (attack intruder, attack mate, approach intruder) and parental care behaviors (egg cleaning) as well as duration measures of parental care behavior (time spent in the nest). The total number of observations for each pair varied based on when eggs were laid but ranged from 3 to 11 observations. Due to the differences in the number of observations, a generalized linear model was used with the lme4 package in R (http://cran.r-project.org/web/packages/lme4/index.html) to test for differences in male and female behavior in each species and for overall differences in behavior between species. Two post-spawning behavioral observations were conducted for each pair, one on the day of egg laying and one 24 h later.

Twenty-four hours after a pair had laid eggs, the male and female were simultaneously removed from the observation tank, anesthetized using MS-222 (160 mg of MS-222 in 500 ml cichlid salt water, pH 8.0) and killed according to protocol (IACUC #2007.103). This time point was chosen in order to ensure synchrony of reproductive cycles and capture parental care. After the fish had been anesthetized, two researchers simultaneously dissected the two fish. Following decapitation, the brain was immediately removed and stored in 1 ml RNAlater (Ambion, Austin, TX, USA). Brains were stored at 4°C for 24 h and then at –20°C. Gonads were removed and stored in 0.5 ml paraformaldehyde at 4°C. Samples were collected from four pairs of each species.

Brains were homogenized (Tissue tearor, Biospec Products, Bartlesville, OK, USA), and total brain RNA was isolated using a standard Trizol (Invitrogen, Carlsbad, CA, USA) protocol adapted for use with phase-lock gel tubes (Eppendorf, Westbury, NY, USA). RNA was quantified with a Nanodrop 1000 (Thermo Scientific, Wilmington, DE, USA). Each total RNA sample, ∼1000 ng, was linearly amplified with a single round of Amino Allyl MessageAmp II aRNA amplification (Ambion), according to manufacturer's protocol (yield 31- to 74-fold). The antisense amplified RNA samples, 4 μg each, were then dye-coupled with Cy3 and with Cy5 and quantified with the Nanodrop 3300 fluorospectrometer (Thermo Scientific) prior to storage at –80°C. One limitation of using whole brain RNA is the possibility that this method will fail to detect certain changes in gene expression, due to discrete anatomical localization or counter regulation in differing brain regions; therefore, subtle regulation of gene expression may be missed by our technique (Filby et al., 2010).

In each hybridization, equal quantities of Cy3 and Cy5 samples (1.5–2.5 μg) were combined in hybridization buffer (Ambion SlideHyb buffer #1) and applied to the second-generation cDNA array (GEO platform ID: GLP6416) designed for the African cichlid A. burtoni (Renn et al., 2004; Salzburger et al., 2008), which contains ∼20,000 features, representing ∼16,000 unique sequences (including ∼35% anonymous features) derived from brain, skin and mixed organ libraries from individuals of all ages, sexes and phenotypes, and has been validated for a variety of cichlid species (Renn et al., 2004; Aubin-Horth et al., 2007; Machado et al., 2009). Cross-species comparative genomic DNA hybridization for these Julidochromis species identified fewer than 10 microarray features that are substantially influenced by sequence or copy number differences between these species (data not shown). Given the within-species design, even these features should not affect the interpretation of results. For each species, the eight individuals were compared between sex and within species (Fig. 1) in a modified loop design (Churchill, 2002) incorporating dye swaps for greater statistical power. Each sample was used for 2–3 hybridizations, one of which was to their mate.

Microarrays were incubated at 48°C for 16–18 h and then rinsed through three successive wash buffers and spun dry prior to scanning on a GenePix 4000B scanner (Molecular Devices, Sunnyvale, CA, USA) with GenePix 5.1 imaging program. Poor spot morphology and technical artifacts were flagged manually. These raw data were imported into R software (v.2.12) (R Development Core Team, 2006) for analysis with the Linear Models for Microarray Data package (LIMMA) (Smyth, 2004). Array features were filtered for poor spot morphology, hybridization artifacts and low intensity (<2 s.d. above local background). Background subtracted intensities were normalized within array using the print tip LOESS method, and replicate features representing the same DNA sequence were averaged and weighted. Intensities were grouped for analysis using the lmFit function adapted to calculate exact P-values in the presence of missing data (A. Jones and S.C.P.R., unpublished), and evaluated for sex-specific bias using the empirical Bayes method (eBayes) (Smyth, 2004). Only genes surviving quality filters in sufficient replicates for analysis in both species were considered for the analysis of differential regulation. Genes deemed significant at P<0.05 were investigated using a gene index annotation developed for the A. burtoni cDNA array (Renn et al., 2004) (http://compbio.dfci.harvard.edu/tgi/cgi-bin/tgi/gimain.pl?gudb=a_burtoni); no multiple hypothesis testing correction was used due to small sample size. Lists of significantly regulated genes for each phenotype were compared using one-line perl scripts provided by FAS Center for Systems Biology (http://sysbio.harvard.edu/csb/resources/computational/scriptome/). All between-species comparisons were done through gene lists rather than statistical comparison, due to the fact that no between-species hybridizations were performed in order to prevent hybridization biases caused by species-specific sequence differences. To test for significant overlap between different gene lists, a Fisher's exact test was used. All raw and processed data are available at the GEOdatabase [(www.ncbi.nlm.nih.gov/projects/geo/) Series: GSE23094 arrays, GSM569290–GSM569309].

Gene-specific primers were designed for four genes (Table 1) using sequence information from A. burtoni: parvalbumin (DY626384), isotocin precursor (DY631081), microtubule-associated protein-1B (MAP-1B, DY628922) and glyceraldehyde 3-phosphate dehydrogenase (GAPDH). Primers were designed using Primer3Plus and optimized using IDT Oligo Analyzer (http://www.idtdna.com). Primer specificity (single band of predicted size) was verified for both Julidochromis species using genomic DNA (gDNA) as a template. RNA samples were available in sufficient quantity for technical validation by RT-qPCR using all eight J. transcriptus individuals and seven of the eight J. marlieri individuals used in the original hybridizations. DNA template was synthesized using oligo (dT) for total RNA and random hexamers for amplified antisense RNA (iScript, Biorad, Hercules, CA, USA). Per reaction volumes were: 2 μl 5× iScript reaction mix, 1 μl primer, 0.5 μl iScript reverse transcriptase, 25 ng μl^–1 of RNA from each individual, and water (added to bring the final reaction volume to 10 μl). For an oligo (dT) primed reaction, the mix was incubated at 42°C for 90 min and then 85°C for 5 min. For a random hexamer primed reaction, the mix was incubated at 25°C for 5 min, followed by an incubation at 42°C for 30 min and finally at 85°C for 5 min. Reaction conditions for RT-qPCR were optimized by creating a standard curve using different primer and template concentrations; efficiency for all primers ranged from 90% to 100%. RT-qPCR (Opticon; Biorad) was performed on cDNA using GAPDH for normalization. Per reaction volumes were: 5 μl 2× ImmoMix (Bioline, Taunton, MA, USA), 0.15 μl 50× SYBR Green, 5 μmol l^–1 of forward primer and reverse primer (Integrated DNA Technologies, San Diego, CA, USA), 2.85 μl nanopure water, and 1 μl template cDNA. All qPCR reactions were performed as follows: 1 cycle of 95°C for 10 min, 35 cycles of 94°C for 10 s, 57–60°C for 20 s, 72°C for 15 s, followed by a hold at 72°C for 2 min. At the end of the reaction a melt curve analysis was performed for each sample. Each sample was run in triplicate and each plate contained no RT controls. A linear model was generated using the lm function in R to look for significant correlation in gene expression between qPCR and microarray results.

Behavioral observations confirmed that J. transcriptus pairs show conventional sex-biased behavior whereas in J. marlieri this behavior is reversed. This allowed us to assess our central question: is aggressive behavior regulated by similar or divergent gene expression mechanisms in both Julidochromis species despite reversal of sex-specific behavior? Overall, microarray results support the conclusion that gene expression in J. marlieri has diverged from J. transcriptus. However, we also found a core set of genes associated with the aggressive phenotype in both species, some of which are confirmed by qPCR. We found evidence of sex-biased gene expression and phenotype-specific gene expression that included a number of candidate genes previously implicated in aggression or parental care in a variety of species (Table 2). The territorial/aggressive phenotypes of both species showed the greatest number of differentially expressed genes. These dominance-biased genes tended to be those of greater fold change, and overlap between the gene lists for the two aggressive phenotypes of each species was highly significant. Functional inference can be drawn from the available annotations (Table 2) (http://compbio.dfci.harvard.edu/tgi/cgi-bin/tgi/gimain.pl?gudb=a_burtoni), and regulated genes with no known homolog represent novel candidates for the molecular control of behavior.

The pairs that bonded conformed to expected sexual size dimorphism for both species. Julidochromis transcriptus males weighed more (2× mean; P=0.006) and were larger (1.3× mean; P=0.008) than their mates whereas J. marlieri females weighed more (2.6× mean; P<0.001) and were larger (1.5× mean; P<0.001) than their mates. There was no significant difference in gonadosomatic index (GSI) between species.

In order to confirm species-specific patterns of sex-biased behavior, approach conspecific, attack conspecific and attack mate behaviors were totaled to quantify aggression. Due to the different numbers of observations per pair, a generalized linear model was used for each species to assess whether variation in behavior was best explained by sex, mate or individual. Based on this analysis, J. transcriptus males were found to be significantly more aggressive than females (P<0.0001, Fig. 2). Conversely, J. marlieri females were significantly more aggressive than males (P=0.00083, Fig. 2).

We examined the post-spawning time point based on the assumption that sex-specific behavior would be more pronounced after eggs were laid, because maximizing effective parental care is one of the hypothesized functions of sex-specific behavior. As predicted, more parental care was provided by J. transcriptus females (P=0.057) and by J. marlieri males (P=0.042) as measured by the time in the nest after spawning. This suggests that the subordinate individual of each species provides the majority of parental care. The measure of direct parental care, egg cleaning, which occurred infrequently, did not show a statistically significant sex bias as had been predicted. While J. marlieri post-spawning females tended to increase aggression (Fig. 2), as predicted, post-spawning J. transcriptus males tended to reduce aggression. Therefore, the sampling time point for gene expression may not have measured the peak in aggressive behavior in J. transcriptus. However, as post-spawning J. transcriptus males stayed closer to the nest, there was reduced opportunity to capture aggressive behavior with our ethogram. Although untested, we predict that behavioral assays including intruder introductions and territory challenges might reveal an increase in post-spawning aggression for J. transcriptus males.

Based on a generalized linear model including both species, J. marlieri females were significantly more aggressive than J. transcriptus males in both control (P=0.0484) and post-spawning conditions (P=0.0046) while J. transcriptus individuals spent more time in the nest with eggs compared with J. marlieri individuals (P=0.0669). Overall, the behavioral results confirmed that J. marlieri females and J. transcriptus males are more aggressive while J. marlieri males and J. transcriptus females spend more time in direct contact with eggs.

Males and females of the same species share the majority of their genome but can differ substantially in their sex-specific behavior, suggesting an important role of gene expression. In some studies, sex-specific gene expression represents the majority of within-species variation in gene expression (Yang et al., 2006; Jin et al., 2001) while other studies find few sex-specific differences in gene expression, even when adults are sexually and behaviorally dimorphic (Santos et al., 2008). Similar to the latter example, the present study found that despite the sexual dimorphism in behavior observed in both Julidochromis species, males and females differed in expression for a minority of the genes analyzed by microarray (∼7.5% in J. transcriptus, ∼9.3% in J. marlieri). Although few in number, in part owing to small sample size, the genes that were differentially regulated between males and females in each species present good candidates for the molecular basis of behavioral and physiological differences between males and females in J. transcriptus and J. marlieri.

We found a number of genes that were differentially regulated by sex in each species (P<0.05). Out of 9528 genes analyzed in both species (5111 with sequence information), 301 showed female-biased expression and 412 genes showed male-biased expression in J. transcriptus. In J. marlieri, 529 genes showed female-biased expression and 356 genes showed male-biased expression. Approximately 15% of the genes analyzed were differentially regulated by one of the four phenotypes. The species difference in the number of regulated genes is not likely due to statistical power as witnessed by the similarity in gene expression level 50 (GEL50) values (Clark and Townsend, 2007) for these two experiments (J. transcriptus 0.2397; J. marlieri 0.2129). Gene lists from males and females of both species were compared to identify any genes that were sex specific in both species. Genes within these sex-specific modules may be related to reproductive physiology. There was a trend towards more genes in the female module (24 genes) than expected by chance but not in the male module (21 genes) (Fisher's exact test: female, P=0.07198; male, P=0.1426, Fig. 3). The 45 sex-biased genes consistent between the two species are much fewer than the genes that were found to be phenotype specific (see below). Similarly, a recent study of monogamous and polygynous cichlids found that most differences in gene expression were both sex and species specific, with few genes showing sex-specific gene expression across species (Machado et al., 2009). In a previous review on the topic, Ellegren and Parsch (Ellegren and Parsch, 2007) reported that many of the genes that show sex-specific expression within species are those that are most sequence diverged between species, suggesting functional and therefore likely expression divergence. Therefore, based on previous research, it is not surprising that so few genes were sex specific across species.

A number of the sex-biased genes have a known role in reproduction and mating behavior. Gonadotropin-releasing hormone (GnRH), which is involved in sexual differentiation and mating behavior in organisms from fish to mammals, was strongly male biased in both Julidochromis species. While GnRH has been implicated in dominance because it is upregulated in A. burtoni territorial males compared with non-territorial males (Renn et al., 2008), elevated GnRH in territorial male A. burtoni may be related to gonad growth (see Soma et al., 1996) rather than behavior per se. In order to fully understand the role of GnRH in Julidochromis, it would be necessary to assess sex-specific expression throughout the reproductive cycle. In many species, including A. burtoni females, GnRH varies across reproductive state (White and Fernald, 1993); nonetheless at the time point assayed here, it is interesting to note that GnRH is expressed at higher levels in males. Prolactin, another male-biased gene in both Julidochromis species, has been linked to parental care behavior in a range of species [e.g. birds (Garcia et al., 1996); and rats (Bridges et al., 1990)]. This male bias is at first surprising given that females show greater parental care than males in J. transcriptus. However, despite the fact that subordinate individuals in Julidochromis spent the most time with the eggs, it is known that males and females of both species provide some parental care. Therefore, the male bias in expression may be consistent with the observation that in fish prolactin is particularly important for paternal rather than maternal care (Schradin and Anzenberger, 1999). The bias in prolactin expression in our present study of monogamous species is inconsistent with findings from a related cooperative breeding cichlid species (Bender et al., 2008).

The present study is among the first to consider whether similar phenotypes in related species are the product of similar gene expression profiles or whether these phenotypes are regulated by different mechanisms. Specifically, we asked whether aggressive behavior is regulated by similar or divergent gene expression mechanisms in both Julidochromis species despite reversal of sex-specific behavior. Other studies addressing similar questions have reached varying conclusions. While two limnetic salmonid species show similar activation of genes important for survival in that niche (Derome and Bernatchez, 2006), few genes show consistent activation according to mating system among another cichlid clade (S.C.P.R. and H. A. Hofmann, unpublished).

In order to address the question of phenotype-biased gene expression, we compared gene lists instead of using direct hybridization to reduce problems associated with interspecific hybridization on the microarray. Dominance-related genes were identified as those that were highly expressed in the aggressive male phenotype in the conventionally sex-biased J. transcriptus and that were also highly expressed in the aggressive female phenotype in the reversed sex-biased species J. marlieri (Fig. 3). Aggressive J. transcriptus males and J. marlieri females shared increased expression of 76 dominance-related genes (34 with sequence information), significantly more than expected by chance (Fisher's exact test: P<2.2e^–16). Conversely, 24 subordinate-related genes were identified as the intersection of the gene lists for highly expressed genes in J. marlieri males and in J. transcriptus females (16 with sequence information), also significantly more than expected by chance (P=0.00049). The fact that only a proportion of the genes that are associated with either social phenotype in each species are shared between species (Fig. 3) suggests that, since J. marlieri and J. transcriptus diverged from their common ancestor, other changes have also evolved in the gene expression patterns regulating these behaviors. Although few in number, it is notable that both of the phenotype-specific modules had highly significant overlap whereas the sex-specific modules did not. This leads to the exciting interpretation that the few genes that do show phenotype-specific gene expression across species are in fact major regulators of behavior. If the genes that are associated with dominate and subordinate phenotypes in both species are major regulators of behavior, this would support the hypothesis that the sex-specific phenotypes in J. transcriptus and J. marlieri are controlled by conserved mechanisms. The dominance module (i.e. the overlapping list of genes that are both upregulated in J. marlieri females and in J. transcriptus males) contained more genes than the opposing subordinate related module or sex-related module, and the overlap between gene lists for the dominance module was more significant than that for any other module (P<2.2e^–16). Furthermore, an examination of gene expression patterns (Fig. 4) reveals that many of the most significantly regulated genes in J. transcriptus males and J. marlieri females are part of the dominance module. These patterns of shared regulation in J. transcriptus and J. marlieri suggest that the functional mechanisms regulating dominance are conserved between species.

Following from the overarching hypothesis regarding gene expression patterns, the individual genes within the dominance- and subordinate-related modules are of great interest as potential key regulators of behavior. Two such genes, highly expressed in both male J. transcriptus and female J. marlieri, are already known to be involved in brain function. Isotocin, a paralog of arginine vasotocin (AVT) and a homolog to mammalian oxytocin, is a candidate for the regulation of dominant and aggressive behavior in Julidochromis. Although the few studies that have investigated a role for isotocin in aggressive behavior have found no correlation (Santangelo and Bass, 2006), its expression pattern here suggests otherwise. Parvalbumin, another identified gene in the dominance module, is expressed in GABAergic neurons (de Almeida et al., 2005), which have been linked to aggressive behavior in a number of mammalian species. The mechanism by which parvalbumin expression might alter γ-aminobutyric acid (GABA) activity to influence aggression is not clear.

AVT was a gene predicted to be upregulated in both J. transcriptus males and J. marlieri females based on extensive evidence supporting a role for AVT in aggressive behavior. However, AVT was significantly upregulated only in J. marlieri females (P=0.0164) but not in J. transcriptus males (P=0.349). Although there are dominant phenotypes common to the two Julidochromis species, there are also species-specific differences in aggressive behavior; in both pre- and post-spawning conditions J. marlieri females were significantly more aggressive than J. transcriptus males. AVT is upregulated in dominant A. burtoni males (Greenwood et al., 2008; Renn et al., 2008) and dominant individuals of both sexes in N. pulcher (Aubin-Horth et al., 2007), a relative of Julidochromis. Furthermore, manipulation of AVT alters territorial behavior in the tropical damselfish (Santangelo and Bass, 2006). These studies establish a clear relationship between AVT and aggressive behavior, and support the hypothesis that AVT plays an important role in regulating aggressive behavior in J. marlieri, the most aggressive phenotype in our study.

We selected genes for verification of array results using qPCR on both RNA and amplified RNA from the Julidochromis individuals. There was good correlation between results from RT-qPCR and microarray data for parvalbumin (F=4.611, P=0.05732) and isotocin (F=51.3, P=3.06e^–05). However, RT-qPCR results for MAP-1B differed from microarray results substantially. The correlation found between the amplified RNA and total RNA RT-qPCR results for isotocin (F=6.728, P=0.0268) and parvalbumin (F=4.511, P=0.05962) shows that the amplification process does not significantly change the detection of relative expression. Given the lack of correlation between amplified RNA and total RNA for MAP-1B, in addition to the lack of correlation with microarray results, we suspect that either the primers or array features or both are showing non-specific effects.

Perhaps the most exciting future direction for comparative microarray studies is the potential to investigate the identified gene modules in light of other studies and to use these results to infer biological function and identify novel candidates. Nearly half of the genes found in the sex-specific or phenotype-specific modules have no sequence information, while others have no known homolog in GenBank (Table 3). Array studies such as this one provide a tentative functional annotation for many of the genes of unknown function that are on the array (Landry and Aubin-Horth, 2007). The current study identified a number of new genes that are likely to be important in aggressive behavior and parental care in Julidochromis and other cichlids, which previously had not been identified (see Table S1 in supplementary material). To see whether any of the putative dominance-related genes in this study had been previously identified as aggression-related genes, we compared our data with a previous study conducted with a smaller microarray comparing A. burtoni territorial males, non-territorial males and females (Renn et al., 2008). We predicted that because the polygynous A. burtoni has conventional sex-specific behavior, its pattern of sex-specific gene expression may align more closely with that of the conventional J. transcriptus rather than the reversed J. marlieri. Unexpectedly, we found that J. marlieri females showed a trend towards sharing significant gene regulation with A. burtoni males (P=0.0966), while no other phenotype had significantly greater shared regulation with dominant A. burtoni males than expected by random chance. In addition to AVT, 14 genes with no functional information were upregulated in these two aggressive phenotypes. In addition to gonadotropin releasing hormone (GnRH), prolactin and somatotropin, male J. transcriptus expressed four genes with no functional information that were also upregulated in dominant A. burtoni. Two of these were part of the dominance module also being upregulated in J. marlieri females. Genes that are associated with multiple aggressive phenotypes are likely to be important in aggressive behavior and will be especially exciting to pursue in future investigations. These genes may have key regulatory roles in the behavior of cichlids and other organisms but, because they have not been previously characterized, they have not been formally investigated. The results of this analysis reiterate how many genes of ecological importance are unknown, and how microarray studies can contribute to providing tentative annotations for these genes (Landry and Aubin-Horth, 2007).

It is interesting to note that J. marlieri females share more gene expression with A. burtoni territorial males than the male aggressive phenotype in J. transcriptus. This latter fact leads to little overlap between the Julidochromis dominance module and genes upregulated in territorial A. burtoni males. One interesting possibility to explain the patterns of overlap between these modules involves the temporal control of territoriality in A. burtoni. Male A. burtoni cycle between territorial and non-territorial states, requiring behavioral plasticity in aggression and dominance. Significant overlap between A. burtoni males and J. marlieri females may suggest that aggression in J. marlieri females was initially a plastic response, more similar to that seen in A. burtoni males, that has become fixed over evolutionary time. In contrast, ancestral male dominance as seen in J. transcriptus may be controlled by different mechanisms.

The phylogenetic relationship of species in the genus Julidochromis has been the subject of recent debate. While Julidochromis was previously considered to be a monophyletic group, recent phylogenetic hypotheses, based on both mitochondrial and nuclear sequence data, suggest paraphyly (Sturmbauer et al., 2010) as well as hybridization leading to the preservation of ancestral alleles. While the exact evolutionary relationship of Julidochromis species remains unresolved the most recent phylogeny suggests that J. transcriptus and J. marlieri are ∼3.7 million years diverged (Sturmbauer et al., 2010) but it remains clear that female aggression in J. marlieri is the derived characteristic, raising the issue of sex-role reversal. Aggressive and territorial behaviors in biparental species function to defend valuable territories, protect young from predators and prevent extra-pair copulations. Male territory defense in most cichlids may be adaptive due to the fact that sexual size dimorphism is typically male biased (Erlandsson and Ribbinck, 1997). Increased female size in J. marlieri could have evolved as a result of increased sexual selection on females (Blanckenhorn, 2005) and female aggression may have been a facultative response. Regardless of the particular driving forces, it is clear that J. marlieri has undergone significant changes in sex-specific behavior since J. marlieri and J. transcriptus diverged. Our results, which show limited conserved sex-specific gene regulation between J. transcriptus and J. marlieri, support the conclusion that major changes in sex-specific gene expression have accompanied this behavioral divergence.

Despite divergence in sex-specific gene expression, there was significant overlap in phenotype-related genes between species, particularly those associated with aggressive behavior. Our identification of a core set of genes shared among aggressive phenotypes fits with the recent evidence leading to the theory of a ‘genetic toolkit’ for social behavior. This theory, similar to ideas concerning the evolution of development, suggests that there may be a small set of conserved genes specialized or co-opted to produce certain behaviors (Toth et al., 2007). Complex behaviors, such as social behavior, could theoretically be constructed from altering this underlying set of behaviorally specialized genes. Evidence for this concept has been found in comparisons of a similar but independently evolved foraging behavior in honeybees and paper wasps that are ∼100 million years diverged (Toth et al., 2010). Researchers found that relatively few genes showed phenotype-specific expression across species but suggested that these genes may be major regulators of behavior. Although J. transcriptus and J. marlieri are more closely related than the taxa discussed above, the core set of genes strongly associated with aggressive behavior in J. transcriptus and J. marlieri conforms to the genetic toolkit theory. The fact that the aggressive phenotype in J. marlieri females is independently evolved reinforces the concept of a genetic toolkit for aggressive behavior, and presents a candidate group of genes important in regulating aggression in Tanganyikan cichlids. Although comparisons with A. burtoni reveal that this genetic toolkit may not include all Tanganyikan cichlids, such comparisons are complicated by aggressive plasticity in A. burtoni. Further comparative studies and additional time points will be needed in order to determine whether the identified core set of aggression genes is in fact used consistently to regulate aggressive behavior in a range of species.

Another theme common in developmental biology that can influence the study of behavior is the concept that changes in traits can often result from changes not only in genes but in gene expression (Carroll, 2000; Toth et al., 2007) (but see Hoekstra and Coyne, 2007). Given that there are species-level differences between J. transcriptus and J. marlieri in aggressive behavior it will be interesting to directly compare levels of aggression-related genes between Julidochromis species.

While J. transcriptus follows the ancestral behavioral pattern of male-biased dominance and female-biased parental care behavior, J. marlieri is reversed in these sex-specific behaviors. We predicted that there would be a high degree of overlap in sex-specific and phenotype-specific gene expression between J. marlieri and J. transcriptus due to their relatedness and similarities in phenotype-specific behavior. However, our investigation revealed that there has been significant divergence in gene expression patterns between J. marlieri and J. transcriptus, perhaps due to the reversal of sex-biased behavior in J. marlieri. The degree to which these differences orchestrate sex-specific behavior as opposed to other divergent phenotypic traits needs to be addressed in future functional work. While most comparisons revealed small overlap between species, overlap between the dominant phenotypes was larger and highly significant. Closer examination of gene expression patterns (Fig. 4) shows that a substantial proportion of the most highly regulated genes in J. transcriptus males and J. marlieri females are part of the dominance module. The large overlap in dominance genes between J. transcriptus and J. marlieri suggests that the gene expression patterns underlying aggressive behavior are shared between the two species, raising the possibility of a genetic toolkit for aggressive behavior in Tanganyikan cichlids. In addition, we found that a number of candidate genes previously associated with aggressive behavior were strongly associated with the aggressive phenotypes in each species. We further provided a tentative functional annotation for a number of unknown genes that are associated with aggressive phenotypes. Investigating the roles of the unknown genes in these phenotype-specific modules is a crucial next step in understanding the mechanisms of aggressive, subordinate and sex-specific behavior in Tanganyikan cichlids.

 
• AVT

  arginine vasotocin

 
• GABA

  γ-Aminobutyric acid

 
• GAPDH

  glyceraldehyde 3-phosphate dehydrogenase

 
• GEL50

  gene expression level 50

 
• GHRH

  growth hormone-releasing hormone

 
• GnRH

  gonadotrophin-releasing hormone

 
• GSI

  gonadosomatic index

 
• LIMMA

  Linear Models for Microarray normalization data

 
• MAP-1B

  microtubule-associated protein-1B

 
• RT-qPCR

  real-time quantitative PCR

We thank Victoria Zero, Kelsey Wood and Heather Machado for insightful discussion on the design and execution of the project. They, and members of the Aubin-Horth lab, provided comment on earlier drafts of this manuscript.