Introduction

For most animals that possess sex chromosomes, there also exists some global mechanism to balance out the effects of gene dosage in males and females. For example, even though mammalian females have two X chromosomes and males only one, the vast majority of genes on the mammalian X show statistically identical expression levels in both sexes. This is achieved through both the inactivation of a single X chromosome in females, which balances the X dosage in the sexes, and the subsequent upregulation of the active X chromosome in both sexes to balance the overall expression level of X-linked genes with autosomal regions (Lyon, 1999; Nguyen and Disteche, 2006; Lin et al., 2007; Johnston et al., 2008). Similarly, Drosophila species equalize gene expression on the X chromosome by unilaterally upregulating X-chromosome expression in males, which simultaneously equalizes dosage effects between males and females, and between the single male X and the diploid autosomes (Lucchesi et al., 2005; Gupta et al., 2006). Caenorhabditis also has a unique dosage compensating mechanism, upregulating X-linked genes in both sexes to achieve X:autosome parity for males, and downregulating the two (now overexpressed) X chromosomes in hermaphrodites to balance the X and autosomes (Meyer and Casson, 1986; Meyer et al., 2004).

The Drosophila, Caenorhabditis, and mammalian X chromosomes, though they share similar evolutionary patterns and pressures, have arisen independently (Vicoso and Charlesworth, 2006), so their associated dosage compensating mechanisms are independent as well. This repeated origin of global, complex methods of dosage compensation suggests that it is important for lineages with sex chromosomes to evolve some mechanism to balance out sex chromosome dose in gene transcription. This is probably due to the complex and interactive nature of gene expression, where multiple components form pathways and networks, often creating epistatic and complex connections (Lynch, 2007). Changes in the expression levels of one or a few genes can affect dozens or hundreds of others through direct or indirect interactions (Boone et al., 2007), and these pathway perturbations could have serious effects on overall viability and survival. In animals with large differences in the gene content between the X and Y, dosage could potentially affect hundreds of genes. In humans, most large aneuploidy and trisomy events are lethal, and of those that can be tolerated by the developing embryo, many produce severe congenital birth defects (FitzPatrick, 2005; Hassold et al., 2007) or sterility (Doswell et al., 2006; Visootsak and Graham, 2006). It is consequently logical that as the degenerative process of the Y chromosome reduces its gene content and makes a larger proportion of X-linked genes subject to gene dosage effects (Charlesworth and Charlesworth, 2000; Charlesworth et al., 2005), some mechanism evolves to balance the dosage of the X chromosome complement in both sexes.

It was therefore surprising when recent studies showed evidence that birds seem to lack an overall sex chromosome dosage compensating mechanism (Ellegren et al., 2007; Itoh et al., 2007). In birds, the female is the heterogametic sex, with one Z and one W chromosome, and the male is homogametic with two Z chromosomes. The lack of dosage compensation of the chicken Z chromosome produces a pronounced overall pattern of male bias for most genes along the Z, which is in marked contrast to the autosomes, which, when averaged over all loci, have a roughly equal expression level between males and females (Ellegren et al., 2007; Itoh et al., 2007). The chicken Z chromosome is large, and contains more than 700 known or predicted protein coding genes, and is greatly diverged from the W in terms of gene content (www.ensembl.org). This suggests that the lack of dosage compensation directly affects a large proportion of the genome, and this becomes even larger when the network structure of the transcriptome is considered.

