Sexually dimorphic gene expression and transcriptome evolution provide mixed evidence for a fast‐Z effect in Heliconius

Abstract Sex chromosomes have different evolutionary properties compared to autosomes due to their hemizygous nature. In particular, recessive mutations are more readily exposed to selection, which can lead to faster rates of molecular evolution. Here, we report patterns of gene expression and molecular evolution for a group of butterflies. First, we improve the completeness of the Heliconius melpomene reference annotation, a neotropical butterfly with a ZW sex determination system. Then, we analyse RNA from male and female whole abdomens and sequence female ovary and gut tissue to identify sex‐ and tissue‐specific gene expression profiles in H. melpomene. Using these expression profiles, we compare (a) sequence divergence and polymorphism; (b) the strength of positive and negative selection; and (c) rates of adaptive evolution, for Z and autosomal genes between two species of Heliconius butterflies, H. melpomene and H. erato. We show that the rate of adaptive substitutions is higher for Z than autosomal genes, but contrary to expectation, it is also higher for male‐biased than female‐biased genes. Additionally, we find no significant increase in the rate of adaptive evolution or purifying selection on genes expressed in ovary tissue, a heterogametic‐specific tissue. Our results contribute to a growing body of literature from other ZW systems that also provide mixed evidence for a fast‐Z effect where hemizygosity influences the rate of adaptive substitutions.

