Slow evolution of sex‐biased genes in the reproductive tissue of the dioecious plant Salix viminalis

Abstract The relative rate of evolution for sex‐biased genes has often been used as a measure of the strength of sex‐specific selection. In contrast to studies in a wide variety of animals, far less is known about the molecular evolution of sex‐biased genes in plants, particularly in dioecious angiosperms. Here, we investigate the gene expression patterns and evolution of sex‐biased genes in the dioecious plant Salix viminalis. We observe lower rates of sequence evolution for male‐biased genes expressed in the reproductive tissue compared to unbiased and female‐biased genes. These results could be partially explained by the lower codon usage bias for male‐biased genes leading to elevated rates of synonymous substitutions compared to unbiased genes. However, the stronger haploid selection in the reproductive tissue of plants, together with pollen competition, would also lead to higher levels of purifying selection acting to remove deleterious variation. Future work should focus on the differential evolution of haploid‐ and diploid‐specific genes to understand the selective dynamics acting on these loci.

Although sexual dimorphism is generally more extensive in animal species, male and female dioecious flowering plants also undergo conflicts over trait optima and are subject to natural and sexual selection leading to a range of phenotypic sexual differences (Barrett & Hough, 2013). Studies of differential male and female gene expression patterns in plants (Muyle, Shearn, & Marais, 2017) indicate that sex-biased gene expression plays a role in the evolution of sexual dimorphism in morphological (e.g., anther and ovule development pathways in asparagus, Harkess et al., 2015), physiological (e.g., salinity tolerance in poplars, Jiang et al., 2012) and ecological traits (e.g., response to fungal infection in Silene latifolia, Zemp, Tavares, & Widmer, 2015).
Extensive studies in plants and animals have shown that genes with sex-biased expression vary in abundance across different developmental stages and tissues (Grath & Parsch, 2016;Perry, Harrison, & Mank, 2014;Robinson et al., 2014;Zemp et al., 2016;Zluvova, Zak, Janousek, & Vyskot, 2010). Evolutionary dynamics analyses also indicate that different evolutionary pressures impact the rate of sequence evolution of sex-biased genes; for example, sex-biased genes in reproductive tissues tend to have different rates of protein evolution compared to unbiased genes (Dean et al., 2017;Lipinska et al., 2015;Mank, Nam, Brunstr€ om, & Ellegren, 2010a;Perry et al., 2014;Sharma et al., 2014). In animal systems, where the rates of sequence divergence of sex-biased genes have been studied more widely, male-biased genes in many species, including Drosophila and adult birds, tend to be more numerous and to have higher expression and divergence rates (Assis, Zhou, & Bachtrog, 2012;Grath & Parsch, 2016;Harrison et al., 2015;Khaitovich et al., 2005) compared to female-biased and unbiased genes.
This has often been interpreted as the signature of sexual selection, particularly sperm competition (Ellegren & Parsch, 2007). However, studies in other organisms have reported elevated rates of evolution in female-biased genes (Mank et al., 2010a;Whittle & Johannesson, 2013), leading to questions about the relationship between rates of evolution and sexual selection. In Arabidopsis, genes expressed in pollen have lower rates of evolution (Gossmann et al., 2014). Moreover, nonadaptive evolutionary processes have been shown to drive the fast rates of sequence evolution observed in sex-biased genes in some systems (Gershoni & Pietrokovski, 2014;Harrison et al., 2015) perhaps related to relaxed purifying selection (Hunt et al., 2011).
Sexual selection in flowering plants is also thought to be strong (Moore & Pannell, 2011), acting on gene expression patterns predominantly through pollen competition. Male gametophytic tissue in Arabidopsis thaliana and rice has been shown to express a higher proportion of recently evolved genes compared to other tissues (Cui et al., 2015). Some of these young genes possess essential pollen-specific functions, suggesting a role for pollen competition in facilitating de novo gene development.
As male-biased mutation is thought to be strong due to the elevated numbers of germ cell divisions in male cells (Whittle & Johnston, 2003), pollen competition, in this case, was suggested to counteract the potentially negative effects of higher mutation rates present in male gametophytes (Cui et al., 2015). Similarly, younger genes in the gametophyte of A. thaliana, rice and soya bean were also found to have higher rates of evolution compared to genes in the sporophytic tissue, however to varying degrees in males and females (Gossmann, Saleh, Schmid, Spence, & Schmid, 2016). Suggested reasons for these findings concerned the lower tissue complexity, and hence lower genetic interaction, in the gametophyte as well as differences between the sexes.
Plants additionally differ from animals in having a longer haploid phase in their life cycle, suggesting that haploid selection may act more forcefully to remove mildly deleterious recessive variation in pollen-expressed genes. Previous work on A. thaliana showed that the predominance of selfing, and similarly the intragametophytic selfing in moss species (Sz€ ov enyi et al., 2014), leads to the more effective purging of mildly deleterious recessive variation (Gossmann et al., 2014). In the obligate outcrossing plant Capsella grandiflora, pollen-specific genes, but not spermenriched genes, evolve under both stronger purifying and positive selection compared to genes from sporophytic tissues (Arunkumar, Josephs, Williamson, & Wright, 2013). These findings are indicative of a potential combined effect of haploid selection and pollen competition acting in pollen-specific cells, whereas selective pressures are expected to be more relaxed for sperm-specific genes as there is no competition between them (Arunkumar et al., 2013).
These studies make it increasingly clear that many evolutionary forces shape the sequence evolution of sex-biased genes, including sexual selection through sperm competition (Ellegren & Parsch, 2007), haploid selection and natural selection (Ingvarsson, 2010).
Particularly in plants, in order to understand the relative contribution of these forces, it is important to study rates of evolution in species with different levels of gamete competition, motivating studies on outcrossing dioecious species.
Willow and poplar species have reproductive structures characterized by clusters of unisexual inflorescences referred to as catkins ( Figure 1). Flowers in male willow catkins present a reduced number of stamens with anthers and filaments; however, they lack a vestigial DAROLTI ET AL.
Flowers in female willow catkins contain pistils with style, stigma and an ovary. However, they also show a high degree of floral reduction as there is an absence of staminodes and, similarly to male catkins, they lack a perianth with petals and sepals (Cronk et al., 2015;Fisher, 1928;Karrenberg, Kollmann, & Edwards, 2002), potentially with a role in facilitating wind pollination (Karrenberg et al., 2002).
Our study of gene expression patterns in male and female S.
viminalis individuals begins to explore the selective forces acting on sex-biased gene evolution in dioecious plants. We analysed sex-biased gene expression patterns in S. viminalis from two different tissues, vegetative (leaf) and sex-specific reproductive (catkin) tissue. We found the reproductive tissue to be more transcriptionally dimorphic and identified overall higher expression levels for male-biased genes than for female-biased genes, consistent with previous studies (Grath & Parsch, 2016). Interestingly, however, we found that in catkin, male-biased genes on the autosomes and the pseudoautosomal region have significantly lower rates of sequence divergence than both unbiased and femalebiased genes. Similarly, female-biased genes show lower rates of sequence evolution in comparison with unbiased genes; however, the difference is not significant. We could not detect any significant differences in the proportion of genes evolving under positive selection between either male-biased or female-biased genes and unbiased genes. The low rates of male-biased sequence evolution could be partly explained by the higher rate of silent mutations in male-biased genes resulting from lower codon usage bias.
However, haploid selection would also be expected in this tissue to exert a stronger purifying force to remove deleterious recessive mutations.

