Patterns of Z chromosome divergence among Heliconius species highlight the importance of historical demography

Abstract Sex chromosomes are disproportionately involved in reproductive isolation and adaptation. In support of such a “large‐X” effect, genome scans between recently diverged populations and species pairs often identify distinct patterns of divergence on the sex chromosome compared to autosomes. When measures of divergence between populations are higher on the sex chromosome compared to autosomes, such patterns could be interpreted as evidence for faster divergence on the sex chromosome, that is “faster‐X”, barriers to gene flow on the sex chromosome. However, demographic changes can strongly skew divergence estimates and are not always taken into consideration. We used 224 whole‐genome sequences representing 36 populations from two Heliconius butterfly clades (H. erato and H. melpomene) to explore patterns of Z chromosome divergence. We show that increased divergence compared to equilibrium expectations can in many cases be explained by demographic change. Among Heliconius erato populations, for instance, population size increase in the ancestral population can explain increased absolute divergence measures on the Z chromosome compared to the autosomes, as a result of increased ancestral Z chromosome genetic diversity. Nonetheless, we do identify increased divergence on the Z chromosome relative to the autosomes in parapatric or sympatric species comparisons that imply postzygotic reproductive barriers. Using simulations, we show that this is consistent with reduced gene flow on the Z chromosome, perhaps due to greater accumulation of incompatibilities. Our work demonstrates the importance of taking demography into account to interpret patterns of divergence on the Z chromosome, but nonetheless provides evidence to support the Z chromosome as a strong barrier to gene flow in incipient Heliconius butterfly species.