Patterns of molecular evolution on sex chromosomes are particularly influenced by gene expression patterns. Sexually dimorphic expression is often caused by natural and/or sexual selection favouring phenotypes that influence the fitness of one of the sexes (Grath & Parsch, 2016). In species with genetic sex determination, the majority of sexually dimorphic traits results from the differential expression of genes present in both male and female genomes (Ellegren & Parsch, 2007). Sex-biased expression is common across taxa from mammals (Rinn & Snyder, 2005) to Diptera (Assis, Zhou, & Bachtrog, 2012), reptiles (Cox et al., 2017), birds (Mank, Nam, & Ellegren, 2010;Mank, Vicoso, Berlin, & Charlesworth, 2010) and Lepidoptera (Rousselle, Faivre, Ballenghien, Galtier, & Nabholz, 2016). For example, in Drosophila melanogaster, 57% of genes have been categorized as sex biased (Assis et al., 2012), and, in Heliconius melpomene, analysis of two different tissues identified up to 29% of expressed genes as sex biased (Walters, Hardcastle, & Jiggins, 2015). The vast majority of genes that exhibit sexually dimorphic expression are active in reproductive tissues and tend to also have distinctive rates of molecular evolution compared to genes without dimorphic expression (Avila, Campos, & Charlesworth, 2015;Parisi et al., 2003Parisi et al., , 2004. Ultimately, the identification of sex-biased genes, and subsequent analysis of patterns of molecular evolution, will contribute to a better understanding of the evolutionary forces shaping sex chromosome and autosome evolution (Assis et al., 2012;Kirkpatrick & Hall, 2004;Zhang, Hambuch, & Parsch, 2004).
Empirical studies of the fast-X effect typically measure two different metrics: (a) the ratio of nonsynonymous to synonymous substitution rates (dN/dS) and (b) the amount of adaptive evolution (α) using the McDonald-Kreitman (MK) test (McDonald & Kreitman, 1991). Studies measuring dN/dS usually test for "faster-X divergence." Although this approach may be useful for comparing sex chromosome and autosomal divergence, measuring the relative rate of nonsynonymous substitutions captures the effects of both adaptive and neutral (or slightly deleterious) mutations. Estimates of α can better test for an excess of adaptive substitution in the sex chromosome ("faster-X adaptation") by combining measures of withinspecies polymorphism and between-species divergence, but α is still sensitive to the rate of accumulation of slightly deleterious mutations and demography (Fay, 2011). For instance, an increase in N e is expected to result in decreased dN/dS and increased α even when the rate of adaptive substitutions remains unchanged. To overcome this problem, extensions of the MK test such as ω a were developed to estimate the rate of adaptation by calculating the frequency distribution of polymorphism after correcting for demographic history and distribution of deleterious effects at functional sites (Galtier, 2016).
The analysis of evolutionary rates between sex and autosomal genes, however, has produced mixed evidence in support of fast-X evolution (Meisel & Connallon, 2013). In some taxa, there is strong evidence for faster-X divergence but not faster-X adaptation, or vice versa (Meisel, Malone, & Clark, 2012). For example, the first calculations of faster-X divergence were carried out in Drosophila where support for elevated dN/dS in X genes has been mixed. Studies that used autosometo-X translocations to control for gene content effect did not reach a consensus on the existence of faster-X divergence Thornton, Bachtrog, & Andolfatto, 2006;) but X-linked duplicate genes have elevated dN/dS compared to autosomal duplicates (Thornton & Long, 2002).
Signals of faster-X sequence divergence in Drosophila have been shown to affect noncoding regulatory regions as well, and might be at least partly explained by differences in gene composition on the X versus the autosomes (Hu, Eisen, Thornton, & Andolfatto, 2013). However, faster-X divergence in other taxa has received stronger support. For example, in humans, chimpanzees and rodents, dN/dS is higher for X genes (Mank, Vicoso et al., 2010;Nielsen et al., 2005).
In contrast, whole-genome analyses of adaptive substitutions have resulted in stronger evidence for faster-X adaptation in Drosophila (Mackay et al., 2012), whereas support for faster-X adaptation in vertebrates is less clear. McDonald-Kreitman tests support faster-X adaptation in wild mouse populations (Baines & Harr, 2007) but, for the European rabbit (Oryctolagus cuniculus), a clear faster-X adaptation signal is only present in populations with large effective population sizes (Carneiro et al., 2012).
Taxa with ZW sex determination provide an interesting contrast.
For female heterogametic taxa such as birds, females only have one copy of the Z chromosome. A fast-Z effect may be expected to result from the expression of recessive mutations on the Z chromosome as Z genes are immediately exposed to selection in females (Charlesworth, 2012). In birds, fast-Z divergence has been reported, but Z male-biased genes were not less accelerated than unbiased genes or female-biased genes (Wright, Zimmer, Harrison, & Mank, 2015). This would not be expected if the fast-Z effect was driven by recessive beneficial mutations, and so, it was suggested that fast-Z in birds does not reflect positive selection (Mank, Nam, et al., 2010;Wright et al., 2015).
For Lepidoptera, results have also been mixed. Sackton et al. (2014) reported that faster-Z evolution was driven by positive selection in silkworms. But, in satyrine butterflies, there were no significant differences in adaptive evolutionary rates between the Z and the autosomes (i.e., no fast-Z adaptation). However, the comparison of male-biased, female-biased and unbiased Z genes in satyrine butterflies revealed increased purifying selection against recessive deleterious mutations in female-biased Z genes (Rousselle et al., 2016). Therefore, considerable uncertainty remains regarding the prevalence and magnitude of the fast-X/Z effect on divergence and adaptation.
Here, we investigate the effects of hemizygosity on the rates of adaptive substitution in the neotropical butterfly genus Heliconius, a ZW sex determination system, by analysing polymorphism, divergence and gene expression genomewide. We test whether there is a fast-Z effect in Heliconius using two species from the H. melpomene and H. erato clades which diverged 12 million years ago (synonymous divergence = 0.15) (Kozak et al., 2015;Martin et al., 2016). Previous analyses of Heliconius transcriptome data have focused on the evolution of dosage compensation and the impact of sex-specific dosage on the levels of gene expression (Walters et al., 2015). In this study, using the same transcriptome data, we first compute sex-biased expression.
Then, accounting for sex-biased gene expression, we (a) calculate coding sequence divergence and polymorphism in H. melpomene and (b) assess the strength of positive and negative selection, and rates of adaptive evolution between H. melpomene and H. erato. We then analyse newly generated female transcriptome data from H. melpomene ovary and gut tissue in order to investigate whether genes expressed in the reproductive tissue of the heterogametic sex have higher rates of adaptive evolution than those expressed in somatic tissues.
Moreover, BRAKER1 predictions that had 1-to-1 overlaps with Hmel2 names were replaced by their original Hmel2 name. For many-to-1 mapping between the BRAKER1 predictions and Hmel2, Hmel2 names were reused and a suffix of g1/g2/g3/etc. was added. The rest were renamed from HMEL030000 onwards.
For these 25 samples, H. m. rosina females were reared in insectaries in Gamboa, Panama. Passiflora platyloba potted plants were monitored daily and 5th-instar caterpillars were removed and taken to the laboratory in large individual containers where they were allowed to pupate and emerge at a constant temperature (24-25°C). The pupating containers in the laboratory were monitored several times a day. Sexually antagonistic pressures are expected to be greater in mature adults, and genes exhibiting extreme sexually dimorphic patterns should be expressed in mature individuals (Gibson, Chippindale, & Rice, 2002). In addition, mating induces behavioural and physiological changes and has been shown to trigger regulatory changes in sex-biased genes (e.g., Immonen, Sayadi, Bayram, & Arnqvist, 2017).
By collecting gene expression data from ovaries at two different lifehistory time points, we aimed to increase the number of sex-specific genes identified. When a female emerged, it was either (a) returned to the insectaries to be mated to a H. m. rosina male (Treatment: old, Supporting Information Table S1) or (b) dissected 1 hr after eclosion under controlled laboratory conditions (Treatment: young, Supporting Information Table S1). Mated females were kept in individual 1 m × 1 m × 2 m cages for 20 days until dissection.
Guts and ovaries were dissected in RNAlater (ThermoFisher, Waltham, MA) at 24-25°C, and tissue was stored in RNAlater at 4°C for 24 hr and −20°C thereafter. Total RNA was extracted with a combined guanidium thiocyanate-phenol-chloroform and silica matrix protocol using TRIzol (Invitrogen, Carlsbad, CA), RNeasy columns (Qiagen, Valencia, CA) and DNaseI (Ambion, Naugatuck, CT). mRNA was isolated from total RNA via poly-A pull-down, and directional cDNA library preparation and sequencing (Illumina HiSeq 2500,
We filtered the results as in Walters et al. (2015) with FDR < 0.05 (alpha = 0.05) (Walters et al., 2015). Ovary-and gut-biased genes have a log 2 -fold change significance threshold >1.5 (option: lcf-Threshold < 1.5, altHypothesis = "greaterAbs"). We defined male, female and unbiased genes as in Rousselle et al. (2016). First, we calculated reads per kilobase per million (RPKM) as: where N c is the number of reads mapped to the genic feature, N tot is the total number of reads mapped in the sample, and L c is the length of the genic sequence in base pairs (Mortazavi, Williams, McCue, Schaeffer, & Wold, 2008). RPKM i is the mean RPKM of gene i across the 10 individuals. Genes for which RPKM f /RPKM m > 1.5 were classified as female biased, genes for which RPKM f /RPKM m < 0.66 were classified as male biased, and the others were classified as unbiased (Rousselle et al., 2016).

| Extraction of orthologous genes, coding sequence alignment and SNP calling
OrthoFinder was used to identify orthologous groups of genes in the H. melpomene and the H. erato transcriptomes (options: -t 48 -a 6). 1-1 orthologous gene sequences were selected for use in subsequent analysis (Supporting Information Table S2). Using Gff-Ex, a genome feature extraction package (Rastogi & Gupta, 2014), we

| Calculation of diversity and selection statistics for 1-1 ortholog alignments between H. melpomene and H. erato: Classic approach
The adaptive substitution rate was estimated by comparing synonymous and nonsynonymous variation in the polymorphism and divergence compartments, as first proposed by McDonald & Kreitman, 1991; see also Bustamante et al., 2005, andMacpherson, Sella, Davis, &Petrov, 2007). We first used the original MK test (referred to as Classic approach hereafter) to estimate the rate of adaptive substitution for all genes found to be orthologous between H. melpomene and H. erato. We calculated (a) synonymous polymorphism (P s ) and (b) nonsynonymous polymorphism (P n ) in H. melpomene, as well as (c) synonymous fixed divergence (dS), and (d) nonsynonymous fixed divergence (dN) between H. melpomene and H. erato. We estimated the rate of adaptive molecular evolution (α) between the two species as: α assumes that nonsynonymous mutations are either adaptive, neutral or strongly deleterious (McDonald & Kreitman, 1991), with −∞ > α ≥ 1, where α = 0 represents the null hypothesis that nonsynonymous mutations are neutral (dN/dS = P n /P S ). α > 0 corresponds to dN/dS > P n /P S and indicates positive selection, whereas α < 0 corresponds to dN/dS < P n /P S and indicates negative selection. These values were calculated using the EggLib C++ function polymorphismBPP (v2.1.11) (De Mita & Siol, 2012) and Bio++ (v2.2.0) (Dutheil & Boussau, 2008) in python (v2.7.5) using scripts adapted from https://github.com/tatumdmortimer

