A greater number of male-enriched genes
Our results are consistent with the predictions of the "male sex drive" hypothesis. Three lines of evidence from our study provide reason to believe that gene expression in the zebrafish lineage is "masculinized." First, we discovered a larger total number of male-enriched than female-enriched genes (Table 2), consistent with other animal studies. A recent study, for example, documented this asymmetry in five Drosophila species (D. melanogaster, D. simulans, D. yakuba, D. ananasse, and D. virilis) using species-specific microarrays , and additional investigations have reported similar findings in Drosophila [11, 19, 40]. Rin et al. also identified a substantially greater number of male-enriched genes, especially within higher fold change classes, based on a transcriptomic comparison of testis and ovary in mice . In two closely related frog species (Xenopus laevis and X.muelleri), Malone et al. revealed a greater overall number of male-enriched genes and demonstrated an even more pronounced male-biased asymmetry among genes that are also differentially expressed between species . Indeed, others have described a related phenomenon, in which male-enriched genes are greatly overrepresented among groups of genes that demonstrate intra- and inter-specific expression polymorphism, relative to female-enriched and sex-unbiased loci [44–46].
Interestingly, recent studies of sex-biased gene expression in Danio rerio have not yielded the same observation of more male-enriched than female-enriched transcripts. In fact, Santos et al. compared ovary and testis transcriptomes in adult zebrafish and reported 1370 male-enriched genes and 1570 female-enriched genes , which contrasts with our finding that more genes are male-enriched. One possible source of the discrepancy might be that the experimental animals were treated quite differently in our study. Santos et al. sampled individuals from a "breeding colony" of six males and six females, and histological analysis of experimental ovaries revealed great variation in oogenic stage among individual females . Females in our study spawned on the same day, and were then isolated from males for five days before being sacrificed. Separation of males and females may not reflect conditions zebrafish experience in nature, but our design allowed us to prevent re-mating and standardize reproductive cycles among experimental individuals. Still, a five-day absence of any stimuli produced by the opposite sex might result in significant behavioral and physiological consequences for males and females, and these could explain the differences between the studies. For example, significant changes in gene expression over a very short time period as a consequence of courtship exposure have been documented in Drosophila . Additional studies should be conducted to assess the potential for plasticity of sex biases in the transcriptome due to behavioral, environmental, developmental [48, 49], and temporal factors.
Differences in array platform and analysis might also explain the discrepancy between studies. Santos et al. employed microarrays constructed from the Sigma-Genosys (Cambridge, UK) Zebrafish OligoLibrary™, which represents approximately the same number of unique transcripts (15,806) as the Affymetrix arrays (14,900), but not necessarily the same transcripts. Furthermore, the expression detection algorithms tailored for Affymetrix GeneChips® are unique, and we applied four of these in this study. It is also worth noting that the microarray fold change estimates from the Santos et al. study are substantially lower (up to 2 orders of magnitude) than the corresponding real-time qPCR fold change estimates, which the authors attribute to spot saturation . Our microarray fold change estimates appear to be more consistent with the real-time qPCR estimates (Table 6), suggesting that array feature saturation is less of a problem in our study. Despite the discrepancy, however, there is agreement between the two studies at the level of expression patterns for individual genes, as nine out of ten top sex-biased genes identified by Santos et al.  also appear in our sex-biased gene list (Table 3).
Two other studies addressed sex-biased gene expression in zebrafish, but neither of them is as relevant to this study as the Santos et al. experiment. Wen et al. conducted a whole body male-female comparison of the zebrafish transcriptome using a cDNA microarray representing 8793 unique EST clusters . The authors identified 383 female-enriched genes in their study; however, they make no mention of male-enriched transcripts, and gonads were not analyzed separately. Another microarray study, by Sreenivasan et al., did separate the gonads, in addition to the brain and kidney, from the "rest-of-body," for males and females . They employed cDNA microarrays containing 6370 unique genes derived from zebrafish gonad EST libraries. Sreenivasan et al. reported 881 genes enriched by ≥ 1.5 fold in the testis relative to the common reference control, and 1366 genes enriched by ≥ 1.5 fold in the ovary relative to the common reference control . The report does not provide details regarding the total numbers of male- and female-enriched genes for each organ comparison, so a direct comparison between this study and ours is difficult.
Another surprising result is that we did not identify genes that, according to our strict consensus criteria, demonstrate sex-biased expression at the level of the zebrafish body. A recent study of sex differences with respect to hepatic gene expression, which also utilized the Affymetrix platform, revealed 1249 sex-biased genes (792 male-enriched, 650 female-enriched) in the adult zebrafish liver . Another study, which examined sex differences of the zebrafish brain transcriptome, identified 42 sex-biased genes (18 male-enriched, 24 female-enriched) . This is in stark contrast to Sreenivasan et al. , who report 3080 genes as differentially expressed between male and female brains, so it is clear that major differences exist among the other zebrafish studies as well. Our study did not involve a direct organ-to-organ comparison (except for gonads), so it is possible that organ-specific signals of sex-biased gene expression were obscured by background gene expression in other somatic tissues. The lack of sexually dimorphic body gene expression in our study could also be a consequence of high among-individual variance in body gene expression, although we took many steps experimentally to reduce this. Furthermore, our statistical criteria for differentially expressed genes were very conservative, so we likely missed some differentially expressed genes, especially if the differences were small. If we relax our criteria and consider a gene differentially expressed if it appears significant in at least one of the four absolute expression comparisons, then we find 112 body sex-biased genes (78 male-enriched, 34 female-enriched). Of these genes, 26 (9 male-enriched, 17 female-enriched) were consistent with the liver results from Robison et al. , but none were consistent with the brain study . The list of 112 genes, and corresponding fold change estimates from the four absolute expression comparisons are included as Additional file 3.
More gonad-soma differences in males
A third result of our study related to reproductive processes and sex-specific gene expression patterns is simply that adult male zebrafish demonstrated many more gonad-soma differences in transcript abundance than females. We detected 5340 genes as differentially expressed between testicular and male body tissue (3002 testis-upregulated, 2338 testis-downregulated). In comparison, only 2380 genes were identified as being differentially expressed between ovarian and female body tissue (981 ovary-upregulated, 1399 ovary-downregulated). These striking transcriptional differences at a tissue-specific level are likely reflections of fundamental reproductive differences between males and females. A microarray study of D. melanogaster adults revealed a similar sex disparity in gonad-biased gene expression and also reported that the expression magnitude of testis-upregulated genes is substantially greater than that for ovary-upregulated genes . Because none of the 981 ovary-upregulated genes identified in our study demonstrated fold change values greater than four, whereas fold change values for 554 testis-enriched genes exceeded six, zebrafish may also conform to this pattern. A general interpretation of this trend might be that there are more specific transcripts essential to processes that take place in the testes, relative to specific transcripts in ovarian tissue.
A small comparison of testis-upregulated or testis-specific genes from other zebrafish studies [28, 42] to those identified as testis-upregulated in our study indicates a high level of agreement (see "testis-upregulated" section of Table 3). In contrast, many of the top ovary-specific or ovary-upregulated genes identified consistently in these studies are absent from our list of top ovary-upregulated genes (Table 5). Why our study differs from the others in this respect remains an open question. Again, the fact that we separated males from females five days prior to sample collection may partially explain the discrepancy, especially if females experience major changes in hormone profiles in the absence of males. High body gene expression variance among females in our samples could also explain why ovary-upregulated genes from the other studies did not demonstrate statistically different expression levels in our study. Additional file 7, a more detailed version of Table 3, includes ten reportedly ovary-upregulated genes and the relevant expression value means, standard errors, and fold change estimates from our data set.
A particularly important class of female reproductive genes, which correspond to members of the zona pellucida egg coat glycoprotein superfamily, demonstrate ovary-specific expression patterns according to several zebrafish studies (zp1 ; zp2 [43, 50]; zp3 [50, 51]). We, however, identified none of the zona pellucida homologs represented on the zebrafish GeneChip® as significantly ovary-upregulated (See Additional file 8 for a list of zp genes, expression value means, and standard errors for each absolute expression analysis). This result is surprising, and the expression values in Additional file 8 indicate high female body zp expression in addition to expectedly strong expression in ovaries. Contamination of the body sample with ovarian tissue could produce this result but is unlikely since we completely removed all visible ovarian tissue from each individual. Even if a dissection left as much as half of the total ovarian tissue inside a body sample, one would not expect equal or greater body transcript abundance (for a truly ovary-upregulated gene), because the contaminating ovary signal would be greatly diluted by the female body RNA. Furthermore, if the female body samples were contaminated with ovarian tissue, we would expect many false positives with respect to male and female body differences, which is clearly not the case. We, therefore, maintain that high female body zp expression in our experiment is either real or a reflection of problematic zp array probesets. In general, there seems to be some disagreement across studies with respect to tissue specific patterns of zp gene expression. For example, significant expression of zp1 and zp2 has been documented in ovary-excised females , and expression of zp3 in female skeletal muscle has also been described . Furthermore, a recent study (which also used Affymetrix zebrafish arrays) of sex-biased gene expression in the liver of zebrafish reported that zp2.2, zp3, zp3a.1, zp3b, and zpcx are all expressed at high levels and are all female-enriched . Based on an estimate by Liu et al., there are likely 10 - 15 zp2 and 17 - 21 zp3 paralogs alone distributed throughout the zebrafish genome , so assaying expression of individual paralogs may not be as straightforward as is assumed. We cannot say for certain that our results reflect this specific problem, but across-study differences in zp probe composition might explain some of the inconsistencies in tissue-specific expression patterns of zona pellucida genes.