Alternative reproductive tactics and sex-biased gene expression: the study of the bulb mite transcriptome

The sexes experience different selective pressures, which can lead to highly divergent phenotypes that are achieved via sex-biased gene expression. The effect of sexual dimorphism on the degree of sex-bias in gene expression can be studied in species characterized by sexually selected alternative male phenotypes. We analyzed gene expression in the bulb mite Rhizoglyphus robini (Acari, Acari-dae) , in which more sexually dimorphic, aggressive ﬁghter males, possessing thickened legs of the third pair which are used to kill rivals, coexist with unar-mored scrambler males. We sequenced transcriptomes of adult females and both types of males and de-novo assembled 114,456 transcriptome-based gene models (TGMs). Signiﬁcantly more TGMs had male-biased expression than female-biased expression. Among TGMs that were over expressed in one, but not both, male morphs (compared to expression in females), we found about four times more ﬁghter-biased genes than scrambler-biased genes. This demon-strates that the degree of expression bias reﬂects the degree of sexually selected dimorphism. However, the number of sex-biased genes was much higher than the number of genes differentially expressed between male morphs, and most male-biased genes were shared between morphs, suggesting that selection pressures act similarly on males irrespective of their morph. Furthermore, we found that male-biased genes evolved at a faster rate than female-biased genes, as evidenced by a higher rate of both gene-turnover and amino acid substitution, indicating that sexual selection, acting more strongly on males, accelerates the rate of molecular evolution. Interestingly, gene turnover was relatively higher, but amino acid substitution rate relatively lower among ﬁghter-biased genes, suggesting that different components of sexual selection may have different effects on the evolution of sex-biased genes.