There are several potential explanations for the lack of dosage compensation in birds, none of which has yet gained consensus. First, it may be that genes are compensated locally on a gene-by-gene and tissue-by-tissue basis. This would mean that critical genes would be individually regulated to a similar expression level in both sexes, and would result in a highly variable local sex-bias pattern across the Z chromosome that would vary somewhat among tissues. There is some support for this, with global transcription studies of sex-biased expression in birds recovering an overall male bias along the Z chromosome, with a wide variance in sex bias locally (Ellegren et al., 2007; Itoh et al., 2007; Melamed and Arnold, 2007). It is also possible that sex chromosome dosage occurs at critical developmental time points (Wilkins, 2002), or shows some other temporal dynamic such as that has been recently noted in mice (Lin et al., 2007). This would mean that the overall level of sex bias on the Z chromosome would differ at different points in the life cycle. There is not much evidence of this in birds, as Z-linked genes show persistent sex bias in early (Scholz et al., 2006) as well as late (Ellegren et al., 2007; Itoh et al., 2007) embryonic development. However, this is a difficult theory to reject at this point, as it requires a unified time-series transcription profile stretching through embryogenesis into adulthood. Finally, a nascent dosage compensating mechanism might be present, but affect only a small crucial portion of the Z chromosome, not having yet succeeded in spreading across the chromosome. This would produce a limited neighborhood on the Z chromosome where gene expression is similar in both sexes across tissues throughout the life cycle. Some evidence for this has been shown (Melamed and Arnold, 2007), but is not yet conclusive.

To disentangle these alternative hypotheses, we created a time-series transcription profile for chicken Z-linked genes in both brain and gonad tissue stretching from embryogenesis to sexual maturity. Furthermore, in order to differentiate sex-biased gene expression from potentially incomplete dosage compensation, we employ two different measures of expression difference. We first calculate fold-change difference in gene expression between males and females, which indicates how much and in which direction transcription varies between the sexes. Fold-change can be problematic in assessing dosage compensation above the level of individual genes however. Using fold-change metrics alone, strong but conflicting sex bias could be mistaken for neighborhoods of dosage compensation when averaging across genes in a defined region. We therefore used an additional metric of sex-bias amplitude in order to differentiate contradictory sex-biased gene expression from potential dosage compensation. Amplitude is only a measure of the degree to which genes are sex biased, but ignores the direction of bias, therefore dosage compensation can be inferred in regions where both amplitude and fold-change approach zero. Using this methodology, we demonstrate that there is no local region of nascent dosage compensation on the chicken Z, but rather gene expression is regulated on an individual level that varies between embryos and adults and by tissue type.

Materials and methods

All details about sample collection, RNA preparation, microarray analysis, quality control, will be described elsewhere (Mank et al., in preparation) and details regarding these methods are available in Supplementary material S1.

Statistical analysis

We calculated sex bias by determining the log2 value of the average male:average female expression level for each probe-set that mapped to the Z or autosomes based on the ENSEMBL map of the chicken genome (WASHUC2, assembled in May 2006). This produces a positive value for loci that are expressed more in males, a negative value for genes expressed more in females and a value near zero for regions that are similarly transcribed in the sexes. We computed the average when multiple probe-sets mapped to the same gene. This resulted in expression data for 635 genes across the Z chromosome. We calculated average fold-change for all Z-linked genes for which we had expression data, and estimated 95% confidence intervals with bootstrapping (1000 repetitions). Statistical difference between adult and embryonic time points within each tissue was determined with t-testing (two-tailed, assuming unequal variance between samples).

We also calculated an amplitude metric, which is the absolute value of the log2 male:female (m/f) expression ratio. The metric effectively estimates degree of expression bias without regard to directionality (male- versus female-bias), and the higher the amplitude value, the larger the degree of sex bias. We calculated average amplitude over the Z chromosome, and 95% confidence intervals and statistical difference between adult and embryonic time points were determined in the same fashion as fold-change.