| Sample collection and sequencing
We obtained RNA-seq data from leaves and catkins from three  (Lohse et al., 2012) to remove adaptor sequences and trim the reads, removing regions where the average Phred score in sliding windows of four bases was <15 as well as reads for which the leading/trailing bases had a Phred score <3. Following trimming, we removed paired-end reads where either read pair was <50 bp (Table S1), resulting in an average of 30 million paired-end reads per sample.

| Expression analysis
We mapped RNA-seq reads to the de novo male genome assembly (Pucholt et al., 2017)  using BLASTN and an e-value cut-off of 1 9 10 À10 .
We extracted read alignments for transcripts in each sample and tissue separately from the filtered transcriptome reference using STRINGTIE and obtained read counts using HTSEQ v.0.6.1 (Anders, Pyl, & Huber, 2015). RPKM values were estimated using EDGER (Robinson, McCarthy, & Smyth, 2010) in R (R core team 2015) and transcripts filtered for a minimum expression threshold of 2 RPKM in at least half of the individuals in one sex (in this case, at least two of the three individuals per each sex) as per previous similar studies Pointer et al., 2013). We only retained transcripts with positional information on annotated chromosomes (Pucholt et al., 2017) for further analysis and normalized separately for each tissue using TMM in EDGER.
We performed hierarchical clustering of average gene expression for genes expressed in both tissues with bootstrap resampling (1,000 replicates) in the R package PVCLUST v.2.0 (R Core Team, 2015;Suzuki & Shimodaira, 2006). We generated a heatmap of log 2 average male and female expression in the two tissues using the R package PHEAT-MAP v.1.0.7 (Kolde, 2012;R Core Team, 2015). We identified sex-biased expression based on a minimum of twofold differential expression (log 2 M:F RPKM > 1 for male-biased expression and < À1 for female-biased expression) and a significant p value (p adj < .05 following FDR correction for multiple testing (Benjamini & Hochberg, 1995)) in EDGER.
The longest transcript for each gene was identified in all species, and a reciprocal BLASTN with an e-value cut-off of 1 9 10 À10 and a minimum percentage identity of 30% was used to identify orthologs.
To ensure the accurate calculation of divergence estimates, poorly aligned regions were masked with SWAMP v 31-03-14 (Harrison, Jordan, & Montgomery, 2014). We employed a two-step masking approach, first using a shorter window size to exclude sequencing errors causing short stretches of nonsynonymous substitutions and then a large window size to remove alignment errors caused by variation in exon splicing . Specifically, we first masked regions with more than seven nonsynonymous substitutions in a sliding window scan of 15 codons, followed by a second masking where more than two nonsynonymous substitutions were present in a sliding window scan of four codons. To choose these thresholds, we imposed a range of masking criteria on our data set and conducted the branch-site test on these test data sets. We manually observed the alignment of genes with the highest log likelihood scores to choose the most efficient and appropriate masking criteria. We subsequently removed genes where the alignment (after removal of gaps and masked regions) was < 300 bp, which likely represent incomplete sequences. This resulted in 7,631 1:1 orthologs.
We tested the robustness of the 1:1 orthologs data set (Supporting Information) by separately inferring orthologous groups using ORTHOMCL (Li, Stoeckert, & Roos, 2003), an approach with higher specificity (Altenhoff & Dessimoz, 2009). As ORTHOMCL relies on the Markov Clustering algorithm, it is useful in identifying cases of coorthology (a duplicate of a gene in one species that is orthologous to a gene in another species) within the total 1:1 orthologous groups identified. By excluding these co-orthologous groups, we recovered fewer 1:1 orthologs (1,346 after filtering for polymorphism and divergence data); however, the divergence results were consistent with our broader data set based on reciprocal BLAST (Table S2). As such, we concluded that the reciprocal best-hit approach was appropriate to use in this case.
We further used branch model 2 (model = 2, nssites = 0, fixomega = 0, omega = 0.4) from the CODEML package in PAML v4.8 (Yang, 2007)  | 697 Based on their genomic location in the S. viminalis genome (Pucholt et al., 2017), we divided orthologs into two groups, orthologs on the autosomes (including the pseudoautosomal region of the Z chromosome) and orthologs on the Z-linked nonrecombining region. Because genes on sex chromosomes can exhibit accelerated rates of evolution (Charlesworth, Coyne, & Barton, 1987), and this may be more often due to nonadaptive processes on Z chromosomes (Mank, Vicoso, Berlin, & Charlesworth, 2010b;Wright et al., 2015), we analysed rates of evolution separately for autosomal and