| Calculation of diversity and selection statistics for 1-1 ortholog alignments between H. melpomene and H. erato: Modelling approach
The Classic approach to the MK test is robust to differences in mutation rates and variation in coalescent histories across genomic locations (McDonald & Kreitman, 1991). Inference of positive selection using the Classic approach of the MK test is not robust, however, to the occurrence of slightly deleterious mutations and demographic change. To account for these confounders, we used a Modelling approach to estimate the strength of positive and purifying selection in addition to the Classic approach described above, using the method of Eyre-Walker and Keightley (2009) and the per-mutation rate of nonadaptive substitutions is calculated as:

| Gene expression level and π n /π s ratios
To test whether gene expression level and chromosome type have a significant effect on π n /π s ratios, we used a multiple regression analysis. We established the linear model: using r (v3.2.5). Chromosome type is either autosome or sex chromosome as assigned in Hmel2 reference genome (Davey et al., 2016).
477 genes with no polymorphism were removed from the analysis.
We plotted diagnostic plots of residuals versus fitted values.

| RNA-sequencing and read mapping
Analysis of gene expression profiles in the data retrieved from  Figure S1).  Table S1). We analysed data from two different time points and from non-sexand sex-specific tissue separately (Treatment: Young and Old).
There is a clear separation of the 25 samples by tissue when we compare gene expression profiles between them. In total, 51% of the total variance is explained by the two-first principal components. PC1 separates the samples by tissue and explains 40% of variance. PC2 explains 11% of total variance and separates samples by age (Supporting Information Figure S3). Ovarian tissue clusters by age more tightly than non-sex-specific tissue (Gut) (Supporting Information Figure S4A,B).
In order to test whether there was a fast-Z effect in Heliconius, we calculated the following statistics for autosomal and Z genes: dN/dS, α, ω a and ω na . We also calculated π sZ /π sA and investigated the relationship between π n and gene expression. dN/dS tests for "faster-Z divergence." If there is a faster-Z effect in a female heterogametic system, genes with female-biased expression will have higher dN/dS than male-biased or unbiased genes. α, ω a and ω na test for "faster-Z adaptation." In a female heterogametic system, faster-Z adaptation would predict that the proportion of adaptive nonsynonymous substitutions (α) in female-biased genes will be significantly higher in the Z chromosome compared to the autosomes. ω a measures the per-mutation rate of adaptive substitutions and, if there is "fast-Z adaptation," ω a is also predicted to be significantly higher for female-biased Z genes. However, in a female heterogametic system with reduced purifying selection on the Z, the per-mutation rate of nonadaptive substitution (ω nα ) would also be higher for female-biased Z genes. Another metric that can be used as an indicator of the strength of purifying selection is the π sZ /π sA ratio. In a female heterogametic system, stronger purifying selection on the Z chromosome would lead to a π sZ /π sA ratio <0.75 due to background selection. Overall, patterns of diversity tend to be associated with gene expression levels. π n is expected to be negatively correlated with expression levels if there is increased purifying selection on highly expressed genes.