We then drew an expression map of the Z chromosome by plotting fold-change as a function of location for each gene individually. We averaged log2 fold-change over adjacent genes, using 3 Mb overlapping windows with a 1Mb shift. This type of analysis could potentially create perceived regions of dosage compensation due to a few genes with strong female bias. To differentiate this situation from true regions of dosage compensation, we also calculated average amplitude, as this is a measure of the overall degree of sex bias (both male and female). Theoretically, if a neighborhood of dosage compensation exists on the Z chromosome, the overall level of sex bias would decrease, and so the average amplitude of that region would approach zero. This metric therefore differentiates between a strongly female-biased gene or set of genes pulling the average sex bias of a region down, which is not evidence of dosage compensation, and reduced overall sex bias in a region, which would be potential evidence of neighborhood dosage compensation. We determined the 95% confidence intervals for these 3 Mb overlapping windows by permuting the underlying fold-change data sets for each tissue/time-point treatment (1000 repetitions).

Gene Ontology

To obtain some understanding of the types of genes that are dosage compensated on the chicken Z chromosome, we calculated statistical overrepresentation of Gene Ontology (GO; Gene Ontology Consortium, 2000) terminology with Ontologizer 2.0 (Robinson et al., 2004), using Parent-Child analysis with the Westfall-Young-Single-Step method to correct for multiple comparisons (Grossman et al., 2007). However, because of the large number of comparisons that these analyses entail, we report all overrepresentations that are significant (P<0.01) in the test before multiple comparison correction. For the gene population set, we used all significantly expressed genes mapped to the Z chromosome. To increase the power of these tests, and to get a generalized pattern of functionality in all our treatments, we only looked at sex in our GO analysis, and removed both tissue and age from our identification of genes that were unbiased, male biased and female biased. This means that we assigned genes to lists based on the overall expression pattern throughout the organism and its lifetime, rather than searching for fine-scale patterning that, given the current gaps in the GO annotation, would likely be meaningless. Our analysis was based on three study lists: unbiased genes (genes for which the difference between male and female expression in all tissue/time points was less than 10%); male-biased genes; female-biased genes. Both male- and female-biased genes were identified with a 1.5 fold-change cutoff and a 10% FDR correction.

Results and discussion

We collected gene expression data for male and female brain at embryonic day (ed) 10, ed15, ed19, and adults, and for the gonad of both sexes at the later three time points. Previous analysis of cDNA array profiling of brain tissue during earlier chicken development indicated that the expression pattern for most Z-linked genes does not vary greatly, and found no evidence of early dosage compensation (Scholz et al., 2006), therefore we focused on late embryonic development (ed10, ed15, and ed19) and adulthood.

We calculated several expression metrics at each tissue/time-point treatment in order to create a dynamic picture of gene expression patterns in the gonad, which is the most sex-biased tissue, and the brain, which is generally among the least (Parisi et al., 2004; Yang et al., 2006). This allows us to study temporal, regional, and local patterns of sex bias along the roughly 630 Z-linked genes with significant expression in our data, and to test alternative theories of dosage compensation.

Age and sex were both strongly correlated to Z chromosome expression patterns in our data, regardless of tissue type (ANOVA, P=1.94 e−16 for age, 2.7 e−15 for sex). When averaged over the length of the chromosome and over age categories, Z-linked genes show a pronounced male bias in expression pattern, with a mean 1.24-fold higher expression in males than in females (log2 m/f=0.308) for brain and 1.32-fold higher in males than in females (log2 m/f=0.396) for gonad. This is not due to an overall male bias in chicken gene expression as average fold-change for autosomal genes does not differ statistically from zero in any tissue or time point. Because males possess two copies of the Z whereas females only possess one, Z-linked gene transcription can be expected to be higher in males than females for all genes that do not have other factors acting on their rate of transcription. Although it may initially seem logical that a twofold change in copy number results in a doubling of transcription, copy number variation studies indicate that gene duplication actually produces a less than twofold increase in expression, with an average of about 1.5-fold change in expression seen in studies of the human genome (Hyman et al., 2002; Pollack et al., 2002; Yao et al., 2006). Possibly, this is due to network buffering (Oliver, 2007). From this, we can tentatively conclude that genes on the Z with this level of fold-change between males and females are that way due simply to dosage effects alone.