| INTRODUCTION
Comparisons between genomes of diverging populations or species have revealed elevated differentiation on the sex chromosomes in several animals, such as flycatchers (Ellegren et al., 2012), crows (Poelstra et al., 2014), Darwin's finches (Lamichhaney et al., 2015), ducks (Lavretsky et al., 2015) and Heliconius butterflies (Kronforst et al., 2013;Martin et al., 2013;Van Belleghem et al., 2017). These patterns of elevated sex chromosome divergence are sometimes readily interpreted as the result of increased reproductive isolation and reduced admixture on the sex chromosomes and, thus, ascribed to a large-X effect (Box 1). However, it remains unresolved whether such elevated sex-linked divergence actually results from more rapid accumulation of isolating barriers on the sex chromosome or could be explained by differences in effective population size between the sex chromosomes and the autosomes (Meisel & Connallon, 2013;Pool & Nielsen, 2007;Wolf & Ellegren, 2017).
When comparing divergence between genomic regions, such as sex chromosomes versus autosomes, measures of population divergence are influenced by within-population diversity (Charlesworth, 1998;Cruickshank & Hahn, 2014) (Box 2). This is explicitly the case for relative measures such as F ST , but also less directly for absolute measures of divergence such as d XY . For absolute divergence measures, this is because the genetic divergence between two alleles sampled from two species includes both divergence accumulated postspeciation, but also diversity already present in the ancestral population before the split. The latter is strongly dependent on effective population size. In a population under equilibrium conditions where the two sexes have an identical distribution of offspring number, the X chromosome-effective population size and genetic diversity is expected to be three-quarters that of the autosomes. Deviations from this ratio can result from multiple unique features of the sex chromosomes (Box 1), and population size changes in particular can have strong differential influence on sex chromosome compared to autosomal diversity (Pool & Nielsen, 2007). Previous studies attempted to control for differences in effective population size on the sex chromosome, for instance among recently diverged duck species from Mexico, but such studies generally do not account for population size changes (Lavretsky et al., 2015). In order to interpret both relative and absolute measures of divergence on the sex chromosomes as evidence of a disproportionate contribution to species divergence and/ or reduced admixture, we need to also account for demographic changes that can influence diversity of the sex chromosomes.
Here, we explore diversity and divergence on the Z chromosome relative to the autosomes among populations of the Heliconius erato and Heliconius melpomene butterfly clades, using these different measures. The H. erato and H. melpomene clades diverged about 13 million years ago (Kozak et al., 2015) and represent unpalatable and warningly coloured butterflies that have independently radiated into many divergent geographic races and reproductively isolated species. Within both clades, speciation has been accompanied by shifts in M€ ullerian mimicry , and where populations come into contact, hybrid phenotypes usually have reduced survival rates due to strong frequency-dependent selection against intermediate colour pattern phenotypes (Jiggins, Mcmillan, Neukirchen, Mallet, & Nw, 1996;Mallet & Barton, 1989;Merrill et al., 2012;. Two species, H. himera and H. e. chestertonii, are geographic replacements of H. erato in dry Andean valleys. They are partially reproductively isolated, but individuals of hybrid ancestry make up about 10% of the population in narrow transition zones between forms (McMillan, Jiggins, & Mallet, 1997;Merrill, Chia, & Nadeau, 2014;Muñoz, Salazar, Castaño, Jiggins, & Linares, 2010). Similarly, H. cydno and H. timareta are geographic replacements of each other and both are broadly sympatric with H. melpomene. Here, both species are reproductively isolated from H. melpomene by a combination of pre-and postmating isolation (M erot, Salazar, Merrill, Jiggins, & Joron, 2017). Species integrity does not seem to involve structural variation such as chromosomal inversions (Davey et al., 2017). Instead, reproductive barriers include strong selection against hybrids, mate choice and postzygotic incompatibilities ( Figure 1a). Assortative mating has evolved in both the H. erato and H. melpomene clades (Jiggins, Naisbit, Coe, & Mallet, 2001;McMillan et al., 1997;Merrill et al., 2014;Muñoz et al., 2010). In the H. erato clade, sterility and reciprocal-cross asymmetry of hybrid sterility have been reported in crosses between H. erato and H. e. chestertonii (Muñoz et al., 2010), but hybrid sterility is absent between H. erato and H. himera (McMillan et al., 1997). In the H. melpomene clade, female sterility (Haldane's rule; Box 1) and reciprocal-cross asymmetry of hybrid sterility occur in crosses between H. melpomene and H. cydno (Naisbit, Jiggins, Linares, Salazar, & Mallet, 2002), H. melpomene and H. heurippa (Salazar et al., 2005) and H. melpomene and H. timareta (S anchez et al., 2015), as well as between allopatric H. melpomene populations from French Guiana and those from Panama and Colombia (Jiggins, Linares et al., 2001). In support of a large-X effect, sterility in these crosses (H. melpomene 9 H. cydno, H. melpomene 9 H. heurippa and H. melpomene 9 H. timareta) was found to be Z-linked.
The presence of incipient species pairs with different levels of reproductive isolation allows us to examine the relative rate of autosomal and Z chromosomal evolution and the factors that are likely influencing patterns of divergence. We take advantage of a large genomic data set composed of 224 whole genomes representing 20 populations of the H. erato clade and 16 populations of the H. melpomene clade. We also use simulations to evaluate the effect that demographic changes have on the estimate of relative rates of divergence on the Z versus the autosomes and demonstrate that in many comparisons, demography can explain much of the observed elevated divergence on the Z relative to the autosomes. However, by taking into account geographic distance or autosomal divergence as a proxy for gene flow, we show that there is evidence for increased divergence on the Z chromosome for species pairs with known postzygotic reproductive barriers. These rates of increased divergence likely reflect reduced admixture on the Z chromosome and provide support for the Z chromosome being a greater barrier to gene flow in some incipient Heliconius butterfly species. VAN BELLEGHEM ET AL.

Large-x or z Effect and What can Cause It
Sex chromosomes have been repeatedly shown to have a disproportionate role during speciation (Coyne & Orr, 2004), demonstrated by three widespread intrinsic postmating effects (Johnson & Lachance, 2012;Turelli & Moyle, 2007): (i) Haldane's rule, (ii) reciprocalcross asymmetry of hybrid viability and sterility and (iii) the large-X effect. Haldane's rule states that where only one sex of the hybrids has reduced viability or fertility, that sex is most commonly the heterogametic sex (Haldane, 1922). Asymmetry of hybrid viability and sterility refers to the situation where reciprocal crosses often differ in their incompatibility phenotype (Turelli & Moyle, 2007). Finally, the large-X effect highlights the disproportionate contribution of the sex chromosomes to the heterogametic and asymmetric inviability/sterility in hybrids (Coyne & Orr, 1989).
Haldane's rule can generally be explained by between-locus "Bateson-Dobzhansky-Muller incompatibilities" (BDMs) (Dobzhansky, 1935;Muller, 1942;Orr, 1996), in which divergent alleles at different loci become fixed between populations and cause inappropriate epistatic interactions only when brought together in novel hybrid allele combinations (Coyne & Orr, 2004). If interacting loci include recessive alleles on the sex chromosome, only the heterogametic sex will suffer incompatibilities in the F1 generation (Turelli & Orr, 1995). Additionally, BDMs between autosomal loci and the sex chromosomes can be specific to a particular direction of hybridization due to their uniparental inheritance and, thus, also explain asymmetric reproductive isolation (Turelli & Moyle, 2007). Hence, hemizygous expression of recessive alleles on the sex chromosome has been put forward as a cause for the disproportionate role of the sex chromosomes during speciation (dominance theory) (Turelli & Orr, 1995).
In contrast to the dominance theory, there is, however, a large body of observations and theory that propose alternative or additional explanations that can cause a large-X effect and/or explain Haldane's rule (Presgraves, 2008;Wu & Davis, 1993). These factors include faster male evolution resulting from intense sexual selection among males (Wu & Davis, 1993), meiotic drive (Frank, 1991), dosage compensation (Jablonka & Lamb, 1991) and faster-X evolution (Charlesworth, Coyne, & Barton, 1987;Vicoso & Charlesworth, 2006). Faster male evolution targets genes with male-biased effects, which may be sex-linked or autosomal, in which case it can cause Haldane's rule only when males are the heterogametic sex. By contrast, faster-X evolution does not necessarily depend on selection related to sex and is predicted if adaptive new mutations are on average partially recessive and thus are presented more readily to selection on the hemizygous sex chromosomes (Charlesworth et al., 1987). In Drosophila, faster-X evolution has been studied extensively. Although it is not ubiquitous, there is clear evidence for faster-X divergence and adaptation (Counterman, Ort ız-Barrientos, & Noor, 2004), particularly for X-linked genes expressed in male reproductive tissues (reviewed in Meisel & Connallon, 2013). In Lepidoptera (butterflies and moths), Haldane's rule and the large-X (or Z) effect have been reported for numerous species (Presgraves, 2002;Prowell, 1998;Sperling, 1994). As lepidopteran females are heterogametic ZW, while males are ZZ, the Z is equivalent to the X. However, as females are heterogametic, faster male evolution is insufficient to explain Haldane's, but faster-X evolution remains a viable explanation (Sackton et al., 2014). Moreover, in Lepidoptera, the large-X effect extends beyond intrinsic isolating barriers and there are differences in many traits and behaviours that map disproportionately to the Z chromosome (Prowell, 1998;Sperling, 1994). These observations are consistent with the faster accumulation of differences on the Z chromosome (faster-X evolution).

Factors Affecting Sex/Autosome Diversity Ratios
Apart from population size changes, factors that can result in deviations from the expected three-quarter X/autosome (X/A) diversity ratio, and could thus potentially affect divergence measures, include (i) sex-biased demographic events leading to different effective population sizes of males and females (Charlesworth, 2001), (ii) selective sweeps and background selection differently affecting the sex chromosomes (Charlesworth, 2012) and (iii) differences in mutation rates between sexes or between the sex chromosomes and the autosomes (Johnson & Lachance, 2012;Sayres & Makova, 2011).
First, different population sizes of males and females can influence the X/A diversity ratio because two-thirds of the X chromosome population is transmitted through females. A male-biased population would thus decrease the X/A diversity ratio below threequarters, whereas a female-biased population would increase the ratio. This effect would be opposite in female heterogametic sex systems (ZW).
Second, the hemizygous expression of the sex chromosome could result in both higher purifying selection and more efficient selection of beneficial recessive mutations (~selective sweeps) and result in a decrease in the expected X/A diversity ratio (Charlesworth et al., 1987). Additionally, differences in recombination rates can lead to different extent of loss of variation through linked selection and thus background selection (Charlesworth, 2012). In Lepidoptera, meiosis is commonly achiasmatic (no recombination) in the heterogametic sex (females) (Suomalainen, Cook, & Turner, 1973;Turner & Sheppard, 1975). A reduction in recombination rate on the (Continues)

| Sampling
We used whole-genome resequenced data of a total of 109 butterflies belonging to the Heliconius erato clade and 115 from the sex chromosomes compared to autosomes, which is commonly found in Drosophila (Vicoso & Charlesworth, 2009), should thus not be expected to decrease Z/A diversity ratios through increased background selection in Heliconius. On the other hand, it has been suggested that effective recombination should be higher, and thus background selection lower, for the Z chromosome when recombination is absent in females (Charlesworth, 2012). This is because the Z chromosomes spend two-thirds of their time in recombining males, whereas autosomes only spend half of their time in recombining males.
Third, because the male germ line generally involves more cell divisions and thus opportunities for replication errors, sex-linked genes may have different mutation rates. Because X-linked genes spend only one-third of their time in males and two-thirds of their time in females, the X chromosome may be subjected to a lower mutation rate. Conversely, the Z chromosome spends two-thirds of its time in males and may therefore become enriched in genetic variation compared to the autosomes (Johnson & Lachance, 2012;Sayres & Makova, 2011;Vicoso & Charlesworth, 2006). Such increased mutation rates on the Z chromosome could also increase the rate of divergence between Z chromosomes (Kirkpatrick & Hall, 2004).

FIGURE B1
The effect of population size on the coalescent and measures of diversity and divergence. The branches represent two populations, 1 and 2, that have split at a certain time (grey dashed line). This branching event occurs on two chromosomes that have a different population size, such as the autosomes (grey) and X chromosome (green). The black lines show the coalescent of two alleles in each population. The branches show that the coalescence process is influenced by the split time as well as the population size. Population size affects the nucleotide diversity within each population (p), the total nucleotide diversity (p T ) and absolute divergence d XY , but not d a as indicated by the vertical coloured lines. For d a , the average of the within-population nucleotide diversity (p S ) is used as the estimate of the ancestral nucleotide diversity (p anc ). The influence of population size on F ST can be seen as resulting from a decrease in the denominator (p T ), but not in the numerator (p T and p s change proportionately) BOX 1 (Continued) VAN (Davey et al., 2016) reference genomes, respectively, using BWA v0.7 (Li, 2013). PCR duplicated reads were removed using PICARD v1.138 (http://picard.sourceforge.net) and sorted using SAMTOOLS (Li et al., 2009). Genotypes were called using the genome analysis tool kit (GATK) Haplotypecaller ( Van der Auwera et al., 2013). Individual genomic VCF records (gVCF) were jointly genotyped using GATK's genotype GVCFs. Genotype calls were only considered in downstream analysis if they had a minimum depth (DP) ≥ 10, maximum depth (DP) ≤ 100 (to avoid false SNPs due to mapping in repetitive regions), and for variant calls, a minimum genotype quality (GQ) ≥ 30. The W chromosome has not been identified in Heliconius, but read depth comparisons between Z and autosomes in males and females (e.g., see supplementary material Martin et al., 2013) support the hypothesis

BOX 2 Measures of Divergence Depend on Population Size
The mutational diversity in present-day samples is directly related to population size, structure and age. This is because population size determines the rate of coalescence within and between populations ( Figure B1). This relationship can be seen with F ST , which was developed to measure the normalized difference in allele frequencies between populations (Wright, 1931). The dependence of F ST on population size can be understood by interpreting F ST as the rate of coalescence within populations compared to the overall coalescence rate (Slatkin & Voelm, 1991). Comparing pairs of populations with different effective population sizes will therefore show distinct F ST estimates even when the split time is the same (Charlesworth, 1998). Absolute divergence d XY is the average number of pairwise differences between sequences sampled from two populations (Nei & Li, 1979). The measure d XY is not influenced by changes to within-population diversity that occur after the split, but does depend on diversity that was present at the time the populations split (Gillespie & Langley, 1979). Therefore, population pairs that had a smaller population size at the time they split will, consequently, have smaller d XY estimates.
To compare pairs of populations that had different ancestral population sizes, d a has been proposed, which subtracts an estimate of the diversity in the ancestral population from the absolute divergence measure d XY (Nei & Li, 1979). An approximation of ancestral diversity can be obtained by taking the average of the nucleotide diversity observed in the two present-day populations. Such a correction should result in the "net" nucleotide differences that have accumulated since the time of split.
To show how these different divergence measures perform, we simulated a simplified population split with varying degrees of migration (m) ( Figure B2). As expected, the values F ST , d XY and d a all increase when migration between populations decreases. F ST and d XY are clearly influenced by population size. While for d XY , this results from the variation present at the time of the split, F ST does not show a simple linear relationship with population size. Only d a represents the net accumulation of differences that can be compared between populations of different sizes, such as the sex chromosomes versus autosomes (but see section 3.3 in Results and discussion).

FIGURE B2
Simulated effect of population size differences on divergence measures F ST , d XY and d a . Simulations were performed for two populations that split 4 million generations ago and vary in their degree of migration (m). A lower effective population size, such as for the X chromosome (green) compared to autosomes (grey), results in higher F ST and lower d XY estimates, but has no effect on d a under these assumptions F I G U R E 1 Diversity and sampling of the Heliconius erato and Heliconius melpomene clade. (a) Heliconius erato chestertonii (green) is reproductively isolated from H. erato (pink) by spatial separation (parapatry), mate choice and (asymmetric) reduced hybrid fertility of both sexes (i.e., no Haldane's rule). Heliconius himera (blue) is reproductively isolated from H. erato by spatial separation and mate choice, but hybrids show no reduced fertility. Heliconius cydno (green) and H. timareta (blue) occur sympatrically with H. melpomene (pink) populations, but are both reproductively isolated by strong mate choice and (asymmetric) reduced fertility of F1 hybrids (i.e., Haldane's rule).  that there is no significant mapping of W-linked reads to the Z and the W is, thus, unlikely to interfere with genotyping. The absence of mapping of W-linked reads to the Z is likely due to the degenerate sequence and highly repetitive nature of the W chromosome. The data set contained 31 and 11 female samples (ZW) randomly distributed among the H. erato and H. melpomene clade populations, respectively (Tables S1 and S2). These samples had lower read and, consequentially, lower genotyping coverage for the Z chromosome, but using the stringent filter thresholds, this does not affect variant and nonvariant sites differently and does not affect measures of nucleotide diversity and divergence.

| Population structure and historical demography
To discern population structure among the sampled H. erato and H. melpomene clade individuals, we performed principal component analysis (PCA) using EIGENSTRAT SmartPCA (Price et al., 2006). For this analysis, we only considered autosomal biallelic sites that had coverage in all individuals.
We inferred changes in the historical population size from individual consensus genome sequences using pairwise sequentially Markovian coalescent (PSMC) analyses as implemented in MSMC (Schiffels & Durbin, 2014). This method fits a model of changing population size by estimating the distribution of times to the most recent common ancestor along diploid genomes. When used on single diploid genomes, this method is similar to pairwise sequentially Markovian coalescent (PSMC) analyses (Li & Durbin, 2011). Genotypes were inferred from BWA v0.7 (Li, 2013) mapped reads separately from previous genotyping analysis using SAMTOOLS v0.1.19 (Li et al., 2009) according to authors' suggestions. This involved a minimum mapping (-q) and base (-Q) quality of 20 and adjustment of mapping quality (-C) 50. A mask file was generated for regions of the genome with a minimum coverage depth of 30 and was provided together with heterozygosity calls to the MSMC tool. MSMC was run on heterozygosity calls from all contiguous scaffolds longer than 500 kb, excluding scaffolds on the Z chromosome. We scaled the PSMC estimates using a generation time of 0.25 years and a mutation rate of 2e-9 as estimated for H. melpomene (Keightley et al., 2014).

| Population genomic diversity and divergence statistics
We first estimated diversity within populations as well as divergence between parapatric and sympatric populations in nonoverlapping 50-kb windows along the autosomes and Z chromosome using python scripts and egglib (data presented in Figure 3, 7 and S4-6) (De Mita & Siol, 2012). We only considered windows for which at least 10% of the positions were genotyped for at least 75% of the individuals within each population. For females, haploid was enforced when calculating divergence and diversity statistics. Sex of individuals was inferred from heterozygosity on the Z. F ST was estimated as in Hudson, Slatkin, and Maddison (1992), as with nucleotide diversity in a population (p i ) calculated as and average within-population nucleotide diversity (p S ) calculated as the weighted (w) average of the nucleotide diversity (p i ) within each population l and k, as Total nucleotide diversity (p T ) was calculated as the average number of nucleotide differences per site between two DNA sequences in all possible pairs in the sampled population (Hudson et al., 1992), as Between-population sequence divergence d XY was estimated as the average pairwise difference between sequences sampled from two different populations (Nei & Li, 1979), as The relative measure of divergence d a was calculated by subtracting d XY with an estimate of the nucleotide diversity (p S ) in the ancestral populations (Nei & Li, 1979), Tajima's D was calculated as a measure of deviation from a population evolving neutrally with a constant size, with negative values indicating an excess of rare alleles (~selective sweep or population expansion) and positive values indicating a lack of rare alleles (~balancing selection or population contraction) (Tajima, 1989). To overcome the problem of nonindependence between loci, estimates of the variance in nucleotide diversity (p) and Tajima's D within populations along the genomes were obtained using block-jackknife deletion over 1-Mb intervals along the genome (chosen to be much longer than linkage disequilibrium in Heliconius ) (K€ unsch, 1989 (Mantel, 1967).
Mantel's tests are commonly used to test for correlations between pairwise distance matrices and were performed using the R package VEGAN (Oksanen et al., 2016). Pairwise distances between populations were calculated from the average of the sample coordinates obtained for each population (Table S3, S4).

| Simulations
To compare patterns in our data to expectations, we simulated genealogies in 50-kb sequence windows under certain evolutionary scenarios. The simulations were performed with a population recombination rate (4N e r) of 0.01 using the coalescent simulator msms (Ewing & Hermisson, 2010). Subsequently, from the simulated genealogies, we simulated 50-kb sequences with a mutation rate of 2e-9 a Hasegawa-Kishino-Yano substitution model using seq-gen (Rambaut & Grass, 1997).
In a first set of simulations, we considered one population that that maintain their integrity, despite ample opportunity for hybridization and gene flow (Arias et al., 2008;Jiggins et al., 1996;McMillan et al., 1997

| Z chromosome divergence in Heliconius erato and Heliconius melpomene
To compare rates of divergence between the Z chromosome and autosomes, we calculated three commonly used measures of diver-   (Pool & Nielsen, 2007). To explore this, we performed coalescent simulations of sequences from populations that underwent a single size change in the past, varying the time and magnitude of this event (Figure 4). In these simulated populations, the decrease in nucleotide diversity that follows population contraction occurs much faster than the increase in diversity that follows an expansion of the same magnitude (Figure 4a). This is because an increase in diversity requires mutation accumulation, whereas drift can more rapidly remove variation to reach a new equilibrium.
Additionally, population size changes have proportionately stronger effects on diversity on the Z chromosome compared to the autosomes ( Figure 4b). This results from populations with a smaller effective population size, such as the Z chromosome, converging faster to their new equilibrium after a population size change (Pool & Nielsen, 2007).
The result of population size change differently affecting the Z chromosome is that divergence measures are also differentially affected by population size change on the Z chromosome compared to the autosomes. The Z chromosome to autosome (Z/A) diversity ratio will be larger than expected in populations that experienced a recent expansion and smaller than expected in those that experi- All the simulations were run for timescales relevant to Heliconius divergence and, therefore, demonstrate that a return to equilibrium values is unlikely after a population increase during the history of these species. In particular, a return to equilibrium Z/A diversity ratios after population size increase can be slow and longlasting during the evolutionary history of a population.

| Demography and its influence on Z/A diversity ratios in Heliconius
To explore how population size changes might have affected Z/A diversity ratios and thus Z/A divergence comparisons within Heliconius clades, we used the behaviour of Tajima's D as a way to assess likely population size changes within species. Tajima's D is a population genetic measure commonly used to detect whether a locus is evolving neutrally in a population (Tajima, 1989). . This branching event occurs on two chromosomes that have a different population size, such as the autosomes (grey) and X chromosome (green). The black lines show the coalescent of two alleles in each population. This coalescence process is influenced by the split time as well as the population size (see Box 2). Moreover, red arrows show the disproportionate effect of population size change on diversity and divergence measures on the sex chromosome (or populations with a smaller size). Note that the effect size will depend on the magnitude as well as the timing of population size change (see Figure 4). Exact expectations, including more complex size changes, can be calculated as demonstrated in Pool and Nielsen (2007) [Colour figure can be viewed at wileyonlinelibrary.com] VAN (Figures 2 and 7). Therefore, the patterns among these Heliconius populations suggest that differences in nucleotide diversity as well as differences in the Z/A diversity ratios are likely driven at least in part by population size changes. Given that samples assigned to a population were collected in close proximity, it is unlikely that estimated Tajima

| Sex-linked incompatibilities increase Z/A absolute divergence ratio
Despite the difficulties in directly comparing divergence on sex chromosomes and autosomes, it may be possible to detect enhanced barriers to migration on sex chromosomes (i.e., reduced effective migration) by comparing population pairs with different levels of absolute migration due to physical isolation, but that otherwise share the same common history. This can be achieved by comparing pairs of populations from the same two species that differ in their extent of geographic isolation. Indeed, previous analyses of sympatric and F I G U R E 6 Simulated effect of population size change on Tajima's D values. We simulated a single population that underwent a population size change of magnitude x, ranging from 0.01 to 100, at a certain moment backwards in time (t). Triangles indicate population size contractions, and circles indicate population size increase. Colours and x-axis represent the time at which the population size change occurred in generations (log 10 (t), going from 1000 to 4,000,000 generations ago) with population size (N e ) equal to 3e6 and 2.25e6 for the autosomes and Z chromosome, respectively. The colours and time at which population size change occurred correspond to colours in the simulations in Figure 7. F I G U R E 7 Relation between Tajima's D, average nucleotide diversity on the autosomes (p A ) (upper panels) and the ratio of nucleotide diversity between the Z chromosome and autosomes (p Z /p A ) (lower panels) for populations of the Heliconius erato and Heliconius melpomene clade and simulated data. Points represent average nucleotide diversity measures obtained from autosomes (p A ) and the Z chromosome (p Z ) in 50-kb windows. Grey vertical bars represent 95% confidence intervals estimated from block-jackknifing (note that these are too small to see in the Tajima's D versus p A plots). Schematics in the upper right panel represent the simulated population size changes. We simulated a single population that underwent a population size change of magnitude x, ranging from 0.01 to 100 (right panels), with population size change occurring between 1000 and 4,000,000 generations (t) ago (colours  (A)) for populations of the Heliconius erato and Heliconius melpomene clade and simulated data. Boxplots represent pairwise measures over 500-km bins. Schematics in the upper right panel represent the simulated populations with different rates of migration (m) expressed as a proportion of the effective population size. We considered pairs of populations that split between 1 and 2 million generations ago (colours) and were connected through migration (m) ranging from 0 to 1e-6 and for which migration was reduced with a factor d on the Z chromosome (size of circles  Figure 1).
Next, to account for geographic barriers, we also carried out a similar comparison using absolute divergence on the autosomes instead of geographic distance, which might reflect a more direct relationship with migration. Using the absolute divergence on the autosomes as a proxy for gene flow, we find a pattern of increased Z/A divergence ratios for species pairs with known postzygotic reproductive barriers (Figure 9). Z/A divergence ratios are signifi- showing that the former and not the latter pairs experience hybrid sterility and Haldane's rule (McMillan et al., 1997;Merrill et al., 2012;Muñoz et al., 2010;Naisbit et al., 2002;Salazar et al., 2005;S anchez et al., 2015) and also agrees with the previous observation of reduced shared variation between H. melpomene and both H. timareta and H. cydno on the Z chromosome .
Note that our simulations suggest that to explain the trend observed in our data, there must be a very strong reduction in migration on the Z chromosome relative to the autosomes (~60% or greater).

| Alternative factors affecting Z/A diversity ratios in Heliconius
Factors other than population size change could result in deviations from the expected Z/A diversity ratio (Box 1). In Heliconius, there is no empirical data on sex chromosome-biased mutation rates. While higher mutation rates on the Z could explain increased Z/A diversity ratios and increased rates of divergence (Kirkpatrick & Hall, 2004;Sayres & Makova, 2011;Vicoso & Charlesworth, 2006), it is unlikely that closely related populations would differ in their mutation rate and that this could explain the observed variation in Z/A diversity ratios among the Heliconius populations. Alternatively, in Heliconius, male-biased sex ratios have been reported in the field, which could result in increased Z/A diversity ratios. However, it has been argued that these malebiased sex ratios are most likely explained by differences in behaviour, resulting in male-biased captures rather than effective sex ratio differences (Jiggins, 2017 (Charlesworth, 2001). However, the frequency of remating in polyandrous Heliconius species is estimated to be only 25-30% higher than in monandrous species  and adult mating is likely still prevalent in presumed pupal-mating species (Thurman, Brodie, Evans, & McMillan, 2018). Correspondingly, we did not find any clear difference in Z/A diversity ratios between the pupal-mating H. erato and adult-mating H. melpomene clade populations ( Figure 7). Finally, it has been suggested that effective recombination may be higher for the Z chromosome when recombination is completely absent in females (Charlesworth, 2012). This is because the Z chromosomes spend two-thirds of their time in recombining males, whereas autosomes only spend half of their time in recombining males. This could lead to less of a reduction in diversity on the Z compared to autosomes than the 25% null expectation. However, this cannot explain the correspondence of Tajima's D and the PSMC inferences and the observed Z/A diversity ratios. Similarly, the pattern of increased Z/A divergence that results from reduced admixture on the Z in population comparisons of geographically closely located Heliconius species should not be affected by sex ratio or recombination and mutation biases. Overall, in Heliconius, the largest influence on variation in Z/A diversity ratios is likely to be demographic changes.

| Consequences for other study systems
Extensive genomic sampling is available for a number of other natural systems that have recently diverged, particularly for birds that also have ZW sex chromosomes, such as flycatchers (Ellegren et al., 2012), crows (Poelstra et al., 2014) and Darwin's finches (Lamichhaney et al., 2015). In these systems, increased coalescence rates (~lineage sorting) on the Z and/or W chromosome have been accredited to the smaller effective population sizes of the sex chromosome. However, it remains unclear whether elevated measures of divergence could indicate elevated rates of between species divergence on the sex chromosomes, resulting from increased mutation or reduced admixture.
In the adaptive radiation of Darwin's finches, there is no evidence for Haldane's rule nor for reduced viability of hybrids due to postmating incompatibilities (Grant & Grant, 1992) and the VAN | 3867 maintenance of isolating barriers is best explained as resulting from ecological selection and assortative mating (Grant & Grant, 2008). In crows, the divergence between hooded and carrion crows seems to be solely associated with colour-mediated assortative mating even in the apparent absence of ecological selection (Poelstra et al., 2014;Randler, 2007). The populations of both Darwin's finches and crows can be characterized by distinct demographic histories (Lamichhaney et al., 2015;Vijay et al., 2016). Therefore, in these species, deviations in divergence measures from neutral expectation on the Z chromosome are potentially also explained by demography. In the divergence of pied and collared flycatchers, species recognition and species-specific male plumage traits are Z-linked (Saether et al., 2007) and female hybrids are completely sterile compared to only low levels of reduced fertility in males (Veen et al., 2001). In agreement with the large-X effect and disjunct rates of admixture between the sex chromosomes and autosomes, genome scans have found signals of increased relative divergence on the Z and W chromosomes (Ellegren et al., 2012;Smeds et al., 2015). The demographic history of these populations is, however, characterized by a severe decrease in population size since their divergence (Nadachowska-Brzyska et al., 2013). In particular, for the W chromosome, the reported excessive decrease in diversity and the high values of relative divergence can thus likely be partly explained by demography (Smeds et al., 2015). However, the excess of rare alleles (~negative Tajima's D) on the W chromosome does contrast with these inferred demographic histories and provides support that the reduced diversity and increased F ST measures result from selection (Smeds et al., 2015).

| CONCLUSION
The disproportionate role of sex chromosomes during speciation has been well documented based on genetic analysis. However, it is less clear how this influences patterns of divergence in natural populations. In Heliconius, we find much of the observed increased absolute divergence on the Z chromosome relative to neutral expectation can be explained by population size changes. This cautions against highlighting increased sex chromosome divergence alone as evidence for a disproportionate role in species incompatibilities or as evidence for faster-X evolution. Although relative measures of divergence are most prone to demographic changes, absolute divergence measures can also be strongly influenced by population size changes. Absolute measures do not therefore provide a solution to the problems inherent in using relative measures to compare patterns of divergence across genomes (Cruickshank & Hahn, 2014). Despite these difficulties, we do find patterns consistent with decreased effective migration on the Z for species pairs with known postzygotic reproductive barriers, in agreement with hybrid sterility and inviability being linked to the Z chromosome in F I G U R E 9 Relation between pairwise autosomal d XY (a) and the ratio of d XY between the Z chromosome and autosomes (d XY (Z)/d XY (A)) for populations of the Heliconius erato and Heliconius melpomene clade and simulated data. Boxplots represent pairwise measures over 1.25e-3 d XY bins. Schematics in the upper right panel represent the simulated populations with different rates of migration (m) expressed as a proportion of the effective population size. We considered pairs of populations that split between 1 and 2 million generations ago (colours) and were connected through migration (m) ranging from 0 to 1e-6 and for which migration was reduced with a factor d on the Z chromosome (size of circles). Population size (N e ) was equal to 2e6 and 1.5e6 for the autosomes and Z chromosome, respectively [Colour figure can be viewed at wileyonlinelibrary.com] these cases (Jiggins, Linares, et al., 2001;Naisbit et al., 2002;Salazar et al., 2005;S anchez et al., 2015). Successfully disentangling the influence of a large-X effect and faster-X evolution on relative rates of divergence will require modelling of the demographic history of each population, including changes that may have occurred before the split of the populations. Such modelling would allow us to better contrast (i) expected within-population Z/A diversity ratios with hypotheses of increased mutation rates, selective sweeps, background selection and mating system and (ii) expected between population Z/A divergence ratios with hypotheses of increased mutation rates or adaptive divergence on the Z chromosome. Additionally, our strategy of contrasting d XY (Z)/d XY (A) ratios with geographic distance provides opportunities for testing reduced admixture between sex chromosomes in systems for which tree-based approaches and/or crossing experiments are unfeasible.

DATA ACCESSIBI LITY
Genome assemblies are available on lepbase.org. Sequencing reads are deposited in the Sequence Read Archive (SRA). See Tables S1 and S2 for accession numbers.