| Coding sequence divergence does not support a significant fast-Z effect
We first compared rates of Z and autosomal sequence divergence log π nij ∼ log π nsj + chromosome_type j + log RPKM More highly expressed genes are more exposed to selection, so in a female heterogametic system with a fast-Z effect, genes with female-biased expression are expected to have higher rates of amino acid substitution if dN/dS is driven by positive selection.

| π sZ /π sA diversity ratio is lower than 0.75
Next, we explored patterns of within-species diversity as an indicator of the strength of purifying selection. In a population at equilibrium with a 1:1 sex ratio, the π sZ /π sA diversity ratio is expected to be 0.75, but stronger purifying selection on the Z chromosome would lead to a reduction in this ratio due to background selection. The π sZ /π sA ratio for H. melpomene is approximately 0.44 (Table 2), which might indicate purifying selection on the Z. However, this ratio can also be influenced by a biased sex ratio (Vicoso & Charlesworth, 2006), differences in recombination rates (Charlesworth, 2012), sex-biased mutation rates (Vicoso & Charlesworth, 2009)

| Increased strength of purifying selection on highly expressed genes
Patterns of diversity were, however, strongly associated with expression levels. Using a multiple regression approach, we found that functional genetic diversity, π n, was significantly negatively correlated with expression level for both autosomal and Z genes (p < 0.01) consistent with increased purifying selection on highly expressed genes (Supporting Information Figure S2) (Figure 1).

