Hemizygosity Enhances Purifying Selection: Lack of Fast-Z Evolution in Two Satyrine Butterflies.

The fixation probability of a recessive beneficial mutation is increased on the X or Z chromosome, relative to autosomes, because recessive alleles carried by X or Z are exposed to selection in the heterogametic sex. This leads to an increased dN/dS ratio on sex chromosomes relative to autosomes, a pattern called the "fast-X" or "fast-Z" effect. Besides positive selection, the strength of genetic drift and the efficacy of purifying selection, which affect the rate of molecular evolution, might differ between sex chromosomes and autosomes. Disentangling the complex effects of these distinct forces requires the genome-wide analysis of polymorphism, divergence and gene expression data in a variety of taxa. Here we study the influence of hemizygosity of the Z chromosome in Maniola jurtina and Pyronia tithonus, two species of butterflies (Lepidoptera, Nymphalidae, Satyrinae). Using transcriptome data, we compare the strength of positive and negative selection between Z and autosomes accounting for sex-specific gene expression. We show that M. jurtina and P. tithonus do not experience a faster, but rather a slightly slower evolutionary rate on the Z than on autosomes. Our analysis failed to detect a significant difference in adaptive evolutionary rate between Z and autosomes, but comparison of male-biased, unbiased and female-biased Z-linked genes revealed an increased efficacy of purifying selection against recessive deleterious mutations in female-biased Z-linked genes. This probably contributes to the lack of fast-Z evolution of satyrines. We suggest that the effect of hemizygosity on the fate of recessive deleterious mutations should be taken into account when interpreting patterns of molecular evolution in sex chromosomes vs. autosomes.

Introduction

Sex chromosomes of sufficiently ancient origin are hemizygous, i.e., effectively haploid in one sex (males in X/Y systems, females in Z/W systems) and diploid in the other one. A consequence of hemizygosity is that recessive beneficial mutations occurring in sex-linked genes are immediately exposed to selection in the heterogametic sex. In contrast, in autosomes a recessive mutation is only exposed to selection at homozygous state, which occurs rarely as far as recently appeared mutations are concerned. Recessive beneficial substitutions are therefore expected to accumulate at a faster rate on the X (respectively, Z) chromosome than on autosomes due to a more efficient positive selection in males (respectively, females; Haldane 1924; Rice 1984; Charlesworth et al. 1987). This should result in an increased evolutionary rate of sex chromosomes relative to autosomes, a phenomenon called the «fast-X» («fast-Z») effect. The theoretical prediction of a faster evolution on the X has been empirically corroborated in several species of mammals (Carneiro et al. 2012; Hvilsom et al. 2011; Kousathanas et al. 2014; Veeramah et al. 2014; Nam et al. 2015) and fruit flies (Betancourt et al. 2002; Thornton and Long 2002; Counterman et al. 2004; Mackay et al. 2012; Avila et al. 2014; Campos et al. 2014; Garrigan et al. 2014). This large body of literature is globally consistent with the hypothesis that recessive adaptive mutations are sufficiently common to significantly accelerate molecular evolution on the X.

A couple of recent studies, however, suggest that the situation could be more complex than suggested above. First, Nguyen et al. (2015) analyzed the effect of gene expression level and recombination rate on X-linked genes evolution in mammals—particularly, the elevated ratio of non-synonymous, dN, to synonymous, dS, substitution rates (dN/dS). X-linked genes tend to have a lower expression level than autosomal genes (Marín et al. 2000; Julien et al. 2012), and gene expression is known to be anti-correlated to dN/dS in a broad range of organisms (Drummond et al. 2005; Drummond and Wilke 2008), presumably because of a reduced intensity of purifying selection in low-expressed genes. Nguyen et al. (2015) suggested that this effect is sufficient to explain most of the fast-X effect in primates and murine rodents without invoking positive selection and dominance effects. Second, two studies in the female-heterogametic birds revealed the existence of a fast-Z effect, the dN/dS ratio of the Z chromosome being higher than the dN/dS ratio in autosomes, but found that Z-linked genes predominantly expressed in males (i.e., male-biased) were not less accelerated than unbiased or female-biased genes (Mank et al. 2010; Wright et al. 2015). This is not expected under the hypothesis of a fast-Z driven by recessive beneficial mutations. Being exposed to selection primarily in the homogametic sex, male-biased genes in Z/W systems should not be affected by hemizygosity (Kousathanas et al. 2014; Sackton et al. 2014). For this reason, Wright et al. (2015) did not interpret the increased dN/dS ratio of Z-linked genes as a consequence of a more efficient positive selection, but rather as an effect of a reduced effective population size of the Z chromosome. The difference in effective population size between sex chromosome (NeZ) and autosomes (NeA) in female heterogametic systems is supposed to be larger than in male heterogametic systems, due to a higher variance of reproductive success in males than in females (Mank et al. 2010). Low NeZ/NeA ratios in birds would entail a stronger intensity of genetic drift in the Z chromosome relative to autosomes, leading to a decreased efficacy of purifying selection and consequently an increased probability of fixation of slightly deleterious mutations, and consequently an increased dN/dS ratio (Wright et al. 2015). It was therefore postulated that the fast-X and the fast-Z effects may be driven by distinct evolutionary forces, the former being due to an elevated adaptive substitution rate on the X, and the latter to an elevated rate of slightly deleterious substitutions on the Z (Mank et al. 2010; Wright et al. 2015). However, this interesting hypothesis was challenged by a recent study suggesting that Z chromosomes in silkmoths (Bombyx) experience an accelerated evolutionary rate due to more efficient positive selection (Sackton et al. 2014), similarly to the pattern reported on the X in mammals and Drosophila.

