The unique genomic properties of sex-biased genes: Insights from avian microarray data

Expression properties

There were several distinct expression patterns that delineate sex-biased genes from the remainder of the transcriptome. First, sex-biased genes were expressed in fewer tissues than unbiased genes for both autosomal and Z-linked microarray targets (p < 0.00001 for both Z and autosomal genes), as determined by equality of proportions test for overlapping regions (Fig. 1). Additionally, based on normalized florescence levels from the microarrays, sex-biased autosomal genes within each tissue were expressed at consistently lower levels than unbiased genes (Fig. 2), even for the sex with higher expression. Z-linked genes show a slightly different pattern however, as the average expression level across male replicates was statistically the same for Z-linked male-biased and unbiased genes. The average expression across female replicates for male-biased genes was somewhat lower than for unbiased genes.

Expression noise, or the variance in expression among individuals, has been shown to negatively correlate with gene essentiality, or those genes that are associated with the core functions in an organism [30–32]. We therefore calculated a normalized variance estimate for the expression among within-sex replicates as a proxy for expression noise. Sex-biased autosomal gene categories had higher within-sex variance values than unbiased autosomal genes in all three tissues (Fig. 3). There was far less difference in variance between Z-linked unbiased and male-biased gene sets.

Clustering of genes with similar patterns of expression

Genes with similar breadth and level of expression have been shown to cluster in eukaryotic genomes [33–35], probably due to passive co-regulation of neighboring genes. We tested this possibility by calculating the average difference between both normalized level of expression between neighboring autosomal genes for which we had expression data, taking data from each tissue separately. We calculated the significance of local similarities by comparing these statistics with those derived from 10000 replicates where gene position was randomized on each chromosome. Neighboring genes showed significantly correlated levels of expression in gonad (p = 0.0002) and heart (p = 0.0114) although not in brain (p = 0.1848).

We also looked for clustering of sex-biased genes. Fold-change in expression, based only upon difference in expression between the sexes without a statistical cut-off p-value to define sex-bias, was significantly correlated in gonad (p < 0.0001), heart (p = 0.0175) and brain (p = 0.0131). All these correlations remained significant at the same level upon removal of comparisons of neighboring paralogs. This suggests that tandem duplications are not responsible for the observed local similarities, and that genes with similar expression patterns for some tissues are more likely than chance to be proximately located.

We also examined the genome for neighboring genes that both had significant evidence for biased expression in the same sex using the absolute value of log₂ fold change > 1 with p _adj < 0.05 to denote sex-bias. There were 157 such pairs using expression in gonad, compared with 100.9 expected by chance. When tandem duplicates were excluded, we observed 144 neighbor pairs, compared with 97.8 expected by chance alone. Both these associations are statistically significant (p < 0.0001 in both cases). However, when brain or heart expression was considered, there were no such examples of neighboring sex-biased genes.

GC3 and codon bias