| Z and autosomal rates of adaptive substitution: testing fast-Z adaptation
We next explored patterns of adaptive evolution using (  Notes. π n /π s , dN/dS ratios, α, ω a and ω na were calculated for autosomal and Z male-biased, female-biased and unbiased genes. π n /π s , dN/dS ratios were calculated for H. melpomene samples, and α, ω a and ω na were calculated between H. melpomene and H. erato. Intervals represent 95% confidence intervals obtained by bootstrapping 1,000 times. Bold (A) denotes significant values within either Z or autosomal categories. Bold (B) denotes significant values within both Z and autosomal categories. Significance indicated separately for All and for sex-biased expression (female, male and unbiased).
There are no significant differences in α values between Z and autosomal genes under the Classic approach (Table 1) Rice (1984), where the accumulation of alleles under sex-specific and antagonistic selection on the Z was expected if the alleles were recessive in females and dominant in males (Rice, 1984). The lack of significant differences in α between sex-biased genes is not consistent with the expectations of fast-Z adaptation, which would predict faster evolution of femalebiased genes due to hemizygosity compared to autosomes, but this could also reflect a lack of power to detect the signal when the total number of genes is reduced.

| Hemizygosity and the rate of adaptive substitutions
There was no evidence for reduced purifying selection on the Z chromosome, as the per-mutation rate of nonadaptive substitution (ω nα ) is 95% CI = [0.038-0.044]) genes, which confirms the low π n /π s already reported and would suggest that purifying selection is stronger in female-biased genes.