The literature therefore suggests that positive selection in favor of recessive beneficial mutations is not always the predominant factor responsible for differences in evolutionary rates between X/Z chromosomes and autosomes. Purifying selection, genetic drift and gene expression level are important players too, and their respective impacts are difficult to disentangle. Whether the fast-X and fast-Z effects have similar or distinct evolutionary causes, for instance, is currently unclear. It is noteworthy that hemizygosity is predicted to improve the efficacy not only of positive selection, but also of purifying selection against recessive deleterious mutations (Charlesworth et al. 1987). The effect of homozygosity on the fate of deleterious mutations has been theoretically and empirically investigated in the context of comparisons between selfing and outcrossing species (Glémin 2007, Szövényi et al. 2014), but has so far received limited consideration from the empirical literature on sex chromosome evolution (but see Veeramah et al. 2014 and Charlesworth 2012). Only one study reported a greater codon usage bias on the X chromosome relative to autosomes in Drosophila melanogaster and Caenorhabditis elegans, suggesting that purifying selection has greater efficacy on the X than on autosomes (Singh et al. 2005).

In an attempt to address these issues, we analyzed sex-linked vs. autosomal polymorphism, divergence and gene expression pattern in two closely related species of Satyrinae butterflies, the Meadow Brown (Maniola jurtina) and the Gatekeeper (Pyronia tithonus). Butterflies are Z/W species typically carrying a higher level of genetic polymorphism than birds (Romiguier et al. 2014), thus offering a good opportunity to discriminate between purifying and positive selection through McDonald–Kreitman (McDonald and Kreitman 1991) related approaches. M. jurtina and P. tithonus belong to the same family (Nymphalidae) as the fully sequenced Glanville fritillary (Melitaea cinxia), in which a high-resolution linkage map is available (Ahola et al. 2014).

Surprisingly, we report in these species not a fast-Z, but rather a slight slow-Z effect. The analysis of sex-biased genes allowed us to demonstrate that purifying selection against recessive deleterious mutations is more efficient in the Z chromosome, presumably as a consequence of hemizygosity. Moreover, although positive selection is slightly enhanced in Z-linked female-biased genes, we did not detect an overall higher adaptive substitution rate in the Z chromosome than in autosomes. Taken together, these results explain the lack of fast-Z in Satyrinae butterflies, and reveal the complexity of the interactions between the evolutionary forces affecting substitution rates in sex chromosomes vs. autosomes.

Materials and Methods

De Novo Transcriptome Assembly

RNA-Seq data were obtained from ten living adult individuals of each of M. jurtina and P. tithonus sampled from their natural habitat in several locations spread across France, Germany, Spain, and Portugal (supplementary table S1, Supplementary Material online). Wings and genitalia have been conserved for sexing. The rest of the body was used to extract total RNA using standardized protocols (Gayral et al. 2011). Complementary DNA was sequenced on Illumina HiSeq 2000. We assessed the quality of the reads with FastQC v0.10.1 (www.bioinformatics.babraham.ac.uk/projects/fastqc) and we removed reads containing adapters as well as reads shorter than 50 bp. We cut the end of reads in a way that ensures that the mean per base sequence quality was above 30. We constructed de novo transcriptome assemblies for each species following strategies B in Cahais et al. (2012), using a pipeline involving Abyss (Simpson et al. 2009) and Cap3 (Huang and Madan 1999). Reads were mapped to predicted cDNAs (contigs) with the BWA (Burrow Wheeler Aligner) program (Li and Durbin 2009). Contigs with a per-individual average coverage below 2.5X were discarded. Open reading frames (ORFs) were predicted using the Trinity package (Grabherr et al. 2011). Contigs carrying ORF shorter than 150 bp were discarded.

Orthology and Sex-Linkage Prediction

We downloaded annotated coding sequences of the Glanville fritillary (M. cinxia) using biomart (www.biomart.org, January 2015). To predict orthology between M. jurtina and P. tithonus, respectively, and the reference M. cinxia, we performed BLASTn searches. Only hits with an e-value below 10−20 were considered. When one ORF from M. jurtina or P. tithonus hit several reference sequences in M. cinxia, we selected the couple presenting the highest score. When several ORFs from M. jurtina or P. tithonus hit the same reference sequence, we checked whether hits were overlapping, in which case we only kept the one with the highest score. Otherwise, we kept them all. Eventually, we verified those predictions of orthology using another reference species of the Nymphalidae family, Heliconius melpomene (The Heliconius Genome Consortium 2012). Sequences that hit the same gene in M. cinxia but distinct genes in H. melpomene were excluded.