The fact that the Z chromosome m/f ratio average was significantly less than a twofold dosage effect of transcriptional targets alone may also attests to the power of selection acting on some genes to either balance male and female expression (in effect a local gene-by-gene dosage compensation), and on other genes that are antagonistic in males to produce female-biased expression (Connallon and Knowles, 2005). Both of these processes, if acting on a significant number of genes, would bring the overall log2 m/f ratio of the Z chromosome down. Indeed, many genes showed either a balanced or female-biased expression pattern at some tissue/time-point treatment.

Fold-change varied between the brain and gonad within time points (Figure 1), and was significantly different between tissues for ed15 (two-tailed t-test assuming unequal variance, P=0.002) and ed19 (P<0.0001), but not adults (P>0.1). These results suggest that dosage regulation of Z-linked genes has both local and tissue-specific aspects, with genes being regulated regarding sex-bias based on their essentiality and tissue-specific functions. We next examined our data for both temporal and neighborhood effects of dosage compensation in brain and gonad.

Figure 1
figure 1

Average Z chromosome log2 m/f fold-change (black) and amplitude (white) for both gonad (a) and brain (b). Whiskers indicate 95% confidence intervals based on bootstrapping (1000 repetitions).

Temporal aspects of sex-bias and dosage compensation

The overall pattern of gene expression shows significant variation in expression between sexes and among time points. Adult fold-change did not differ from embryonic time points in either tissue when averaged over all Z genes (Table 1; two-tailed t-tests), however, there was a temporal pattern in the data when assessed by amplitude of expression bias. We used amplitude, which is a measure of sex bias that is not subject to directionality, to identify the difference between regions of conflicting sex bias and regions lacking sex bias. Amplitude, as with fold-change, was more pronounced in gonad than in brain data. Unlike fold-change however, t-testing indicated that amplitude was significantly greater in adult brain tissue (P⩽0.01 in all comparisons) and adult gonad (P⩽0.001 in all comparisons) than embryonic time points in the same tissue when averaged across all significantly expressed Z-linked genes. This suggests that the degree of both female and male sex bias is elevated in both brain and gonad during adulthood.

Table 1 Z chromosome average fold change and amplitude

There are several reasons why developing embryos might need to have more effect dosage compensating mechanisms than adults. This could happen if embryos were less tolerant of perturbations in gene networks due to Z chromosome dosage differences in males and females (Wilkins, 2002). However, the overall pattern in our data suggests that is not the case. Fold-change is statistically consistent across all time points within tissues. This, in conjunction with the increase in amplitude, indicates that sex-biased expression increases in adults as genes become more male- and female biased, but not that Z chromosome dosage is equalized in embryos. For the latter situation to be true, the data would have to show a decrease in both fold-change and amplitude in embryos when compared to adults. An increase in both male and female sex bias would theoretically follow maturity as the production of sex-specific gametes and behaviors associated with reproduction would require an increase in sexually dimorphic gene expression.

Neighborhood effects of sex-biased expression

Gene-by-gene fold-change varied over the length of the Z chromosome (Figures 2 and 3), with a strong dip around 28 Mb. This region has been previously identified as abnormal with regards to sex-biased gene expression (Melamed and Arnold, 2007), where it was hypothesized to be a neighborhood of nascent dosage compensation. To verify this finding, and to look for other neighborhoods of possible dosage compensation, we computed the moving average for both fold-change and amplitude for each tissue/time-point treatment using a 3 Mb window and a 1 Mb shift. Although the average fold-change for the windows around 28 Mb approaches zero, our amplitude data suggest that this is not in fact a region of dosage compensation. Theoretically, if dosage compensation were equalizing the gene expression levels between the sexes in this region, both amplitude and fold-change would approach zero. However, this is not the case, as the amplitude does not dip perceptibly in any of the analyses.

Figure 2
figure 2