| Female ovary-biased and gut-biased genes
Next, we explored the expression of genes in female reproductive tissue.
Overall there were a greater number of genes with gut-biased expression (#Gut Auto = 153) than ovary-biased expression (#Ovary Auto = 40) in the autosomes. However, there was an over-representation of Z ovaryexpressed genes than expected by chance (#Gut Z = 6; #Ovary Z = 6; chisquare test; p < 0.05). However, the number of genes in each category is relatively small so these tests should be treated with caution.
Of the 205 differentially expressed genes between the two tissues, only 9 in the ovaries and 26 in the gut could be used to calculate dN/dS, π n /π s and α. The other genes either do not have a 1-1 ortholog with H. erato or there were too many undetermined characters (gaps or Ns) to be able to estimate the parameters. Of the 35 genes for which molecular evolution statistics could be calculated, all 9 ovarybiased and 25 of 26 gut-biased genes are autosomal; and 1 gutbiased gene maps to the Z (Gut dN/dS Z = 0.365; Gut π n /π s Z = 0.093; Gut α Z = 0.295). We did not detect any significant differences in π n / π s ; dN/dS; or α for autosomal ovary-and gut-biased genes.

| D ISCUSS I ON
Elevated rates of coding sequence evolution on the sex chromosome relative to autosomes have been reported for several species, consistent with the theoretical prediction of fast-X evolution. Here, we find evidence for enhanced rates of adaptation on the Heliconius Z chromosome: Z genes have a significantly higher rate of adaptive TA B L E 2 H. melpomene π s and π sZ /π sA ratio from pairwise alignments for Z and autosomal genes Notes. π s calculated from pairwise alignments for Z and autosomal genes. π sZ /π sA ratio used to estimate N eZ /N eA . Intervals represent 95% confidence intervals obtained by bootstrapping genes (1,000 replicates). Bold (A) denotes significant values within either Z or autosomal categories. Significance indicated separately for All and for sex biased (female, male and unbiased).
F I G U R E 1 Expression level of Z and autosomal genes. Median expression level of Z genes is significantly lower than autosomal genes (p < 0.05). Notches on boxplot display the confidence intervals around the median evolution when all expressed genes are considered. However, fast-X theory predicts that genes highly expressed in the hemizygous sex should be especially prone to fast-X evolution, and this prediction was not satisfied in our data. Female-biased genes did not evolve faster when located on the Z chromosome. The evidence for fast-Z evolution in Heliconius is, therefore, mixed.
In other taxa, there is strongest support for fast-X evolution in groups with complete dosage compensation (Mank, Vicoso, et al., 2010;Meisel & Connallon, 2013). Theory predicts that opportunities for fast-X evolution should increase in species with somatic X-inactivation such as eutherian mammals, as there is effectively haploid expression of the sex chromosome in cells, increasing the chances of recessive beneficial mutations being fixed (Charlesworth, 2012). Groups such as Lepidoptera have been reported to have more complex patterns of sex chromosome dosage compensation. In Heliconius males, expression of Z genes is reduced below autosomal levels, but this dosage compensation mechanism is imperfect, with males showing increased expression relative to females on Z chromosome genes (Walters et al., 2015). However, the apparent incomplete dosage compensation could be a consequence of an uneven distribution of sex-biased genes on sex chromosomes (Gu & Walters, 2017, Huylmans, Macon, & Vicoso, 2017. Regardless, when we compare rates of divergence and adaptation for genes with sex-biased expression, the expectations of fast-Z evolution are not clearly met. Although we might expect faster rates of adaptive evolution for female-biased genes, we observe a weak tendency for faster rates of evolution in male-biased genes. This might mean that the Z could in fact be a hotspot for dominant alleles that benefit males. If this is the case, faster rates of evolution in male Z genes could contribute to a fast-Z effect that is unrelated to hemizygous nature of the Z in females (Rice, 1984). It is important to note that sex-specific genes can evolve rapidly both due to (a) neutral processes (as they experience relaxed selection in the other sex) and (b) the accumulation moderately deleterious mutations (Dapper & Wade, 2016;Gershoni & Pietrokovski, 2014Mank, 2017). Nevertheless, our results of ω α for male-biased Z-linked genes seem to support faster adaptive evolution of such genes: as the Z spends less time in females as compared to autosomes, relaxed selection should be greater for autosomal male-biased genes (Rice, 1984).
Although a fast-Z effect has been observed in Bombyx mori (Sackton et al., 2014), no such pattern was reported in two satyrine butterflies where the dN/dS ratio of Z genes was slightly lower than for autosomal genes (Rousselle et al., 2016). In Heliconius, although dN/dS was not significantly different between autosomal and Z genes, we did find evidence for a faster rate of adaptive substitution.
Interestingly, our data also show that dS on the Z chromosome is higher than dS on the autosomes perhaps indicating a male-biased mutation rate, as Z chromosomes spend more time in males than in females (Miyata, Hayashida, Kuma, Mitsuyasu, & Yasunaga, 1987).
Although hemizygosity is expected to expose beneficial mutations to selection and increase rates of adaptive evolution on the Z chromosome, it is also expected to increase the efficacy of purifying selection, which would act to reduce evolutionary rates. It may be that the balance between these two forces differs across lepidopteran species, leading to the mixed pattern of fast-Z evolution in some taxa but not others. It is important to add, however, that the α values estimated in this study are substantially higher than α reported in Martin et al. (2016). Martin et al. (2016) estimated α using the approach developed by Messer and Petrov (2013), and, in simulations, it has been shown that it is possible that there is an overestimation of DFEα (the method used in this study) in scenarios with strong sweeps or population expansion (Messer & Petrov, 2013). Wright et al. (2015) interpreted the high dN/dS in Z genes of birds as a consequence of reduced effective population size rather than positive selection. The difference in effective population size between sex chromosomes and autosomes in female heterogametic systems is predicted to be larger than in male heterogametic systems due to higher variance of male reproductive success (Mank, Nam et al., 2010). Indeed, we estimate that coding regions on the Z chromosome have an N e 0.44 times that of autosomes. We might therefore predict a considerable reduction in the efficacy of purifying selection on butterfly Z chromosomes. This should lead to higher ω na and π n /π s ratios on the Z due to stronger genetic drift. However, as in satyrine butterflies (Rousselle et al., 2016), in Heliconius, ω na is not higher on the Z relative to autosomes. In Heliconius, dN/dS and π n /π s are higher in the Z compared to autosomes, but this is not significant. This means that, in contrast to birds, the difference in the effective population size of the Z relative to autosomes is not sufficient to reduce the efficacy of purifying selection at a detectable level.
One possible explanation for this difference is the generally much higher effective population sizes of Lepidoptera, which could allow for efficient selection even on sex chromosome (Rousselle et al., 2016). Another is that by not using all genomic sites to estimate the π sZ /π sA diversity ratio, we might be underestimating its true value due to a stronger effect of background selection. The latter is supported by the observation that, in a recently published paper using all genomic sites to estimate π sZ /π sA , Van Belleghem et al. Another factor that might counteract the fast-Z effect is adaptation from standing variation. Larger populations are more polymorphic and, therefore, have an increased probability of adaption from standing genetic variation. Adaptation from standing genetic variation is expected to result in faster autosome evolution, independent of the dominance of beneficial alleles (Orr & Betancourt, 2001), which would counteract the fast-Z effect. This may be especially relevant when overall population sizes are large, as in many Heliconius species, such that standing variation becomes a comparatively important source of adaptive variation compared to de novo mutations.
As sex genes tend to be expressed in sex-specific tissue such as the testis and the ovaries, we aimed to investigate patterns of molecular evolution in ovary-biased genes. Unfortunately, there are no ovarybiased genes with 1-1 orthologs between H. melpomene and H. erato that are in the Z. This means we could not test the effect of hemizygosity on non-sex-specific and sex-specific female expression directly.
The lack of 1-1 orthology may mean that these genes are rapidly evolving. On the one hand, adaptive evolution rates could be higher. And, indeed, autosomal ovary-expressed genes have higher rates of adaptive evolution than gut expressed genes. On the other hand, the lack of 1-1 orthology could also be due to rapid nonadaptive evolution (Gershoni & Pietrokovski, 2014).
It has been recently shown that, like in other taxa, sex-biased genes can experience rapid turnovers in butterflies (e.g., Papilio) (Huylmans et al., 2017 (Kozak et al., 2015) compared to 35MYA for P. xuthus and P. machaon (Zakharov, Caterino, & Sperling, 2004 Council. We thank the editor and two anonymous reviewers for their comments that helped us to improve this manuscript.

DATA ACCE S S I B I LIT Y
Sequencing data are available from the European Nucleotide Archive