Third position GC content (GC3) has been linked to codon bias [36], which has in turn been demonstrated for different classes of sex-biased genes in D. melanogaster [14]. We found that GC3 differed by expression class for autosomal heart and gonad genes (one-way ANOVA, 2 d.f., p < 0.0001 in both cases). For the heart, female-biased genes had lower GC3 content than unbiased genes (Tukey's post hoc test, p < 0.01). GC3 for both male-and female-biased genes expressed in the gonad was significantly different than the unbiased gene class (p < 0.01 in both cases). There was no difference in GC3 among Z-linked expression classes (Fig. 4).

Codon usage, as measured by the effective number of codons (ENC) [37] also differed by expression class for autosomal heart and gonad genes (one-way ANOVA, 2 d.f., p < 0.003 in both cases), though not as starkly as GC3 levels (Fig. 5). Post hoc testing indicated that only autosomal female-biased genes expressed in the gonad showed more codon bias than genes with roughly equal expression in males and females.

Figure 6 shows the relationship between GC3 and codon usage bias (ENC) in the genes in our dataset, with a smoothed spline fitted to the data. Codon bias is greatest when GC3 is skewed away from equal usage of GC or AT base pairs. The same relationship is found if the GC content of surrounding noncoding regions is used instead of GC3. When the ENC values of each sex bias category were compared with those predicted by GC3 from the model fitting, no significant differences were found (one-way ANOVA on residuals), suggesting that GC3 values do not predict codon-bias in our sex-biased expression categories.

Gene Ontology

Expression properties

Sex-biased genes were present in both somatic and gonad tissues, though the gonad showed the highest proportion of sex-biased gene expression across the range of fold-change cutoff values. This likely stems from the fact that the gonad is the most sexually dimorphic organ, and although the testes and ovary share a common precursor, in many ways they represent two distinct tissue sets. It is therefore logical that sex-biased gene expression would be the most pronounced in this region.

Sex-biased genes share several expression properties that have implications to their evolution. First, sex-biased genes were more tissue specific in our analysis (Fig. 1), which is consistent with previous work suggesting sex-biased genes are more narrowly expressed than unbiased genes [5, 39]. Coding regions that are broadly expressed face many, sometimes contradictory, evolutionary constraints on gene expression evolution [40], and it may be that genes with extensive expression throughout an organism are less able to evolve sex-biased expression patterns without disrupting critical functions in other tissues.

In addition to being expressed at lower levels (Fig. 2), sex-biased genes also showed more variance in expression (Fig. 3), at least for genes mapped to autosomal locations. Variance is an indicator of expression noise, and genes with larger amounts of expression noise have been shown to be less important in organismal survival and functionality [30–32]. Somatic organs are functionally identical in males and females, therefore sex-biased autosomal genes expressed in these tissues theoretically do not have important developmental or regulatory roles. The reduced difference between autosomal biased and unbiased genes in the gonad may reflect the sex-specific nature of the ovary and testes, where the pattern of sex-biased expression observed in the gonad is required for the formation of these distinct tissues. The lack of dosage compensation of the Z chromosome means that most genes have a default male-biased expression pattern, and this confounds the analysis of male-biased Z-genes and essentiality. The reduced variance observed across all tissues for male-biased Z-linked genes when compared to male-biased autosomal genes likely stems from the default male-biased expression pattern, as male-biased genes may have more important roles in the organism, but still exhibit male-biased expression simply due to Z-linkage. Overall, these findings suggest that somatic sex-biased genes, particularly those attached to the autosomes, may have less critical roles for organismal survival than unbiased genes.

There are indications [6] that avian sex-biased genes show an altered form of the accelerated rates of protein evolution seen in other animals[4, 9], with female biased genes in birds having the highest rate of functional change. The above-described connections between tissue specificity, expression level, and expression variance and sex-bias may offer explanations for this rapid rate of evolution. Rates of evolution have been positively correlated with tissue specificity [41, 42], and negatively correlated with overall expression level [43, 44]. Additionally, increased variance in expression level may indicate that a gene is not critical [30], and therefore subject to less evolutionary constraint, which could theoretically manifest in an accelerated rate of protein change for genes with high expression variance. This would mean that genes which are tissue specific, have low expression levels, as well as high expression variances, are also likely to exhibit accelerated rates of protein evolution. Sex-biased genes possess all the characteristics, and it seems therefore likely that these characteristics play an important role in the rapid rate of evolution for sex-biased genes.

Genomic distribution

We found significant clustering of similarly expressed genes. Neighboring genes in both the heart and gonad showed significant correlations in sex-averaged expression level. Additionally, when assessed just by fold-change, sex-biased genes in all three tissues showed significant clustering, though only gonad showed significant non-random distribution after adjusted-p-values were used to denote significantly sex-biased genes. All the observed clustering patterns, for both expression level and for sex-biased expression may best be explained by local effects such as bidirectional promoters and transcriptional read-through, which cause the transcription of genes to influence expression of their neighbors [45].

GC3 and codon bias

Work in D. melanogaster has demonstrated different levels of codon-bias among sex-biased expression classes [14], and we recovered a similar pattern in chicken. Specifically, we observed codon-bias in female-biased genes expressed in the gonad. This could be explained be due to selection for specific codon usage patterns, or to differing GC3 levels that we observed across sex-biased categories (Fig. 4), as GC3 patterns can influence codon bias [36]. When we compared the codon-bias values of each sex bias category with those predicted by GC3 from the model fitting, we found no significant differences (one-way ANOVA on residuals), suggesting that there no evidence for any differences in selection for codon-usage between sex-biased and unbiased genes. Our observations therefore suggest that in chicken, as in mammals, codon-usage bias is mainly determined by the GC content of the isochore where each gene is situated, rather than natural selection acting on some aspect of gene expression [46].

Gene Ontology

We found some over-representation of Gene Ontology terms in our sex-biased expression categories, particularly for those sex-biased genes expressed in the gonad. It is of course tempting to speculate about the possible reasons that over-represented GO terms are associated with sex-bias expression classes. While significantly over-represented GO terms from the male-and female-biased lists lack an immediately obvious sex-specific functionality, the unbiased gonad list is striking in the number of basal metabolic terminologies. Housekeeping and basal-functioning genes would not be expected to exhibit sex-biased expression, partly because they are important to both sexes and partly because they are generally broadly expressed and therefore less likely to exhibit sex-biased expression, as described above. This may explain the over-abundance of such terms as cellular metabolic process, mitochondrion, and ligase activity, which are theoretically common to all organs and tissues.

Sample collection and preparation

White Leghorn embryos were euthanized after 18 days of incubation (ed18), and the brain, heart, and left gonad collected from four male and four female individuals. This data was originally collected for an analysis of dosage compensation, and further details about study design, tissue preparation, RNA isolation, and microarray hybridization and analysis have been previously published [17]. We note only major experimental details and where the analysis presented here differs. All microarray data is available at the Gene Expression Omnibus (project GSE 8693).

Annotation of probe sets

Annotations and genomic locations for the microarray probe sets were extracted from Ensembl [50] using version 2.1 of the chicken genome (WASHUC2 May 2006). All genes that were not assigned to a specific chromosome in this build were removed from further analysis, as many of these genes may be sex-linked and show dosage effects in expression. Additionally, as the W chromosome assembly is incomplete, the few W-assigned gene targets were also removed. The remaining dataset was comprised of coding regions that have been definitively mapped to either the autosomes or the Z chromosome. This dataset was then parsed into autosomal and Z components for all subsequent analyses, as the lack of effective and complete dosage compensation on the chicken Z chromosome [17, 28] results in male-biased expression for Z-linked genes due to simple dosage effects, which confounds subsequent studies of expression bias.

Identification of sex-biased genes

For each actively transcribed target, the log₂ (average male/average female) fold-change value was calculated. This treatment results in negative values for genes that are female-biased, positive values for genes that are male-biased, and values approaching 0 for genes with similar expression in both sexes. Genes were further divided into fold-change cutoff categories, based on the Wilcoxon test statistic of the log₂ values, and adjusted for a 5% false discovery rate using the Benjamini Hochberg correction. Therefore, absolute value of log₂ fold change > 1 corresponds to a 2× difference in expression in one sex compared to the other, absolute value of log₂ fold change > 1.5 represents roughly a 3× increase, and absolute value of log₂ fold change > 2 corresponds to a 4× difference. Relative expression was calculated across within-sex replicates for sex-biased (p _adj < 0.05, absolute fold-change >1) and unbiased gene classes for each tissue separately, and 95% confidence intervals were determined via bootstrapping (1000 repetitions).

In order to investigate the relationship between expression variance, an indicator of noise, and fold-change, we normalized the expression level for each gene target such that the averaged expression level was equal to one within each set of sex-specific replicates. This reduces the problem where genes with low expression levels will have associated small variance estimates, which could create a false signal. We then computed average within-sex variance for each sex-bias expression class (male biased log₂ fold change > 1, p _adj > 0.05; female biased genes where log₂ fold change < 1- and p _adj < 0.05; and unbiased probe sets), and computed 95% confidence intervals via bootstrapping (1000 replicates).

Clustering of genes with similar patterns of expression

For all genes that were significantly expressed in our study, we calculated the difference in average level of expression between it and its nearest neighbor for which gene expression had been measured in the same tissue. We used sex-averaged expression to correct for sex-bias effects. We also measured the difference in fold-change between neighboring gene pairs for which we had expression data, and counted the number of times gene pairs both had evidence for significant biased expression in the same sex (absolute log₂ fold change > 1, p < 0.05). We examined the significance of these values by comparing them with datasets where gene order was randomized on each chromosome using 10000 replicates.

GC3 and bias in codon usage

Third codon GC (GC3) content has been linked to codon-bias [36, 47], which has in turn been linked to sex-biased gene expression patterns in Drosophila [14]. Therefore, in order to investigate the relationship between sex-biased gene expression and codon bias, we calculated GC3 and the Effective Number of Codons (ENC) [37] from the Ensembl annotation of the chicken genome. ENC can range from 20, in the case of extreme bias where only one codon is used for each amino acid, to 61, where all alternative codons are equally likely. We calculated the average GC3 and ENC for each sex-bias expression category within each tissue at the fold-change < 1- or > 1, p _adj < 0.05 cutoff. 95% confidence intervals were computed with bootstrapping (1000 replicates). Within each tissue, we compared the GC3 and ENC among autosomal expression categories using one-way ANOVA (2 d.f.), and between Z-linked categories using a one-tailed t-test assuming unequal variance. As ENC is strongly governed by GC content in vertebrates, we also tested whether differences in ENC values between categories were significantly different from that predicted by the observed differences in GC3. To do this we fitted a smoothed cubic spline to the relationship between ENC and GC3, and used this to estimate predicted ENC values from GC3. Significance was tested performing a one-way ANOVA on the residuals of model fitting.

Gene Ontology analysis

We compared the functional role of biased and unbiased gene sets for each tissue using Ontologizer [48] and the Gallus Gene Ontology (GO) [49] database to identify over-represented GO categories. For each tissue, the actively expressed genes were used as the population gene set in term-for-term comparisons with each of the sex-biased expression categories as separate study sets (absolute fold change > 1, p _adj < 0.05 cutoffs to denote sex-bias). Significant over-representation was determined for each term with Fisher's exact test, and the Bonferroni correction was used for multiple comparisons. As the Bonferroni method can be unnecessarily over-conservative, we report over-represented terms that are significant before Bonferroni correction as well.