Z-linked/autosomal annotations of genes in M. cinxia (Ahola et al. 2014, http://www.helsinki.fi/science/metapop/research/mcgenome2_downloads.html; last accessed February 2015) were propagated to their predicted orthologs in M. jurtina and P. tithonus. We removed ORFs predicted to be Z-linked based on M. cinxia annotations but presenting heterozygous sites in females of M. jurtina or P. tithonus in the main analysis.

Divergence and Polymorphism Analyses

Orthologous M. jurtina and P. tithonus coding sequences were aligned using MACSE [Multiple Alignment of Coding SEquences accounting for frameshifts and stop codons (Ranwez et al. 2011)] and the alignments were cleaned for gaps (by removing all sites that presented a gap in at least one species). We used the programs bppml and mapNH [http://biopp.univ-montp2.fr/forge/testnh (Romiguier et al. 2012; Guéguen et al. 2013)] to estimate the synonymous and non-synonymous substitution rate (dS, the number of synonymous substitutions per synonymous site and dN, the number of non-synonymous substitutions per non-synonymous site) by substitution mapping under the Nielsen–Yang model (Nielsen and Yang 1998). M. cinxia and H. melpomene coding sequences were too distant from the Satyrinae to allow proper estimation of dS—they were not analyzed. For any set of genes (Z-linked, autosomal, sex-biased), the number of synonymous and non-synonymous substitutions and the number of synonymous and non-synonymous sites were summed across genes and dN/dS was calculated by taking the ratio of the sums. Ninety-five percent confidence intervals were determined by bootstrapping genes (1000 replicates).

For each individual, diploid genotypes were called according to the method described in Tsagkogeorga et al. 2012 (model M1). This method estimates the sequencing error rate in the maximum likelihood framework and calculates the posterior probability of each possible genotype. Genotypes supported by a posterior probability higher than 95% are retained, otherwise missing data is called. Polymorphic positions were filtered for possible hidden paralogues using a likelihood ratio test based on explicit modeling of paralogy (Gayral et al. 2013, Romiguier et al. 2014). Sites with less than six genotyped individual (over ten) were excluded. Non-synonymous (πn) and synonymous (πs) nucleotide diversity were computed for each set of genes using Bio ++ (Guéguen et al. 2013). Ninety-five percent confidence intervals were determined by bootstrapping genes (1000 replicates).

Estimation of the Strength of Positive and Purifying Selection

The contribution of positive selection to the process of amino-acid substitution was estimated using the method of Eyre-Walker and Keightley (2009) as implemented in Galtier (2016). This method, which elaborates on the McDonald–Kreitman test (McDonald and Kreitman 1991), models the distribution of the fitness effect (DFE) of deleterious non-synonymous mutations as a negative Gamma distribution. The model is fitted to the synonymous and non-synonymous site frequency spectra (SFS) and the expected dN/dS under near-neutrality is deduced. Here we used folded SFS to avoid polarization issues. The difference between observed and expected dN/dS provides an estimate of the proportion of adaptive non-synonymous substitutions, α. The per mutation rate of adaptive and non-adaptive amino-acid substitution, respectively ωa = α(dN/dS) and ωna = (1-α)(dN/dS), were also computed.

Positive selection analysis was conducted separately in M. jurtina and P. tithonus, using each species as outgroup to each other. The number of Z-linked coding sequences (genes) for which both polymorphism and divergence data was available was insufficient for a proper estimation of α. Therefore, we used distinct sets of genes for the polymorphism and divergence analyses: synonymous and non-synonymous site frequency spectra (SFS) were built based on 151 Z-linked genes in M. jurtina and 144 Z-linked genes in P. tithonus, whereas dN and dS were calculated based on 90 Z-linked genes for which a pair of orthologs between the two species was available. This means assuming a common DFE across distinct sets of genes. The high similarity of the DFEs between Z and autosomes, as well as the remarkable similarity of the estimated DFEs in mitochondrial data sets from different species observed by James et al. (2016) comfort us in this assumption. A similar procedure was followed for autosomal genes, in which SFS were built based on 5,636 and 5,394 genes in M. jurtina and P. tithonus, respectively, and divergence was estimated from 5,212 genes. Ninety-five percent confidence intervals were determined by bootstrapping genes (1000 repetitions).

When we analyzed Z-linked genes depending on sex-biased expression, we did not have enough genes in each category to use the above-described method. We instead computed a modified version of the Direction of Selection (DoS) statistic (Stoletzki and Eyre-Walker 2011). DoS is defined as the difference between the proportion of fixed differences that are non-synonymous and the proportion of polymorphisms that are non-synonymous [here we used DoS = dN/(dN + dS) − πn/(πn + πs)]. A positive DoS indicates adaptive evolution, DoS = 0 indicates neutral evolution, and a negative DoS indicates segregating slightly deleterious mutations (Stoletzki and Eyre-Walker 2011). Here DoS values does not go with confidence interval to see if the values differ between the three categories because the DoS is computed from two different sets of genes, one for the divergence data and one for the polymorphism data.

Divergence and Polymorphism Comparison between Sex-Biased Expression Genes

If differences in coding sequence evolution between Z and autosomes were primarily driven by an increased efficacy of selection on recessive mutations in females, we would expect a stronger purifying selection in Z-linked genes that are predominantly expressed in females, compared with Z-linked genes predominantly expressed in males. To test this hypothesis, we defined three categories of Z-linked genes: female-biased (i.e., with a higher level of expression in females than in males), unbiased (expressed at approximately the same level in both sexes) and male-biased (with a higher level of expression in males than in females). To define sex-biased genes in the two species, we estimated the expression level of each coding sequence in each individual using the “idxstats” and “depth” tools of the SAMTOOLS library (Li et al. 2009). We computed the « RPKM » (Reads Per Kilobases Per Million) as follow: RPKM = Nc*109/Ntot*Lc, where Nc is the number of reads mapped onto the focal coding sequence, Ntot is the total number of reads mapped of the focal individual, and Lc is the length of the focal coding sequence in base pair (Mortazavi et al. 2008). We calculated for each gene the mean RPKM in females, RPKMf, and the mean RPKM in males, RPKMm.

For the comparison of πn/πs ratios between sex-biased genes, genes for which RPKMf/RPKMm > 1.5 were called female-biased, genes for which RPKMf/RPKMm < 0.66 were called male-biased, and the other ones were called unbiased genes (supplementary tables S2 and S3, Supplementary Materials online). We also compared dN/dS ratios and DoS measurements between Z-linked sex-specific genes, but number of male-biased genes, as defined above, was not sufficient, so we rather sorted genes according to the difference in expression level between males and females and created three bins of equal sizes (30 genes each) (supplementary table S4, Supplementary Materials online). Ninety-five percentage confidence intervals were determined by bootstrapping genes (1000 replicates) in each category for dN/dS as well as polymorphism comparisons. We conducted the same analyses on autosomal genes as a control.

Expression Level Influence on the dN/dS and πn/πs Ratios

We tested if a link between gene expression level and dN/dS or πn/πs ratios existed in our dataset and if it could explain the difference between Z and autosomes. To do that, we established two linear models using R (v 3.1.2, 2014): where dNij is the dN of the ith coding sequence that is linked to chromosome type j (i.e., Z or autosome). For the polymorphism, we added a species specific effect. RPKMi is the mean RPKM of coding sequence i across all individuals. We excluded the coding sequences with no polymorphism or substitution. We ascertained normality, homoscedasticity and independence of the variables by plotting observed versus predicted values.

Results

Transcriptome Assembly, Genotyping, and Expression Level

RNA-seq data were generated in ten individuals per species (supplementary table S1, Supplementary Materials online). Details of the successive sorting steps are presented in supplementary figure S1, Supplementary Materials online. In M. jurtina, transcriptome assembly yielded 145,564 contigs, based on which 36,864 ORFs were predicted. Among these ORFs we recovered orthologous sequences to 11,768 genes from M. cinxia, i.e., a large fraction of the reference genome. Of these, 7,647 corresponded to autosomal genes in M. cinxia, 378 to Z-linked genes, and 3,743 were unassigned and not considered further. These ORFs were genotyped in the ten individuals. Positions (codon sites) at which less than six individuals were sufficiently covered to be accurately genotyped were discarded, leading to the removal of 119 Z-linked and 2011 autosomal ORFs. We conservatively removed 108 predicted Z-linked ORFs that exhibited at least one polymorphic site in at least one female (females are expected to be haploid for Z-linked genes). These can reflect genotyping errors, or genes that have been translocated from the Z to an autosome since the divergence with M. cinxia (Ahola et al. 2014), or genes located in putative pseudo-autosomal regions. Including these genes of uncertain assignment to the Z-linked data set yielded results qualitatively similar to our main analysis (supplementary table S5, Supplementary Materials online). About 151 predicted Z-linked and 5,636 autosomal ORFs were finally selected for polymorphism analysis. Similarly, in P. tithonus we obtained 110,120 contigs, 32,959 ORFs, of which 10,780 had a predicted ortholog in M. cinxia. 353 genes were predicted to be Z-linked, of which 144 were selected for polymorphism analysis, together with 5,394 autosomal genes. For the divergence analysis, we selected and aligned 5,212 autosomal and 90 Z-linked coding sequences that were predicted to be orthologous between M. jurtina and P. tithonus.

In our dataset, the mean level of expression of the Z chromosome was 23% and 41% lower than the autosomal mean expression level in M. jurtina and P. tithonus, respectively (supplementary fig. S2, Supplementary Materials online). In M. jurtina, the mean expression level of Z-linked genes was very similar in males and females (female Z-linked genes expression level was 96% of the male Z-linked genes). In P. tithonus, the mean expression level of Z-linked genes was, on average, 14% higher in males than females.

Z-linked vs. Autosomal Divergence

dN/dS was computed in pairwise alignment (table 1). To appreciate the significance of the difference between Z and autosomes, we sampled without replacement 90 autosomal pairwise alignments (1000 replicates), thus matching the number of available Z-linked pairwise alignments (fig. 1). In contrast to what has been observed in the other Z/W systems studied so far, we did not detect any fast-Z effect, the mean dN/dS ratio of Z-linked genes being even slightly lower than the autosomal one.

Table 1

dN, dS, and dN/dS Ratio Obtained Using Pairwise Alignments for Z-Linked and Autosomal Genes

	#cds	Mean length	dN	dS	dN/dS
Z-linked	90	726	0.025 [0.019; 0.031]	0.31 [0.26; 0.36]	0.082 [0.065; 0.10]
Autosomal	5212	922	0.025 [0.024; 0.026]	0.26 [0.26; 0.27]	0.094 [0.090; 0.097]

Note—Intervals represent 95% confidence intervals obtained by bootstrapping genes (1000 replicates).

. 1.—

Distribution of the dN/dS ratio obtained by resampling without replacement of 90 autosomal pairwise alignments (1000 replicates). dN/dS ratio of the Z-linked pairwise alignments is indicated in red.

Besides, we observe a higher dS on the Z chromosome relative to autosomes, which is consistent with the existence of a male-biased mutation rate (Miyata et al. 1987).

Using a multiple regression analysis, we found that neither gene expression level nor chromosome type (i.e., Z versus autosome) had a significant effect on dN (P = 0.358 and P = 0.285 for expression level and chromosome type, respectively; supplementary table S6, Supplementary Materials online).

Z-linked vs. Autosome Polymorphism

We compared levels of non-synonymous (πn) and synonymous (πs) polymorphism and πn/πs ratio between Z and autosomes (table 2). The Z chromosome exhibited a lower πs than the autosomes in both M. jurtina and P. tithonus. Using πsZ/πsA ratios, we estimated that the NeZ/NeA ratio is below 0.6 in both species, indicating quite an important difference in effective population size between Z and autosomes. The average πn/πs ratio of Z-linked genes was slightly higher than in autosomal genes. To appreciate the significance of the difference between Z and autosomes, we sampled without replacement 151 (in M. jurtina) and 144 (in P. tithonus) autosomal genes (1000 replicates), thus matching the number of available Z-linked genes. We found that the difference in πn/πs between Z and autosomes is not statistically significant, the observed value in Z-linked genes being well within the autosomal distribution (fig. 2). There are no indications that the Z chromosome experiences a reduced efficacy of purifying selection despite its low Ne relative to autosomes (table 2)—everything else being equal; we would expect a higher πn/πs in Z than autosomes due to increased intensity of genetic drift in the former.

Table 2

πn, πs, and πn/πs Ratios Values Obtained from Samples of 151 and 5,336 Genes in M. jurtina and 144 and 5,396 Genes in P. tithonus for Respectively Z-Linked and Autosomal Genes

		M. jurtina	P. tithonus
	πn	0.0009 [0.0007; 0.0011]	0.0008 [0.0005; 0.0010]
Z-linked	πs	0.010 [0.0088; 0.012]	0.006 [0.005; 0.007]
	πn/πs	0.086 [0.068; 0.108]	0.131 [0.096; 0.171]
	πn	0.0022 [0.0022; 0.0023]	0.0012 [0.0012: 0.0013]
Autosomal	πs	0.031 [0.030; 0.031]	0.0096 [0.0093; 0.098]
	πn/πs	0.073 [0.071; 0.078]	0.126 [0.121; 0.131]
	πsZ/πsA	0.334	0.599

Note—Intervals represent 95% confidence intervals obtained by bootstrapping genes (1000 replicates).

. 2.—

Distribution of the πn/πs ratio obtained by resampling without replacement of 151 and 144 autosomal genes (1000 replicates). πn/πs ratio of the Z-linked genes are indicated in red. (A) M. jurtina and (B) P. tithonus.

Using a multiple regression approach, we found that πn is significantly, negatively correlated to expression level both among autosomal genes (P < 2e-16) and Z-linked genes (P = 0.0016, supplementary table S7, Supplementary Materials online). This result is typically interpreted as reflecting an increased strength of purifying selection on highly expressed genes (Drummond et al. 2005). Z-linked genes have, on average, a lower expression level than autosomes. This is another reason, in addition to their reduced Ne, why a significantly higher πn/πs ratio on Z than autosomes was expected. The lack of a Z-chromosome effect on πn/πs despite reduced expression and smaller Ne suggests that purifying selection is more efficient on the Z chromosome than on autosomes.

Purifying Selection and Sex-Specific Gene Expression

Figure 3 shows that, both in M. jurtina and P. tithonus, the πn/πs ratio is higher in male-biased (60 and 51 genes in M. jurtina and P. tithonus, respectively) than in female-biased genes (26 and 42 genes in M. jurtina and P. tithonus, respectively), unbiased genes (65 and 61 genes in M. jurtina and P. tithonus, respectively) being intermediate. The pattern is clearer in M. jurtina than in P. tithonus, where CI are somewhat overlapping. This result is fully consistent with the existence of an effect of hemizygosity on the efficacy of purifying selection against recessive, deleterious mutations. Such a pattern was not detectable in autosomes (supplementary table S8, Supplementary Materials online). Interestingly, πn/πs in Z-linked male-biased genes was not only higher than the Z chromosome average, but also higher than the autosomal average. This might reflect the increased effect of genetic drift in the Z relative to autosomes, which, promotes the segregation of slightly deleterious, recessive alleles in Z-linked male-biased genes because they are not expressed in the heterogametic sex.

. 3.—

πn/πs ratios of Z-linked genes of the following categories: female-biased expression, unbiased expression and male-biased expression. Error bars represent ninety-five percent confidence intervals obtained by bootstrapping genes (1000 replicates). (A) M. jurtina and (B): P. tithonus.

Z vs. Autosomal Rate of Adaptive Substitution

We assessed the prevalence of adaptive evolution in Z-linked vs. autosomal genes. Following the method of Eyre-Walker and Keightley (2009), we computed the proportion of adaptive non-synonymous substitutions α, as well as ωna and ωa, the per mutation rates of non-adaptive and adaptive substitution, respectively. We sampled without replacement 151 (M. jurtina) and 144 (P. tithonus) autosomal genes (1000 repetitions) to match the same number of genes as for the sex chromosome to establish a SFS, this way generating the expected distribution of α, ωa and ωna in Z-linked genes under the hypothesis that they follow the same process as autosomal genes (fig. 4).

. 4.—

Distribution of α, ωa, and ωna obtained by resampling without replacement 90 autosomal genes (1000 replicates). α, ωa, and ωna values of the Z-linked genes (obtained with 90 genes for divergence and 151 and 144 genes for polymorphism for M. jurtina and P. tithonus, respectively) are indicated in red. (a) α in M. jurtina, (b) ωa in M. jurtina, (c) ωna in M. jurtina, (d) α in P. tithonus, (e) ωa in P. tithonus, and (f) ωna in P. tithonus.

Both species showed a similar α between Z and autosomes and a slightly lower ωa on the Z relative to autosomes (table 3 and fig. 4). The 95% confidence intervals were large for Z linked genes and none of the observed differences between Z and autosomes were significant. This suggests that hemizygosity has no strong effect on the rate of adaptive substitution in Satyrinae butterflies, as the prevalence of positive selection in all Z-linked genes taken together was not increased relative to autosomes. In P. tithonus, we observed a slightly lower ωna on the Z chromosome, and in M. jurtina, ωna was similar between Z and autosomes, confirming the pattern of πn/πs variation between chromosome types.

Table 3

Rate of Adaptation α and Adaptive (ωa) and Non-Adaptive (ωna) Substitution Rates

	Z-linked			Autosomal
Species	α	ωa	ωna	α	ωa	ωna
M. jurtina	0.72 [0.34;1]	0.059 [0.027;0.094]	0.023 [0.00;0.052]	0.78 [0.74;0.81]	0.072 [0.068;0.077]	0.021 [0.018;0.023]
P. tithonus	0.66 [0.24;1]	0.054 [0.019;0.094]	0.028 [0.00;0.064]	0.63 [0.59;0.68]	0.059 [0.053;0.065]	0.035 [0.030;0.038]

Note—α, ωa, and ωna were computed according to the method of Eyre-Walker and Keightley (2009), using all the available genes with polymorphism data and all the available genes with divergence data. Intervals represent 95% confidence intervals obtained by bootstrapping genes (1000 replicates).

Positive Selection and Sex-Specific Gene Expression

We aimed at testing whether adaptive evolution was more prevalent in female-biased and unbiased Z-linked genes than in male-biased Z-linked genes. We had too few genes in each of the three categories to estimate the α and ωa statistics with sufficient accuracy, so we only compared dN/dS ratios and the DoS statistic.dN/dS ratios in Z-linked female-biased genes (M. jurtina: dN/dS = 0.11 ± 0.03, P. tithonus: dN/dS = 0.10 ± 0.04) were always higher than in the unbiased (M. jurtina: dN/dS = 0.07 ± 0.03, P. tithonus: dN/dS = 0.08 ± 0.02) and male-biased (M. jurtina: dN/dS = 0.07 ± 0.02, P. tithonus: dN/dS = 0.07 ± 0.02) categories of Z-linked genes, which might indicate an increased rate of adaptive substitutions associated to hemizygosity, especially knowing that no significant difference between categories of gene expression was detected in autosomes (supplementary table S9, Supplementary Materials online). Nevertheless, confidence intervals were overlapping, and no significant difference between unbiased and male-biased Z-linked genes was detected.

Contrasting polymorphism and divergence patterns, we obtained positive DoS values for Z-linked female-biased genes, which is indicative of the presence of adaptive substitutions (0.060 and 0.0022 for M. jurtina and P. tithonus, respectively). For both the unbiased and male-biased Z-linked gene categories, DoS values were negative, indicating a prevalent effect of purifying selection. The DoS statistic was closer to zero in Z-linked unbiased genes (M. jurtina: −0.0026; P. tithonus: −0.029) than in male-biased genes (M. jurtina: −0.030; P. tithonus: −0.062), as expected under the hypothesis of positive selection being enhanced by hemizygosity.

Discussion

Faster rates of coding sequence evolution on the Z chromosome relative to the autosomes have been observed across a wide range of species (Dalloul et al. 2010; Mank et al. 2010; Ellegren et al. 2012; Sackton et al. 2014; Wang et al. 2014; Wright et al. 2015), but different reasons have been given to explain the phenomenon. Studies in birds point to genetic drift as the cause of the fixation of a substantial proportion of slightly deleterious mutations, as the effective population size of the Z chromosome is expected to be particularly low relative to the effective population size of autosomes (Mank et al. 2010; Wright et al. 2015). However, a study on the silk moth Bombyx mori revealed a fast-Z effect that seems to be due to enhanced positive selection on the Z (Sackton et al. 2014).

Here using another Z/W taxa, we did not find evidence for a fast-Z evolution, and our estimated α was similar between Z and autosomes in the two species. However, we did find evidence for an increased efficacy of purifying selection of the Z relative to autosomes: in spite of a reduced effective population size and expression level, Z-linked genes have a πn/πs ratio and a per mutation rate of non-adaptive substitutions (ωna) that are similar to autosomes. We link this phenomenon to the effect of hemizygosity because female-biased Z-linked genes experience a more efficient purge of slightly deleterious mutations and a slightly higher rate of adaptive evolution than male-biased genes.

Does Sex-Specific Expression Predict Sex-Specific Fitness Effects of Mutations?

By comparing divergence and diversity of Z-linked genes according to sex-specific expression, we intended to investigate the influence of hemizygosity on the strength of selection. Our reasoning is based on the assumption that gene expression level is a predictor of the average intensity of the fitness effect of mutations. We therefore assumed that Z-linked genes expressed at a higher level in males than in females are mostly submitted to selection in males, but weakly or not selected in females, thus escaping the effect of hemizygosity—under this hypothesis, recessive mutations affecting a male-biased Z-linked gene would not affect the phenotype at heterozygous state, whereas recessive mutations affecting a female-biased gene would. The same rationale was used to investigate the effect of hemizygosity in mouse (Kousathanas et al. 2014), silk moth (Sackton et al. 2014), fruit flies (Avila et al. 2015), and birds (Mank et al. 2010).

A link between expression level and the fitness effect of mutations is suggested by the negative correlation between dN/dS and expression level that is found in a wide range of organisms (Drummond and Wilke 2008). In our dataset, we did not observe such correlation; however we did observe a negative correlation between gene expression level and πn (supplementary tables S5 and S6, Supplementary Materials online). This tends to indicate that highly expressed genes are more constrained, justifying the use of expression level as a predictor of the intensity of fitness effects of mutations when comparing diversity among classes of sex-biased expression. Z-linked genes with a higher expression in male than in female are not necessarily genes with a male-specific function: these might, alternatively, correspond to genes that are not dosage-compensated, and consequently, to genes that are not dosage-sensible. Such dosage-insensitive genes are likely to be submitted to a relaxed purifying selection, which can lead to an elevation of their πn/πs ratio irrespective of hemizygosity. To rule out this possibility, we estimated the πn/πs ratio of Z-linked genes with more than a twofold male to female level of expression, for which the difference of expression between male and female exceed the difference in ploidy level, and is likely to reflect a sex-biased function. With this restricted set of Z-linked male-biased genes, we obtained πn/πs ratios of 0.072 in M. jurtina and 0.16 in P. tithonus, still higher than the female-biased ratios (respectively 0.044 and 0.10, fig. 3, supplementary table S9, Supplementary Materials online). We conclude that our report of a lower πn/πs ratio in female-biased than in male-biased genes can safely be interpreted as reflecting the effect of hemizygosity.

Gene Expression Level, Gene Content and Recombination: Three Potential Biases

When comparing evolutionary rates between different categories of genes, several biases can potentially arise. Expression level often differs between sex chromosomes and autosomes in the heterogametic sex. There is a lack of consensus in regards to complete dosage compensation in Lepidoptera species (Zha et al. 2009; Harrison et al. 2012; Walters et al. 2015). Here we show that in M. jurtina and P. tithonus, female and male expression levels are roughly similar in Z-linked genes, thereby limiting the potential bias due to differences in gene dosage on the Z chromosome pattern of diversity. Our conclusions are anyway conservative with respect to this potential bias because we report similar πn/πs ratio in Z-linked and autosomal genes (fig. 2 and table 2), despite a lower average gene expression level in the Z chromosome. A recent study in mammals showed that two genomic features, namely GC-content and gene expression level, are sufficient to explain the higher dN/dS of X-linked genes compared with autosomes (Nguyen et al. 2015). The absence of a detectable effect of hemizygosity on dN/dS in mammals could be explained by a more efficient purging of recessive deleterious mutations on the X, which tends to reduce dN/dS, thus offsetting the increased fixation rate of recessive beneficial mutations, similar to the satyrine situation. This hypothesis could not be explicitly tested in mammals because of the relatively low level of within-species polymorphism in this group.

Another potential bias is a difference in gene content between Z and autosomes, which could lead to differences in the adaptive mutation rate between chromosome types. We do not detect any obvious difference between the estimated DFEs of Z and autosomes (supplementary figs. S3 and S4, Supplementary Materials online), which does not suggest that Z-linked genes are more prone to adaptation than autosomal genes. Moreover, enrichment tests performed in birds and primates revealed no significantly enriched gene ontology terms for Z(X)-linked genes relative to autosomes (Hvilsom et al. 2014; Wright et al. 2015). Nevertheless, differences in gene content between Z and autosomes remains largely unknown in satyrine Lepidoptera, so we cannot totally exclude any influence of gene content on our results.

Finally, when comparing Z-linked versus autosomal genes, one should also consider a potential difference in terms of recombination rate. Indeed, it has been shown that in various families of Lepidoptera females lack recombination (Suomalainen et al. 1973; Turner and Sheppard 1975; Traut 1977; Fisk 1989). We can assume that it is also the case in Satyrinae, and consequently, that the population-effective recombination rate for a given rate of recombination r in females between two loci is r/2 for autosomes and 2r/3 for Z. Assuming that background selection is at work, and everything else being equal, one could thus expect πsZ > πsA, and more efficient purifying selection on the Z due to reduce linkage (Charlesworth 2012), perhaps confounding the effects of hemizygosity. Nevertheless, our data do not meet the prediction of a higher neutral diversity on Z than on autosomes, and in the absence of a direct estimation of recombination rate in Satyrines this remains speculative.

Lack of Fast-Z Effect in Satyrine Butterflies: Where Are the Adaptive Substitutions on the Z Chromosome?

In this study, we obtained results that set the Satyrinae apart from the majority of other species in which sex chromosome evolution has been investigated so far: in spite of the evidence for an effect of hemizygosity on purifying selection, we did not observe an increase in ωa on the Z relative to autosomes (table 3), which, combined with the slight decrease in ωna on the Z, led to a slight slow-Z effect. Hemizygosity is expected to promote the fixation of recessive adaptive mutations, and to facilitate the purge of recessive deleterious mutations. Here we show that the latter effect can be stronger than the former. An effect of hemizygosity on the adaptive rate in Satyrinae is indeed suggested by the higher dN/dS and DoS we report in female-biased than in male-biased genes, but the impact of deleterious mutations dominates, so that the net effect is a slow-down, not an acceleration, of molecular evolution on the Z.

Therefore, why would satyrine butterflies behave differently from the other species in which molecular evolution of sex-chromosomes has been examined (Hvilsom et al. 2014; Kousathanas et al. 2014; Avila et al. 2015; Wright et al. 2015)? The rate of evolution of Z-linked genes is determined by the distribution of selection coefficients of mutations occurring on the Z, the distribution of dominance coefficients, and the relative proportion of male-biased vs. female-biased genes. In principle, any peculiarity of Satyrinae regarding one of these parameters, of which we have no empirical measurement, could contribute to explaining the lack of a fast-Z effect in this group. Additionally, if adaptation uses standing genetic variation, if mutations are recurrent, or if individual bouts of adaptation are restricted to a small amount of genes, hemizygosity is not expected to influence the adaptive substitution rate of the Z chromosome (Meisel and Connallon 2013). Any of these situations might apply in satyrines more often than in other groups, for some yet undetermined reason.

Alternatively, it could be that our Z-linked genes sample, based on RNAseq data, is biased towards genes receiving less adaptive mutations than in other studies. Transcriptome-based analysis of Z-linked genes in birds (Wright et al. 2015) did reveal a fast-Z effect, but failed to detect any conspicuous influence of adaptive evolution, similar to our current analysis. Nevertheless, we report for both autosomal and Z-linked genes a high proportion of adaptive substitution (α) (between 0.63 and 0.78), so it is not likely that we missed the fraction of fast adapting genes. There might also be, finally, a publication bias in favor of the fast-Z, especially in the pre-Next Generation Sequencing era, when data sets were limited in size and failure to detect a fast-Z effect could be attributed to lack of power. The analysis of sex-linked genes evolution in additional taxa appears required to settle this issue and confirm, or not, the generality of the fast-Z effect.

The Impact of Genetic Drift on the Sex Chromosome in Female Heterogametic Taxa

The use of the ratio πsZ/πsA as a proxy of NeZ/NeA is a reasonable approximation assuming that there is no sex mutation bias. Here we report a slightly higher dS on the Z chromosome than on autosomes, suggestive of a male-biased mutation rate. However, the difference in dS between the two compartments is small in absolute terms, so we considered that the πsZ/πsA is still representative of NeZ/NeA . This is in line with the results obtained in Bombyx (Sackton et al. 2014) and in birds (Axelsson 2004), where Z-linked genes evolve only marginally faster than autosomes.

Theoretically, the ratio NeZ/NeA is expected to vary between species due to differences in intensity of sexual selection toward males, which reduces the number of reproductive males in the population and biases the efficient sex ratio (Wright et al. 2015). Our report of NeZ/NeA ratios that are substantially lower than 0.75 in M. jurtina and P. tithonus (table 2) are consistent with Bateman's principle, which states that the variation in reproductive success is greater among males than among females (Oberhauser 1988). These ratios are comparable to the ones reported in birds (0.29 < NeZ/NeA < 0.46, Wright et al. 2015), and according to Vicoso and Charlesworth (2009), such ratios should lead to an elevation of the fixation rate of deleterious mutations on the Z chromosome for all ranges of dominance.

Nevertheless, we do not detect a strong effect of genetic drift on the Z, as we observe neither an increase of the dN/dS and πn/πs ratios nor a significantly higher ωna on the Z relative to autosomes (fig. 2, tables 1 and 3). Only in male-biased Z-linked genes, which are under limited influence of hemizygosity, did we detect a reduced efficacy of purifying selection. The difference in effective population size between Z and autosomes does not seem to be sufficient to lead to a strong decrease in the efficacy of purifying selection on the Z compare to autosomes in M. jurtina and P. tithonus. Overall Ne may be one of the reasons explaining the different patterns observed between birds and Lepidoptera (Sackton et al. 2014). A low NeZ/NeA ratio, may thus have weaker consequences on drift rate at sex-linked genes in large than in small populations.

Conclusions

We compared coding sequence evolutionary rates between the hemizygous Z chromosome and autosomes in two closely related species of butterflies. Combining diversity, divergence and expression data, we assessed the influence of hemizygosity in shaping Z-linked genes evolution. Our results indicate that due to hemizygosity, purifying selection is more effective on the Z chromosome relative to autosomes, preventing slightly deleterious mutations to reach fixation. Adaptive substitutions seem to be rare enough not to enhance the adaptive fixation rate in the Z relative to autosomes. Those two results explain why we do not detect a faster Z evolution in M. jurtina and P. tithonus, unlike what has been previously found in other systems. We suggest that the effect of hemizygosity on the fate of recessive deleterious mutations, which has been largely neglected until now, should be taken into account when interpreting patterns of molecular evolution in sex chromosomes vs. autosomes.

Data Accessibility

Illumina raw reads are deposited under the project PRJNA326910 in the SRA database.

Supplementary Material

Supplementary tables S1–S9 and figures S1–S4 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).