| Polymorphism analysis
We obtained polymorphism data by mapping the RNA-seq reads to the reference genome assembly using STAR aligner v2.5.2b (Dobin et al., 2013) in the two-pass mode and with default parameters, retaining uniquely mapping reads only. We conducted SNP calling using SAMTOOLS mpileup and VARSCAN v2.3.9 mpileup2snp (Koboldt et al., 2012). We ran SAMTOOLS mpileup with a maximum read depth of 10,000,000 and minimum base quality of 20 for consistency with VARSCAN minimum coverage filtering. The base alignment quality (BAQ) adjustment was disabled in SAMTOOLS as it imposes a too stringent adjustment of base quality scores (Koboldt, Larson, & Wilson, 2014 We calculated mean p N (number of nonsynonymous polymorphisms over nonsynonymous sites) and mean p S (number of synonymous polymorphisms over synonymous sites) for each gene category as the ratio of the sum of the number of polymorphisms to the sum of the number of sites (p N = sum P N /sum N; p S = sum P S /sum S).

| Analysis of synonymous codon usage bias
Codon usage bias was estimated using CODONW (http://codonw.sour ceforge.net) through the effective number of codons (ENC) (Wright, 1990 (Wright, 1990). This measure is not biased by the different lengths of the coding regions being analysed, and as such, it has been shown to be more reliable than other commonly used methods of estimating codon usage bias (Comeron & Aguad e, 1998). The effective number of codons was calculated for all the genes with divergence and polymorphism data (Table 2).

| Tests of positive selection
To identify genes evolving under adaptive evolution, we used the  (Andolfatto, 2008;Begun et al., 2007). For genes with significant deviations in D N , D S , P N and P S , a higher nonsynonymous-to-synonymous substitutions ratio relative to polymorphisms ratio (d N /d S > p N /p S ) represented a signature of positive selection. We then tested for significant differences between sex-biased and unbiased genes in the proportion of genes with signatures of positive selection using Fisher's exact test.
For each gene category, we also used the divergence and poly-  We first assessed transcriptional similarity across tissues and sexes using hierarchical clustering of gene expression levels (Figure 2). We found that the reproductive tissue was more transcriptionally dimorphic than the vegetative tissue, consistent with studies in many other species (Jiang & Machado, 2009;Mank, Hultin-Rosenberg, Webster, & Ellegren, 2008;Pointer et al., 2013;Yang, Zhang, & He, 2016). Expression for male catkin clustered most distantly from both male and female expressions in leaf. We identified 3,567 genes (43% of all filtered catkin genes) showing sex-biased expression in catkin (log 2 fold change > 1 or < À1, p adj < .05), compared to expression in the vegetative tissue, where we identified only seven (0.09%) leaf sex-biased genes (Figure 3). Even with a more relaxed fold change threshold for defining differentially expressed genes (log 2 fold change > 0.5 or < À0.5, p adj < .05), we still could not identify any additional leaf sex-biased genes. There were also no shared sex-biased genes between the two tissues.

| Dynamics of catkin sex-biased gene expression
Although female-biased genes (n = 1,820) were slightly more numerous than male-biased genes (n = 1,747), the magnitude of differential expression (log 2 FC) for male-biased genes was significantly greater than that for female-biased genes (Wilcoxon rank sum test p < .001). Average male expression for male-biased genes was significantly higher than average female expression for female-biased genes ( Figure 3, Wilcoxon rank sum test p < .001), although male expression for female-biased genes was significantly lower than female expression for male-biased genes (Figure 3, Wilcoxon rank sum test p < .001).
We The paucity of sex-biased genes in the leaf tissue makes it a useful comparison to further assess the sex-specific changes that give rise to male-and female-biased genes. We therefore used leaf expression as the putative ancestral expression state. For the subset of catkin sex-biased genes that also had expression in the leaf tissue, we determined the difference in expression between catkin and leaf across the same fold change thresholds used in Figure 4. For malebiased genes in the catkin, we found significant differences between catkin and leaf expression in both sexes, although to a lesser extent in females ( Figure S1). On the other hand, for catkin female-biased genes, we also observed large differences in male expression between catkin and leaf samples; however, we found little to no female expression changes between the two tissues ( Figure S1).
We further divided catkin sex-biased genes into autosomal (including the pseudoautosomal region of the sex chromosomes) and Z-linked genes. On the autosomes, we found 3,536 sex-biased genes (1,728 male-biased and 1,808 female-biased genes). On the nonrecombining region of the Z chromosome, we found only 31 sexbiased genes (19 male-biased and 12 female-biased genes); however, considering the narrow region of recombination suppression between the sex chromosomes (Pucholt et al., 2017, 3.5-8.8 Mbp), these sex-biased genes represented 44% of the total identified gene content in the nonrecombining sex-chromosome region.

| Rates of evolution
We compared the overall ratios of nonsynonymous-to-synonymous nucleotide substitutions (d N /d S ) between catkin and leaf and found F I G U R E 2 Heatmap and hierarchical clustering of average male (blue) and average female (red) gene expression in catkin and leaf. The heatmap represents all the filtered genes expressed in both tissues (7,257). Hierarchical gene clustering is based on Euclidean distance with average linkage for log 2 RPKM expression for each gene. Numbers at nodes represent the 1,000 replicates percentage bootstrap results DAROLTI ET AL.
| 699 no significant differences between the two (p = .476, significance based on permutation tests with 1,000 replicates). We also did not find a significant difference in the evolution of unbiased genes between the two tissues (p = .056 from permutation tests with 1,000 replicates), likely influenced by the large overlap of genes between them (97% of catkin unbiased genes represent 58% of the unbiased genes expressed in leaf). We found too few significantly sex-biased genes in the leaf tissue to make any statistical comparisons of rates of sequence evolution between catkin and leaf sexbiased genes.
We also compared the ratio of d N /d S between sex-biased and unbiased genes in catkin to test for differences in the rate of evolutionary divergence. Interestingly, we found that on autosomes, although male-biased genes have more amino acid substitutions than both unbiased and female-biased genes, as shown by significantly higher d N values, d N /d S for male-biased genes was significantly lower, indicating slower rates of functional evolution relative to unbiased (Table 1; Table S2) and female-biased genes (p < .001, significance based on permutation tests with 1,000 replicates). Similar results were obtained when we estimated d N /d S from a data set of 1:1 orthologs that excluded cases of co-orthology (Table S2), as well as from omega values resulting from running CODEML branch model 2 in PAML on concatenated sequences of genes in each sex-bias gene category (Table S3). This lower d N /d S ratio is caused in part by a disproportionate increase in synonymous substitutions compared to nonsynonymous substitutions, causing the relationship between d N and d S in male-biased genes to lie further away from direct proportionality than in the case of unbiased genes ( Figure S2).   Female-biased autosomal loci also showed the same pattern as male-biased genes relative to unbiased genes; however, this result was not significant (Table 1; Table S2). On the nonrecombining Z, male-biased genes also show lower rates of evolution compared to unbiased genes; however, this finding was not significant, likely due to the small sample size of male-biased genes (n = 3). In contrast, female-biased Z-linked loci showed accelerated rates of evolution in comparison with male-biased Z-linked genes (p < .001, significance based on permutation tests with 1,000 replicates).
Highly expressed genes are often observed to exhibit lower d N / d S values (Cherry, 2010;Drummond, Bloom, Adami, Wilke, & Arnold, 2005; P al, Papp, & Hurst, 2001;Slotte et al., 2011); therefore, to determine whether expression level might explain our results, we divided sex-biased and unbiased genes into quartiles based on overall expression. As expected, we found that as gene expression level increases, the rate of sequence divergence decreases and this holds true for both sex-biased and unbiased genes ( Figure S3). To further investigate the effect of expression level on the variation in rates of sequence divergence between sex-bias categories, we used a multiple regression analysis to predict d N /d S results based on expression level and degree of sex-bias. For defining the degree of sex-bias, genes were classed into five groups, highly female-biased genes (FC ≤ À3), lowly female-biased genes (À3 < FC ≤ À1), unbiased genes (À1 < FC < 1), lowly male-biased genes (1 ≤ FC < 3) and highly male-biased genes (FC ≥ 3). We found a significant negative relationship between d N /d S values and both average log 2 RPKM expression level (b = À.03, p < .001) and degree of sex-bias (b = À.04, p = .014). There was no significant effect of the interaction between expression level and degree of sex-bias on d N /d S results, suggesting that any differences in the rates of sequence evolution due to sex-bias are independent of the gene expression level for each sex-bias category. Despite these results, the adjusted r 2 was very low (r 2 = .01), indicating that other factors, such as purifying or haploid selection, largely explain the vast majority of sequence divergence results.
We also estimated average levels of synonymous codon usage bias for sex-biased and unbiased genes to determine whether this could explain the differences in the rates of synonymous substitutions between the gene categories. Stronger codon usage bias has been associated with higher gene expression as selective forces act to increase translational efficiency (Duret, 2002;Ingvarsson, 2010).
We estimated codon usage bias for genes in each category through the effective number of codons (ENC), where stronger codon bias was indicated by lower ENC values. The differences in codon bias between the different gene categories were subtle, and the gene frequency spectra for all categories were distributed towards the higher end of the effective number of codons (ENC), hence lower codon usage bias ( Figure S4). However, male-biased genes had significantly lower codon usage bias than both unbiased (Table 2) and female-biased genes (p < .001, significance based on permutation tests with 1,000 replicates). These findings, together with the higher rates of synonymous substitutions in male-biased genes compared to unbiased and female-biased genes, indicate weaker purifying selection on silent mutations in male-biased genes (Sharp & Li, 1987).
We used polymorphism data to calculate the ratio of nonsynonymous-to-synonymous polymorphisms (p N /p S ). Sex-biased genes on both autosomes and the nonrecombining Z region have significantly higher nonsynonymous and synonymous polymorphism levels compared to unbiased genes; however, the p N /p S ratio was not significantly different in either of the comparisons (Table 1). To distinguish between the selective pressures acting on sequence evolution, we used the McDonald-Kreitman test of selection, comparing the ratios of d N /d S to p N /p S for each gene category. Following filtering, we recovered six unbiased, one male-biased and two female-biased genes showing signatures of positive selection (Table 3). However, there was no significant difference in the proportion of genes evolving under positive selection between either of the gene categories (   (Gossmann et al., 2014).
This could be the result of the greater haploid selection in plants; however, it could also be, at least partially, the result of the selfing mating system in this species, which leads to the purging of recessive deleterious variation. Similarly, in the self-incompatible close relative of A. thaliana, C. grandiflora, a larger fraction of pollen-specific genes was found to evolve under strong purifying selection and to also exhibit faster protein evolution rates compared to sporophytic genes (Arunkumar et al., 2013). This is suggested to be the result of both higher pollen competition and the haploid nature of the pollenspecific tissue.
Here, we investigate the evolution of sex-biased genes in S. viminalis, a perennial dioecious (obligate outcrossing) species with partial wind pollination. Similarly to C. grandiflora (Kao & McCubbin, 1996), S. viminalis theoretically experiences far higher levels of pollen competition than A. thaliana, particularly intermale competition. Although we might expect the high levels of sperm competition in S. viminalis to produce higher rates of protein evolution for male-biased genes, we observed the opposite. Moreover, in contrast to work in C. grandiflora (Arunkumar et al., 2013), we did not find evidence of a high proportion of male-biased genes under positive selection.
The observed dynamics of sex-biased gene expression in S. viminalis is consistent with previous reports in a wide range of species.
Equivalent to studies on somatic and reproductive tissues in animal systems (Mank, 2017;Pointer et al., 2013;Yang et al., 2016), we found that the reproductive tissue was far more transcriptionally dimorphic than the vegetative tissue (Figures 2 and 3). Additionally, in plant species in particular, very few studies have been able to identify any significant sex-biased genes in nonreproductive tissues (Robinson et al., 2014;Zemp, Minder, & Widmer, 2014;Zluvova et al., 2010). We also found that, in catkin, male-biased genes were expressed at significantly higher levels and had a higher magnitude of sex-bias than female-biased genes (Figure 3). The level of sexbiased gene expression found in the S. viminalis reproductive tissue is markedly lower than that in animal species (Jiang & Machado, 2009;Pointer et al., 2013), consistent with the significantly higher degree of sexual dimorphism in animal systems. On the other hand, we found a larger percentage of sex-biased genes compared to several plant and algae species with low levels of sexual dimorphism (Harkess et al., 2015;Lipinska et al., 2015;Zemp et al., 2016). This is indicative of higher intersexual morphological differences in the S.
viminalis reproductive tissue, which is consistent with previous descriptions of the structural differences between male and female catkins (Cronk et al., 2015).
Contrary to findings from the dioecious Silene latifolia (Zemp et al., 2016), however similarly to reports from animal and algae systems (Lipinska et al., 2015;Perry et al., 2014), our results indicate that sex-biased gene expression has likely evolved as an outcome of expression changes in males ( Figure S1). This would also explain why catkin male samples are more transcriptionally different than catkin female samples with respect to leaf samples ( Figure 2). These results suggest that ancestral intralocus sexual conflict may have been more detrimental to males, leading to the evolution of sex-biased gene expression in order to resolve such conflicts.
Additionally, although not statistically significant, we found that male-biased genes had higher p N /p S values compared to both unbiased and female-biased genes, which is in stark contrast to divergence results where we found male-biased genes to have significantly lower d N /d S values. Given that perturbations in population size can alter estimates of polymorphism (Pool & Nielsen, 2007;Tajima, 1989), it is difficult to assess the causes of the contrasting results between d N /d S and p N /p S estimates for sex-biased genes. Nevertheless, divergence estimates are less sensitive to demographic fluctuations and we more strongly rely on this measurement in our analyses of evolutionary rates of sex-biased genes.
Sex-biased genes in willow exhibit higher expression levels than unbiased genes, and highly expressed male-biased and female-biased genes had significantly lower rates of evolution than unbiased and lowly expressed sex-biased genes ( Figure S3). The fact that highly expressed genes evolve more slowly could be due to a range of different reasons, which are still highly debated (Drummond et al., 2005). The structural or functional features of the proteins they encode (Drummond et al., 2005), high pleiotropic constraints acting on the genes (P al et al., 2001) as well as gene conversion events (Petes & Hill, 1988) have all been suggested as potential mechanisms through which highly expressed genes could have lower rates of sequence evolution. Although the high expression of many sexbiased genes in S. viminalis may partially explain their slower rates of evolution, our analysis revealed a very weak correlation between expression level and rate of evolution, indicating that, in this case, expression level does not largely explain the low rates of sex-biased gene evolution.
It is interesting that the lower d N /d S values of male-biased genes are associated with an overall increase in synonymous mutations DAROLTI ET AL.
Haploid selection is far more extensive in plants due to both the larger proportion of the life cycle spent in the haploid phase and active gene transcription, which has been observed in gametes, particularly in pollen (Otto et al., 2015). In addition to haploid selection, male gametophytes in angiosperm species are under strong sexual selection pressures (Erbar, 2003;Snow & Spira, 1996), particularly in outcrossing species. Mechanisms of sexual selection in angiosperms include pollen tube and pistil interactions and pollen competition over ovules, which is exacerbated in outcrossing species (Bernasconi et al., 2004).
It is important to note that the reduced floral structure and microscopic nature of the catkin (Cronk et al., 2015) makes it nearly impossible to separate haploid from diploid reproductive tissue in this species. However, our catkin preparations are highly enriched for haploid cells (Figure 1) when compared to the vegetative samples. We expect that rates of evolution for purely haploid sex-biased tissue would be even lower than what we observe if haploid selection is indeed the primary cause of the slower rates of evolution.
Apart from insect pollen dispersal, willows also have wind-dispersed pollination (Peeters & Totland, 1999) and experience high levels of pollen competition. The observed patterns of gene sequence evolution in S. viminalis support the notion that pollen competition in conjunction with haploid selection produces greater levels of purifying selection on male-biased genes. This would remove deleterious variation and lead to significantly slower rates of functional gene sequence evolution. Interestingly, the algae Ectocarpus, a species where sex-biased genes are subject almost entirely to haploid selection, shows accelerated rates of evolution for both male-and female-biased genes (Lipinska et al., 2015). This suggests that haploid selection may not be the only force that influences the rate of evolution of sex-biased genes in haploid cells. Indeed, data from haploid-specific genes (pollen-specific genes in S. viminalis) would help to more precisely determine the degree to which the currently observed lower rates of evolution of male-biased genes can be explained by haploid selection or other factors such as expression breath (Arunkumar et al., 2013;Gossmann et al., 2014;Sz€ ov enyi et al., 2013).
In summary, our findings are generally consistent with previous reports on the patterns of sex-bias gene expression in plant and animal species. However, different forces may differentiate patterns of evolution between animal and plant systems. The reduction in haploid selection in animals may limit the power of purifying selection to remove mildly deleterious variation, particularly when it is largely recessive. In S. viminalis, we observe reduced rates of evolution for male-biased genes, consistent with increased purifying selection from the extended haploid phase. Even though male-biased genes show relaxed levels of codon bias, this does not seem to be a major driver of the reduced rate of evolution. Future work should focus on investigating the differences in the relative strength of haploid versus diploid selection in dioecious angiosperm species in shaping the evolution of sex-biased genes.

ACKNOWLEDG EMENTS
Sequencing was performed by the SNP&SEQ Technology Platform in Uppsala, Sweden. The facility is part of the National Genomics