Fold-change map of the chicken Z chromosome for ed10 (a), ed15 (b), ed19 (c) and adult (d) brain samples; log2 m/f ratios are plotted for each gene. Black lines represent the average fold-change for overlapping windows, grey lines represent the same for amplitude (both fold-change and amplitude windows were 3 Mb with a 1 Mb shift). The fold-change and amplitude moving average are often overlapping, in which case the amplitude moving average is plotted in front.

Figure 3
figure 3

Fold-change map of the chicken Z chromosome for ed15 (a), ed19 (b) and adult (c) gonad samples; log2 m/f ratios are plotted for each gene. Black lines represent the average fold-change for overlapping windows, grey lines represent the moving average for amplitude (both fold-change and amplitude windows were 3 Mb with a 1 Mb shift). The fold-change and amplitude moving average are often overlapping, in which case the amplitude moving average is plotted in front.

We statistically verified that fold-change, but not amplitude varies in this neighborhood by permuting the individual expression data for Z-linked genes to achieve confidence intervals for overlapping windows. The region 26–29 Mb falls below the 95% confidence range for fold-change for all brain and the embryonic gonad samples. However, this region was not statistically different regarding overlapping amplitude windows. Furthermore, the removal of the two strongly female-biased genes in this region, which are not immediately adjacent to each other, and which are described in more detail below, resulted in fold-change windows that were statistically identical to the remainder of the Z chromosome. This further supports the notion that the region is not dosage compensated, but rather contains some genes of strong female bias.

This permutation procedure identified another unusual region at 58–62 Mb. This area has an unusually high fold-change that is not due to chance (P<0.05) for brain ed10 (average log2 m/f fold-change for genes in this region=0.478) and ed15 (log2 m/f=0.464), and gonad ed15 (log2 m/f=0.468) and ed19 (log2 m/f=0.6248), indicating that the overall male bias of this area is greater than the Z chromosome as a whole. This agrees with previous observations (Melamed and Arnold, 2007), and could be taken as evidence that sex bias may be regulated in large genomic blocks, which has previously been shown at a smaller scale for sex-biased genes (Mank et al., 2008). However, although corresponding amplitude was higher than expected by chance (P<0.05) for the brain ed10 and ed15 samples, closer examination of the fold-change maps (Figures 2 and 3) revealed that these regions do not appear to have an excess of strongly male-biased genes, but rather a deficit of genes that show equalized male and female expression levels.

Female-biased genes in a male-biased world

The two strongly female-biased genes around 28 Mb (Table 2) are clearly visible in the brain panels (Figure 2), but also discernible to a lesser extent in the gonad samples (Figure 3). The first gene, mapped to position 27.96 Mb, is female biased in all tissue time points with expression at least twofold greater in females. This gene bears no homology to known or predicted transcripts in other organisms. The second gene, Pro-relaxin 3, maps to position 28.16 Mb of the Z chromosome, and is at least 1.5-fold female biased in all time points except ed15 in the gonad, when it is significantly male biased. Pro-relaxin 3 is a member of the relaxin protein family, which originated early in vertebrate evolution (Wilkinson and Bathgate, 2007). Although it is known to be involved in mammalian parturition, it is expressed in many somatic tissues as well (Ma and Gundlach, 2007). Its function in avian reproduction and sexual differentiation is as yet uncharacterized.

Table 2 Genes with consistent female bias on the Z chromosome

In all cases where these genes were significantly female biased, it was due to a reduction in the male expression level, when compared to the average male expression level across all Z-linked genes. This could potentially indicate that these genes are incorrectly assigned to the Z chromosome, and are in fact located on the W chromosome, which would result in female-biased or female-limited expression. We therefore PCR verified that both female-biased genes are genetically present in the DNA of five male and five female samples (data not shown), thereby eliminating the possibility that they are W-linked genes.