Introduction
The sexes, which are defined by the types of gametes they produce, experience different selective pressures that often result in highly divergent phenotypes. While the sexes may differ with respect to the loci associated with sex chromosomes, divergent phenotypes are mostly achieved by sex-bias in gene expression patterns, defined as the predominant or exclusive expression of genes in one sex (reveiwed in Ellegren and Parsch 2007;Parsch and Ellegren 2013).
Sexual dimorphism is thought to be driven to a large extent by sexual selection, a process arising from reproductive competition, and leading to the evolution of traits that facilitate access to partners of the opposite sex and their gametes. This process, acting mostly on males, has led to the emergence of a vast diversity of sexually dimorphic traits, such as body size, weapons, sexual ornaments or ejaculate traits (reviewed in Andersson 1994). However, the contribution of sexually selected traits to the evolution of differences in gene expression patterns between sexes is not well documented.
Reproductive competition sometimes leads to the emergence of alternative reproductive phenotypes in males, where large, often armored, territorial males coexist with less sexually dimorphic phenotypes adopting nonaggressive ways to secure access to females (reviewed in Gross 1996;Brockmann 2001;Oliveira et al. 2008).
Alternative mating tactics thus offer the potential to investigate the role of the sexual selection-driven dimorphism in the evolution of sex-biased gene expression, but this potential has started to be realized only recently (Snell-Rood et al. 2011;Pointer et al. 2013). In a recent study that specifically addressed this issue Pointer et al. (2013) found that the less sexually dimorphic male morph is less masculinized, and more feminized compared to the dominant morph in a wild turkey.
Apart from affecting sexual dimorphism, sexual competition is also likely to affect the rate of evolution of genes with sex-biased expression. Sex-biased genes have been shown to evolve at a faster rate than other genes, and this pattern has been observed both as an increased rate of sequence evolution and an elevated gene turnover rate (e.g., Ranz et al. 2003;Zhang et al. , 2007Haerty et al. 2007;Assis et al. 2012). For example, Drosophila melanogaster males gain an advantage in sperm competition, defined as competition between ejaculates coming from different males that have inseminated the same female (Parker 1970), through the use of male accessory gland proteins (Clark et al. 1995). These proteins show a high dN/dS ratio (e.g., Begun et al. 2000;Swanson et al. 2001) and have a variety of effects on the physiology of the inseminated female, including increasing the female's egg-laying rate (Chen et al. 1988;Heifetz et al. 2000). These effects benefit males, but they may harm females; for example, they may reduce female lifespan (Chapman et al. 1995). Such antagonistic effects are a manifestation of interlocus sexual conflict, whereby traits beneficial for members of one sex harm their partners of the opposite sex (Parker 1979). Such processes can lead to an evolutionary arms race, in which male adaptations are matched by female counter-adaptations that neutralize harmful effects, thus accelerating the rate of evolution (Chapman et al. 2003;Arnqvist and Rowe 2005;Ellegren and Parsch 2007;Parsch and Ellegren 2013;Pennell and Morrow 2013).
Here, we sequence and assemble, for the first time, the transcriptome of the male-dimorphic bulb mite, Rhizoglyphus robini (Astigmata: Acaridae), and identify genes with sex-biased expression. The bulb mite is a model species in sexual selection and conflict research (e.g., Radwan 2004Radwan , 2007. As in some other species of acarid mites (reviewed in Radwan 2009), males of this species are dimorphic. While the morphs do not differ in size (Radwan 1995), the third pair of legs in fighter males is greatly thickened (Fig. 1) and used in intraspecific contests as a lethal weapon, whereas in benign scrambler males and in females, these legs are unmodified (Radwan et al. 2000). Male morph is heritable in R. robini (Radwan 1995). Postcopulatory sexual selection also plays an important role in this species, as high female promiscuity induces intense sperm competition among males (Radwan and Siva-Jothy 1996). Finally, experimental evolution under varying levels of sexual selection has demonstrated that male adaptations to sperm competition are harmful to females (Tilszer et al. 2006), providing evidence for interlocus sexual conflict in this species. Thus, both pre-and postcopulatory selection as well as sexual conflict could be expected to increase the rate of evolution of sex-biased genes in the bulb mite. Because some traits associated with strong sexual selection (intrasexual aggression and weapons) occur only in the fighter morph, we expected to find a larger number of genes with higher expression compared to expression in females in the fighter morph and a smaller number of such genes in the scrambler morph. Furthermore, we compared the rate of evolution of male-biased genes with that of female-biased genes, and investigated whether the rate of evolution of malebiased genes is affected by male morph.

Samples and sequencing
The mites used for this study came from a population collected from onions in a garden near Krakow, Poland and maintained on yeast at 24°C at >90% humidity. Before collection, a cohort of mites (AE3 days) were maintained under the common garden conditions. Nymphs were collected in a nonfeeding, quiescent stage, and adults were maintained on filter paper for 24 h before collection in order to clear their digestive systems of yeast food. Total RNA was extracted from adult females, two male morphs, and tritonymphs. About 100 individuals in each group were used for the extraction, therefore expression levels we report can be treated as means over many individuals within each sex and morph.
RNA was extracted using RNAzol-RT; residual DNA was removed by DNAse treatment. Sequencing libraries were prepared from polyA+ RNA with the TrueSeq RNA kit and 2 9 100 bp paired-end (PE) sequencing was performed in a single Illumina HiSeq2000 lane.

Assembly and transcriptome-based gene model
The lack of an annotated reference genome for the bulb mite or any closely related species hampered the analysis of differential expression. Without a reference genome, distinguishing among alleles, close paralogs, and alternatively spliced transcripts becomes difficult or impossible. De-novo transcriptome assemblers reconstruct a large fraction of expressed transcripts, including those derived from recently duplicated genes and alternatively spliced isoforms, although simulations show that many isoforms can be erroneously inferred (Vijay et al. 2013). To avoid problems with redundancy in the assembled transcriptome, we used transcriptome-based gene models (TGMs) as the basis of the differential expression (DE) analysis.
In order to assemble the transcriptome, we first removed low-quality bases with the DynamicTrim program from the SolexaQA package (Cox et al. 2010). Transcriptome was de-novo assembled using all 154.6 million PE reads; 26.3, 26.2 and 20.2 million PE reads were from fighters, scramblers, and females, respectively and 81.8 million of PE reads from tritonymphs; tritonymphs were used only for transcriptome assembly and not for expression analysis. To reconstruct transcripts, we used Trinity de-novo transcriptome assembler (release 2012-06-08; (Grabherr et al. 2011)). Trinity was run with default parameters, except -bfly_opts "-edge-thr = 0.25", which reduces the sensitivity of alternative splice detection. After transcriptome assembly, we created a bioinformatics pipeline ( Fig. S1), consisting of four steps, to build the TGMs.
First, the Trinity assembly was used as input for CAP3 assembler (Huang and Madan 1999) with stringent parameters (-o 40 -p 99); the aim of this step was to obtain an extension of almost perfectly overlapping sequences. Second, CD-HIT (Li and Godzik 2006) was used to cluster sequences with at least 95% identity. The main purpose of this step was to collapse alleles from the same locus into a single, representative sequence. We expected that in the CD-HIT-reduced data, alternatively spliced sequences from the same locus would be retained as separate contigs that would be characterized by high sequence identity over most of their length but possibly contain inserted or deleted sequences corresponding to alternatively spliced exons. To merge these isoforms into a single sequence containing all expressed exons, we applied the following procedure (third step): MegaBLAST (Altschul et al. 1990) was used to group sequences which represented potential isoforms of the same gene. We applied a single-linkage clustering algorithm, where each sequence within a group fulfilled the following similarity criteria with at least one other sequence within the group: (1) the total length of all high-scoring pairs (HSP) for shorter sequences covered at least 70% of the total sequence length; and (2) the alignment identity between sequences was at least 96%. In the last step, we aligned all sequences within a group using the multiple global alignment algorithm implemented in MAFFT (Katoh et al. 2002), which allows the insertion of large (e.g., exon-size) gaps into the sequence alignment. If the overall identity was at least 96%, a single consensus sequence was generated for the group. Otherwise, alignments were divided into smaller groups that fulfilled the identity criteria and separate consensus sequences were generated. All consensus and singleton sequences were included in the TGM database.

Identification of open reading frames (ORFs)
We used the pipeline implemented in Trinity to identify and predict likely coding sequences in the TGM database (open reading frames, ORFs). With this approach, the longest ORFs (at least 100 amino acids long) in each of the six possible reading frames were scored according to a Markov Model based on hexamers. If the score of the putative proper ORF coding frame was positive and was higher than that of the presumed incorrect reading frames, then the ORFs in that frame were reported (Grabherr et al. 2011). Only complete ORFs were retained.

Gene expression
Gene expression analysis was performed with the pipeline included in Trinity. Briefly, for each sample, reads were mapped to TGMs with Bowtie (Langmead et al. 2009), then transcript abundance was estimated by RSEM (Li and Dewey 2011). Next, we identified and analyzed differentially expressed genes and performed all pairwise comparisons between females, fighters, and scramblers, using edgeR (Robinson et al. 2010), a Bioconductor package (Gentleman et al. 2004), with trimmed mean of M-values normalization, which accounts for differences between samples in the number of reads. We used a relatively restrictive 0.01 false discovery rate (FDR) to identify sexbiased genes, in order to avoid diluting trends in downstream analyses with false positives.

Sequence divergence and expression
Tetranychus urticae, which diverged from R. robini about 400 mya, is the most closely related mite species with a sequenced genome. Such an ancient divergence makes an analysis of dN/dS unfeasible, because of saturation issues. Instead, we adopted the strategy of (Snell-Rood et al. 2011) and used the amino acid distance (dA) between T. urticae and R. robini orthologous sequences as our measure of the rate of evolution of protein-coding sequences. To find one-to-one R. robini-T. urticae orthologs, we used BLASTX (e-value ≤ 10 À8 ) to find R. robini ORFs with a single hit in the database of T. urticae protein sequences (ver. 20120618). For each pair of orthologs we used estwise (Birney et al. 2004) to provide a proteinguided translation of R. robini nucleotide sequences. Predicted protein sequences were extracted from the estwise output and aligned in MAFFT. Amino acid (AA) distances were calculated in PAML4 (Yang 2007) as the maximum likelihood estimate of the number of amino acid replacements per site under the empirical substitution model WAG (Whelan and Goldman).

Functional annotation, gene ontology, and enrichment analysis
Gene ontology analyses were performed in Blast2GO (ver. 2.6.0; (Conesa et al. 2005)), a tool for the assignment of functional categories to genes or gene products. We performed a BLASTX search in the SwissProt database (Bairoch et al. 2009) using our transcriptome and imported the results into Blast2GO, where the functional annotations and enrichment analysis steps were done with default settings and an FDR of 0.01.

Results
A single sequencing run of a RNA pool resulted in 154.6 million of 2 9 100 bp paired end reads (Table 1). Raw reads were deposited in BioProject portal (PRJNA213807). We obtained 137,151 transcripts in the original Trinity assembly. These were reduced to 114,456 TGMs (including protein coding and noncoding expressed sequences) after applying the pipeline for merging transcripts likely derived from the same genomic location. 96,826 TGMs were found in fighter males, 100,157 in scrambler males, and 84,600 in females. Majority of our protein-coding TGMs appear to represent nearly complete coding sequences. More than a half (56.0%, 3877 of 6911) of proteins from SwissProt database (Bairoch et al. 2009) showing significant similarity to our TGMs, were covered at >70% of their length by the best-hit TGMs (Fig. S2), and most proteins (78%) were covered in more than 50% of their length.
At the 0.01 FDR, the expression level of 4.0% of TGMs was higher in both male morphs compared to females (male-biased genes), and 1.3% of TGMs had higher expression in females (female-biased genes; Table 2., Fig. 2). The difference between the sexes in the proportion of biased genes was highly significant (v 2 = 1521.4, P < 0.001). In pairwise comparisons between females and each male morph, we found that there were more genes with higher expression in fighters (fighter-only-biased: 2.4% of all TGMs vs. scramblers-only-biased: 0.5% of all TGMs; v 2 = 1595.0, P < 0.001; Fig. 3). Similarly, we found a larger percentage of genes with higher expression in females in comparison with fighters than in females in comparison with scramblers (0.6% vs. 0.1% of all TGMs; v 2 = 470.3, P < 0.001; Table 2, Fig. 3). The unbiased genes, i.e., those with no significant differences (at 0.05 FDR level) in expression between females and either male morph, constituted 83% of all TGMs. In the comparison between male morphs, 297 genes (0.26% of TGMs) showed higher expression levels in fighters, and 76 genes (0.07%) showed higher expression levels in scramblers.
Male-biased genes were significantly overrepresented in many gene ontology (GO) categories that corresponded to catalytic activity and metabolic processes (see Fig. S3, Fig. S4). GO categories associated with catalytic activity also dominated among fighter-only-biased genes (Fig. S5); however, we also found overrepresentation of genes associated with gene expression, RNA, nucleic acid, and chromatin binding, as well as several subcategories within the cellular component (Fig. S5). Among scrambler-onlybiased genes, only peptidase activity and lipid metabolic process categories were overrepresented. In contrast, female-biased genes were mostly overrepresented in the cellular component GO category e.g., microtubule organizing center (Fig. S6). The numbers of TGMs with complete ORFs for which orthologs were identified in the T. urticae genome are shown in Table 3. A significantly lower proportion of orthologs was found among male-biased ORFs than among unbiased ORFs (v 2 = 126.72, P < 0.001). This was also the case for the proportion of orthologs among femalebiased ORFs (Table 3; v 2 = 12.71, P = 0.004), although the difference compared to unbiased ORFs was less pronounced than in the case of male-biased ORFs (Table 3). Indeed, the proportion of hits to T. urticae was significantly lower among male-biased ORFs compared to female-biased ORFs (v 2 = 10.57, P = 0.0012). Interestingly, among ORFs with higher expression levels (compared to females) in only one of the male morphs, we found 2.5 times as many T. urticae orthologs for scrambler-only-biased ORFs than for fighter-only-biased ORFs (v 2 = 32.10, P < 0.001; Fig. 3).
Means and standard deviations of amino acid distances between T. urticae and R. robini orthologs are shown in Table 4. Male-biased genes evolved at a significantly faster rate than unbiased genes (t 1101.1 = 15.83, P < 0.001) and female-biased genes (t 499.8 = 6.25, P < 0.001). We found the same result when we analyzed scrambler-only-biased genes (t 131.3 = 5.39, P < 0.001 and t 216.6 = 3.96, P < 0.001 for comparisons with unbiased and female-biased genes, respectively). In contrast, we found no significant differences when we analyzed fighter-only-biased genes (t 266.1 = 1.81, P = 0.072 and t 496.3 = 0.59, P = 0.7226 for comparisons with unbiased and female-biased genes, respectively). Finally, there was no significant difference in the rate of evolution between female-biased and unbiased genes (t 342.6 = 1.48, P = 0.139).

Discussion
Among the total of 114,456 Transcriptome-based Gene Models (TGMs), 16,630 contained complete ORF. The number of ORFs which are likely to correspond to protein coding genes is comparable to that of other Metazoans (EnsemblMetazoa ver.21; Kersey et al. 2012). It is difficult to tell whether the number of TGMs not containing ORF is unusually high, because good comparative data are not available. However in general a large fraction of eukaryotic genome is expressed and thus we expect a large number of noncoding RNAs (Mercer et al. 2009;Ulitsky and Bartel 2013). Because their expression is often low and/or tissue specific, their number may be inflated due to incomplete assembly.

Morph-specificity of male-biased genes
While we did not perform replicated sequencing of males and females for logistic reasons, each sample consisted of about 100 individuals so that random differences in expression patterns between individuals could be expected to cancel out. Furthermore, about 90% of scramblerbiased genes were also found to be fighter-biased (Fig. 3), indicating that our identification of sex-biased genes was reliable.
Transcriptome-based gene models for which both male morphs had a higher level of expression than females, were three times as numerous as female-biased TGMs (Table 2). This pattern has also been noted in other systems (Malone et al. 2006;Zhang et al. 2007;Small et al. 2009) and is consistent with the hypothesis of maledriven diversification in gene expression levels between the sexes. These genes constituted 57% of all genes that showed higher expression than female genes in at least one male morph. The prediction that more sexually dimorphic fighter morph will be characterized by a larger number of genes with expression level higher compared  to females was confirmed by our results: fighter-onlybiased genes were about four times as numerous as scrambler-biased genes. Fighter-biased genes constituted 35% of the genes that had higher expression levels in at least one male morph. Thus, morphological (thickened legs) and behavioral adaptations to intrasexual combat characterizing the fighter morph have apparently driven differential expression of a considerable proportion of the genome. In male-dimorphic dung beetle Onthophagus thaurus, hornless males were more similar to females, compared to horned males, in gene expression patterns in developing head and brain, but not in leg (Snell-Rood et al. 2011). In the wild turkey (Meleagris gallopavo), less sexually dimorphic subordinate males were more feminized in gonads, but interestingly, also in spleens, compared to dominant males (Pointer et al. 2013). Similar to the result of Pointer et al. (2013), we found that most of the male-biased genes were shared between morphs (Fig. 3). This indicates that even though sexually selected morphology of fighters is associated with specific patterns of gene expression, most of the physiological pathways are shared between morphs.
The number of TGMs whose expression level differed significantly between male morphs was about 15 times lower than those differentially expressed between the sexes (0.26% vs. 4%). This finding contrasts with the result obtained by Snell-Rood et al. (Snell-Rood et al. 2011) for the dung beetle, in which morph-bias was of a similar order of magnitude as sex-bias. Snell-Rood and colleagues interpreted their findings as a support for "developmental decoupling hypothesis" (West-Eberhard 1989) posing that morph or environment-specific gene expression patterns allow independent evolution of divergent phenotypes within species. However, they only analyzed selected tissues (horn, brain, and leg) in the developing pupae, and they did not include genital tissues known to be highly sexually dimorphic in expression in other systems (Parisi et al. 2004;Small et al. 2009;Assis et al. 2012). On the other hand, some genes that are differentially expressed between morphs are probably expressed only during certain developmental stages, and these could not be detected in our study. This is particularly likely during the last nymphal stage, when genes determining sexual dimorphisms and alternative male morphologies are expressed. Morph-bias in gene expression during the last nymphal instar in R. robini awaits further investigation.
Gene orthology analysis indicated that male-biased genes were overrepresented in categories corresponding to metabolic processes and catalytic activities associated with energy metabolism, such as kinase activity. This finding is consistent with the higher motility of males, which spend a fair amount of time searching and competing for females. Selection for efficient energy metabolism is likely stronger in fighter males, as fights may be very long and require fast mobilization of considerable resources. Consistent with this hypothesis, genes in many metabolism-related categories were fighter-only biased. Thus, evolution of sexually selected male morphologies and behaviors appears to be associated with changes in sexspecificity in expression of a large number of genes involved in metabolic processes.

Rate of evolution of sex-biased genes
Male-biased genes, and to a lesser degree female-biased genes, evolved faster than the genomic average, as evidenced by a significantly lower fraction of orthologs identified among the spider mite proteins. Furthermore, male-biased genes evolved faster than female-biased genes. Faster evolution of male-biased genes has also been reported for D. melanogaster and Drosophila pseudoobscura, which diverged 25-50 million years ago (Zhang et al. 2007;Assis et al. 2012). Fighter-only-biased ORFs represent an especially low proportion of the identifiable orthologs in T. urticae. This pattern may, in part, reflect the fact that fighter morphology and behavior, associated with a number of male-biased genes, do not exist in T. urticae. In contrast, we found that the proportion of female-biased ORFs with identifiable orthologs in T. urticae was more similar to that for unbiased genes.
Among ORFs with orthologs in T. urticae, the amino acid substitution rates were higher among male-biased genes than among female-biased genes, but the difference was only 12.5%. This result may reflect the fact that R. robini and T. urticae diverged about 400 mya (Dabert et al. 2010). It therefore seems likely that the orthology we identified was predominantly found in highly conserved genes, whereas for most faster-evolving genes, orthology could no longer be inferred. Thus, if one compares the rate of evolution of male-biased and femalebiased genes using distantly related taxa, a stronger signal could be expected to come from the gene-turnover rate than the rate of molecular evolution. Interestingly, significantly fewer orthologs to T. urticae were found among fighter-only-biased ORFs than among scrambler-onlybiased ORFs (Table 3), but when we analyzed the rate of evolution, the situation was reversed, with scrambleronly-biased ORFs evolving faster, and at a similar rate to male-biased ORFs in general (Table 4). This result may suggest that the evolution of genes associated with fighter morphology and behavior occurred either long ago or at a very high rate, such that orthology to any ancestral genes was blurred. On the other hand, male-biased and scrambler-only-biased genes may keep evolving because of the constant arms race associated with sperm competition and sexual conflict. The slower evolution of fighter- only-biased genes is not consistent with the idea that morph-specific gene expression should result in greater sequence divergence due to relaxed selection (Snell-Rood et al. 2010;Van Dyken and Wade 2010).
Alternatively, faster rate of evolution of male-biased genes may be due to tissue specificity of their expression. Sex-specific genes may often be characterized by tissuespecific expression patterns (e.g., reproductive proteins (Swanson and Vacquier 2002); brain transcripts (Mank et al. 2007)). Because genes with tissue-specific expression are likely to have fewer pleiotropic effects than genes with widespread expression, they may be less evolutionarily constrained and, indeed, there is abundant evidence that they evolve at a faster rate (Duret and Mouchiroud 2000;Zhang and Li 2004;Haerty et al. 2007;Ekblom et al. 2010). The fast evolution of sex-biased genes might thus, in principle, be due to their tissue-specificity. We are unable to resolve this in the bulb mite, but Meisel (2011) has recently shown in D. melanogaster that even among narrowly expressed genes, those that are sex-biased evolve the fastest.

Conclusions
Comparison of male and female transcriptomes in a mite species expressing alternative mating tactics provided support for the hypothesis that dimorphism driven by sexual selection is responsible for sex-bias in expression of a significant part of the bulb mite genome. Our data are consistent with those reported for other species in showing that male-biased genes evolved at a faster rate than female-biased genes. This faster evolution was reflected both in a higher gene-turnover rate (as revealed by the fraction of genes with T. urticae orthologs) and amino acid substitution rate (measured as the AA divergence from T. urticae orthologs). Interestingly, fighter-onlybiased genes showed patterns of evolution different from that of male-biased and scrambler-only-biased genes. A lower fraction of fighter-only-biased genes had orthologs in T. urticae, but these orthologs were characterized by lower amino acid divergence. This suggest that genes associated with the distinct morphology, behavior and physiology of fighter males evolve in different way compared to genes which are involved in processes common to both morphs, such as sperm competition, possibly because the latter are involved in inter-locus conflict, resulting in continual evolution of male-biased genes. Figure S1. Pipeline for the construction of transcriptomebased gene models. Figure S2. Distribution of SwissProt protein coverage (%) by best-hit transcriptome-based gene models. Figure S3. Molecular function domain in gene ontology for significantly overrepresented male-biased genes. Figure S4. Biological process domain in gene ontology for significantly overrepresented male-biased genes. Figure S5. Molecular function domain in gene ontology for genes fighter-only-biased genes exclusively in fighter males. Figure S6. Cellular component domain in gene ontology for significantly overrepresented female-biased genes.