The two unusual genes are located within the male hypermethylated region, and near the DMRT1 locus (at 26.7 Mb), which likely involves in sex determination in birds (Shetty et al., 2002; Smith et al., 2007), and is male biased in our embryonic gonad samples. Although the female-biased genes have no known role in avian sex determination, they must be important in the formation of avian sexual dimorphism, as they are so strongly and consistently female biased through so many time points, and especially in light of the overall male bias of the Z chromosome. They may have an important function downstream in sexual differentiation, and merit further genetic investigation.

Gene Ontology

The GO is a curated database of gene functionality (Gene Ontology Consortium, 2000), and can be used to compare the attributes of different gene lists. To assess whether the expression characteristics of Z-linked genes are correlated to unusual gene functions, we used the GO to assess overrepresented terminologies in three gene lists (155 unbiased genes, 166 male-biased genes and 6 female-biased genes, see ‘Materials and methods’ for criteria that were used to define these lists), each compared to all expressed Z-linked genes. These lists were designed to maximize the presence of sexually dimorphic genes, and so in establishing them, we averaged across all tissue/time-point comparisons.

Both the unbiased and male-biased gene lists returned overrepresented GO terms, as shown in Table 3. These results should be interpreted with some caution, as there is a considerable danger with tests of GO terms to falling into the post hoc ergo propter hoc fallacy, and finding plausible but otherwise unsupported explanations of gene expression characters from the overrepresented functions. In addition, only 28% of the probe-sets on the Affymetrix Chicken GeneChip that have been mapped to chromosomes 1–28 or the Z in the WASHUC2 genome build have GO annotations as of January 2008, leaving much opportunity for omission in interpretation (Rhee et al., 2008). Using the human GO annotation could yield higher coverage, as 54% of the mapped Affymetrix chicken probe-sets have known human orthologs with corresponding GO annotations. However, the hundreds of millions of years separating birds and mammals require an assumption of conserved functionality that seems tenuous at best, therefore we chose to the chicken GO annotation.

Table 3 Over represented Gene Ontology terms

It is reasonable to expect the unbiased gene list to be enriched for genes with important functions that cannot tolerate dosage differences between males and females. It is difficult to know if the most significantly overrepresented term from this list, molecular transducer activity, fits this description. Of the genes that were male biased in both gonad and brain tissue, there were several significant terms, with pore complex and DNA metabolic process both significant after multiple-testing correction. There were no overrepresented terms from the female-biased list, most likely because there were so few female-biased genes on the Z.

Conclusions

We examined a time-series global gene transcription profile in chicken to test whether there were regional (hypothesis I), developmental (hypothesis II) or only local gene-by-gene (hypothesis III) dosage compensating mechanisms for the Z chromosome. We found no evidence of a regional (I) (Melamed and Arnold, 2007) or developmental (II) (Lin et al., 2007) aspects to dosage compensation. Rather, our data suggest that dosage compensation is entirely a local process (III), performed in a gene-by-gene basis independently for each tissue or developmental time point. This has important implications both to the nature of gene regulation and the evolution of sex chromosomes. Specifically, birds are able to tolerate a mild male bias for many Z-linked genes, and only equalize transcription levels for some genes, presumably when and where it is so required. How this gene-by-gene regulation is carried out remains an interesting area for further research, though some insight into possible mechanisms can be gleaned from work on aneuploid mutants. For some loci, trisomic or monosomic variants do not exhibit noticeable dosage effects (Birchler, 1979; Guo and Birchler, 1994), possibly due to network buffering or transcriptional regulators (Birchler et al., 2001). In addition, not all loci are equally sensitive to dosage effects, as some aneuploid mutants with observable dosage effects have no observable phenotype (Veitia et al., 2008). If aneuploidy is a fair analogy for sex chromosome dosage differences between the sexes, this would suggest that those genes on the Z chromosome with balanced transcription are associated with autosomal regulatory factors, whereas unbalanced Z-linked genes may have different, less dosage-sensitive functions.

All this means that sex chromosome evolution need not require a concomitant wholesale mechanism to equalize gene dosage, which raises new questions as to why global dosage compensation is present in other lineages.