Influence of the large‐Z effect during contact between butterfly sister species

Abstract Recently diverged butterfly populations in North America have been found to exhibit high levels of divergence on the Z chromosome relative to autosomes, as measured by fixation index, Fst. The pattern of divergence appears to result from accumulation of incompatible alleles, obstructing introgression on the Z chromosome in hybrids (i.e., the large‐Z effect); however, it is unknown whether this mechanism is sufficient to explain the data. Here, we simulate the effects of hybrid incompatibility on interbreeding butterfly populations using a model in which populations accumulate cross‐incompatible alleles in allopatry prior to contact. We compute statistics for introgression and population divergence during contact between model populations and compare our results to those for 15 pairs of butterfly species interbreeding along a suture zone in central Texas. Time scales for allopatry and contact in the model are scaled to glacial and interglacial periods during which real populations evolved in isolation and contact. We find that the data for butterflies are explained well by an otherwise neutral model under slow fusion conditions. In particular, levels of divergence on the Z chromosome increase when interacting clusters of genes are closely linked, consistent with clusters of functionally related genes in butterfly genomes.

peninsula to the tip of Florida (see Figures S1 and S2). The species sampled by Cong et al. diverged on the order of 1 million years ago  and have, as a result, experienced multiple periods of glacial cooling and interglacial warming. During glacial periods, central Texas was subjected to severe decreases in temperature (Annan & Hargreaves, 2013), which would have caused drastic, if not total isolation of sister species in southeastern and southwestern refugia; during the most recent warming period, sister species migrated into Texas, while major portions of their populations remained in refugial regions, isolated from the suture zone by large distances. To determine the influence of hybrid incompatibility with the Z chromosome during contact, we will at first neglect the effects of isolation by distance and consider a generic model of secondary contact (Geneva et al., 2015;Harris & Nielsen, 2016) in which a population divides, and the resulting sister populations evolve for a period in allopatry while accumulating hybrid incompatibilities and later begin to interbreed. We then compare statistics for introgression and population divergence for gene sequences in our model to those obtained for real populations by Cong et al. (2019).  (Geneva et al., 2015), defined as the fraction of independent sequence windows along a genome with G min ≤ G 0 , where G 0 is a threshold for introgression, and (b) the fixation index, or relative divergence function, F st (Bhatia et al., 2013); I gf measures the fraction of sequence windows where introgression has occurred, while F st measures the degree of genetic difference between populations (details are provided in the Methods section). Multiple genomic samples were collected from each sister population, and separate index values were computed for autosomes and Z chromosomes. The results are shown in Figure 1; data points in this figure describe index values for sister organisms that have been classified as different species in the literature (green), more closely related organisms for which classification is uncertain (yellow), and samples of the same species (red).
When populations are compared through their autosomes (Figure 1a), the data exhibit a continuous pattern across the entire range of index values; however, for the Z chromosome (Figure 1b), the data obtained from samples of the same species (red) are separated from those of closely related species by a gap of "missing" values, suggestive of a sudden transition (Kronforst et al., 2013;Nosil et al., 2017). For different species (green and yellow data points), relative divergence (F st ) values for the Z chromosome are always larger than those for autosomes ( Figure 2). However, the fraction of divergent nucleotide positions in the Z chromosomes of sister species is similar to those for autosomes (see Figure 5a of reference Cong et al., 2019), indicating similar absolute rates of divergence, inconsistent with a "faster-Z" effect (Avila et al., 2014;Meisel & Connollon, 2013). In accord with these results, Cong et al. have argued that the pattern of data in Figure 1 reflects the influence of negative fitness interactions between autosomes and Z chromosomes in hybrids during periods of interbreeding-in other words, a "large-Z" effect (Muirhead & Presgraves, 2016;Van Belleghem et al., 2017).
Our goal in this work is to determine whether this mechanism can explain the data for butterflies-specifically, the gap in Figure 1b, and the large differences, ΔF = F Z − F A , between F st values for autosomes and Z chromosomes shown in Figure 2. To accomplish this, we simulate populations of model butterflies under conditions that scale to those experienced by real populations during glacial isolation and contact. Chromosomes in our model consist of adjacent gene segments of identical length (see Appendix A, Figure A1).
Mutation rates and rates of crossing over within gene segments are scaled to estimates for Drosophila, and Heliconius butterflies, and rates of crossing over between gene segments are varied to reflect the typical separation, or degree of linkage between genes on butterfly chromosomes. Mutations are individually neutral. Hybrid incompatibility in the model occurs as a result of negative fitness interactions between mutant alleles that rise to fixation in different populations during allopatry, similar to the model described by Orr (Orr, 1995). Interactions are pairwise and connect loci in autosomes to loci in the Z chromosome(s). The fitness cost for a pair of interacting loci resembles the "pathway" model described by Lindke and Buerkle (Lindtke & Buerkle, 2015). Depending on the strength of the interactions and the migration rate, the model leads either to fusion or continued divergence during contact. Here, we focus primarily on fusion conditions. For a given set of conditions (i.e., interaction strength, migration rate, time in allopatry, etc.), we conduct multiple simulations in parallel to generate statistical profiles for index values and other quantities of interest. We first show that a purely neutral model, in which hybrid interactions are turned off during contact, is unlikely to explain the data for ΔF (we refer to this situation as the "null model" below). During allopatry, F Z is slightly larger than F A , as expected, due the smaller population size, and hence higher substitution rate for Z chromosomes (Van Belleghem et al., 2017) (populations maintain an approximately 1:1 sex ratio in the model); results for the null model during contact are similar to those for populations evolving in allopatry. However, as hybrid interactions are increased, profiles for ΔF begin to resemble the pattern of data in Figure 2, particularly when the rate of crossing over between genes is small, reflecting closely linked genes on butterfly chromosomes. In this case, which would perhaps correspond to interacting clusters of functionally related genes (Cong et al., 2016;McDonald & Rosbash, 2001), we find that the model can explain the large values of ΔF shown in Figure 2 under realistic conditions. The index of gene flow, I A , agrees (on average) with the data in Figure 1a during both allopatry and contact, and the model also leads to a "statistical" gap in I Z analogous to that Figure 1b

| ME THODS
We simulate model populations in three phases: (a) equilibration of an initial, ancestral population, (b) division of the population and evolution of sister populations in allopatry, and (c) contact between sister populations, subject to hybrid interactions between mutant alleles acquired in allopatry. In each phase, populations evolve by plain Wright-Fisher dynamics with random mating between male and female individuals (Gillespie, 2004). In each generation, mutations occur within genes at a rate per gene per generation. Pairs of male and female individuals are then selected at random according to fitness for mating. Male genomes undergo explicit meiosis, in which chromosomes are duplicated, and the resultant chromatids undergo random crossing over (Veller et al., 2018) with separate rates, r and r ′ , for crossing over within and between gene segments We consider two possible evolutionary scenarios; the scenarios are the same, except for the initial population size and the transition into allopatry: In scenario (i), the initial population is cloned, while in scenario (ii), the initial population divides into equal parts. We first describe scenario (i) and then describe the differences between scenarios (i) and (ii).
In scenario (i), an initial population of size N is equilibrated for Δt E generations. Let g l = g l , g � l denote the allelic state of a diploid locus l in a genome g. All genomes in the initial population begin with g l = 0 uniformly. Mutation events during a simulation act to assign the maternal (g l ) or paternal (g ′ l ) state of a locus to 1. All mutations are individually neutral. After equilibration, loci l that have fixed in the population for the mutant allele type are returned to their initial states, g l = 0. The population is then duplicated, and the resultant sister populations (each of size N) evolve in allopatry for a period Δt A . At the end of this period, loci that have fixed for the mutant allele across both populations are returned to their initial states. Several pairs of loci are then selected to participate in hybrid fitness interactions (see below), and the two populations evolve in contact for a period Δt C subject to fitness costs incurred due to interactions formed by various allele combinations at the selected loci.
In scenario (ii), the entire procedure is the same, except that the initial population has size 2N before dividing into sister populations of size N (in this case, the equilibration period is twice as long). In both scenarios, the W chromosome acts only to determine the sex of an individual. To describe the cost of hybridization, we select a small number of loci in autosomes for which p i,1 ∼ 1 and p i,2 ∼ 0 to interact negatively with loci in the Z chromosome(s) for which p j,1 ∼ 0 and p j,2 ∼ 1 (see Appendix A).
We then repeat this process with the population subscripts interchanged, selecting an equal number of loci in autosomes with p i,2 ∼ 1 and p i,1 ∼ 0 to interact negatively with loci in the Z chromosome(s) for which p j,2 ∼ 0 and p j,1 ∼ 1 (typically, p l, = 0, 1 for selected loci in our simulations, but in general, this will depend on the mutation rate, the amount of time spent in allopatry, and the number of selected loci).
Let f g i , g j denote the log fitness cost for a pair of selected loci, . We define f as follows: If both loci are homozygous for the mu- The model for f is the same as the "pathway" model used by Lindtke and Buerkle (2015) to describe hybrid interactions among autosomal loci, except for the factor f = s∕4 when both loci are heterozygous (in this case, f = 0 in the pathway model).
Genomes in our simulations contain three sets of chromosomes; autosomes carry three genes, and Z chromosomes carry six genes, each of length L loci, as shown in Appendix A, Figure A1. Since interacting loci are required to have p l, ∼ 1 in one of the diverging populations, we are somewhat limited in regard to the number of interacting loci we can select to define w in a given simulation.
Here, we select six pairs of loci in each simulation for which the differences between p l,1 and p l,2 above are largest. Unless otherwise noted, we select loci that connect the first pair of autosomes to the Z chromosome(s).  Halligan and Keightley (2009). Assuming an effective population size for Drosophila of 10 6 individuals Li et al., 1999), and a typical gene length of about 1770 bp (Keightley & Eyre-Walker, 2000), we obtain a scaled mutation rate of N ≃ 0.18 − 1.8 per gene per generation. Here, we adopt a value 2N = 1 in the lower range of these estimates (during publication, we became aware of a significantly larger estimate for Heliconius (Keightley, Pinharanda, et al., 2014)

| RE SULTS
To begin our investigation, we explore the time dependence of the statistics G min and F st for populations evolving in allopatry for comparison with the results of Geneva et al. (Geneva et al., 2015).
To define these objects, let d l = | | | g l − g l | | | denote the difference (Hamming distance) between genomes g and g at (haploid) locus l , and let denote the distance between g and g for a window of loci, [ , + Δ]. Assume that we have sampled a small number of genomes prior to contact is F A ≃ 0.5. The results are similar to those obtained during allopatry in Figure 4, and, accordingly, the null model is unlikely to explain the data in Figure 2; note that I Z is shifted upward from I A by an amount similar to that in Figure 4a.
In Figures 6 and 7

| D ISCUSS I ON
The model seems to capture the data in Figures 1 and 2  However, in these situations, the only distinguishing factor between the dynamics of genes on autosomes and Z chromosomes is population size. Thus, for example, we would expect I Z to approach I A (i.e., which is similar to I Z in Figure 1) if males were more abundant than females in the model.
Another issue is the estimate used for the mutation rate, N . The recent estimate for Heliconius noted above (Keightley, Pinharanda, et al., 2014) is several times larger than the estimate used in our simulations. Larger mutation rates would lead to more rapid divergence of populations in allopatry and different conditions for fusion and continued divergence following contact. However, it is worth recalling that the value of N e used to estimate the length of glacial periods for the model is in the lower range of values for Heliconius. (Van Belleghem et al., 2017). Larger and perhaps more realistic values of N e would lead to shorter glacial periods (i.e., smaller values of in the relation N e = Δ above), which would act to compensate for an increased mutation rate in allopatry. In addition, populations would have less time to interbreed during contact, leading to plots that more closely resemble those for slow fusion in Appendix B.
Finally, it is important to remark that index values computed for a pair of species will depend on where the specimens are collected.
The locations of specimens studied in this work often extend over thousands of kilometers on either side of the suture zone (see, e.g., Figure S1). In these distant regions of the landscape, sister populations evolve in greater isolation and, hence, diverge at a higher rate. side of the suture zone are different, which has probably led to some level of divergent adaptation (e.g., discordant mating cycles (Cong et al., 2016), mate preferences (Kronforst et al., 2013), and environment preferences), limiting the rates of interbreeding and gene flow in some complex way (Edelaar et al., 2008;Flaxman et al., 2014;M'Gonigle et al., 2012). Given the present scale of computing power, and the increasing ease of obtaining genetic information, it would be worthwhile to develop software capable of modeling the properties above for realistic population sizes and genome structures (Haller & Messer, 2019). The present work is the first step toward this goal.

ACK N OWLED G M ENTS
It is a pleasure to thank Jing Zhang and two anonymous reviewers for helpful comments during the completion of this work. This study is supported in part by a grant (to NVG) from the National Institutes of Health (GM127390).

CO N FLI C T O F I NTE R E S T
None declared.

A PPE N D I X A H Y B R I D I NTER AC TI O N M O D EL
To model hybrid interactions, we select a small number of loci in autosomes for which p i,1 ∼ 1 and p i,2 ∼ 0 to interact negatively with loci in the Z chromosome(s) for which p j,1 ∼ 0 and p j,2 ∼ 1; we then repeat this process with the population subscripts interchanged, selecting an equal number of loci in autosomes with p i,2 ∼ 1 and p i,1 ∼ 0 to interact negatively with loci in the Z chromosome(s) for which p j,2 ∼ 0 and p j,1 ∼ 1. Figure A1 illustrates a pair of such interactions in a male F1 hybrid genome formed at the time of contact between diverged populations in our model. Chromosomes from populations 1 and 2 are colored red and green (respectively) with gene segments denoted by shaded blocks, and selected loci in genes denoted by vertical lines. Typically, p l = 1, 0 for selected loci, however, this is not always true for the simulations in Figure 6a due to the shorter time spent in allopatry.

A PPE N D I X B LOWE R M I G R ATI O N R ATE S
As the migration rate N is decreased, a point is reached where plots of I A versus F A and I Z versus F Z explicitly resemble those in Figure 1, concurrent with large values of ΔF ( Figures B1 and B2).

F I G U R E A 1
Hybrid genome formed at the point of contact between model populations. Chromosomes from populations 1 and 2 are colored red and green, with gene segments indicated by shaded blocks. The illustration shows a pair of interactions (dotted lines) connecting selected alleles (solid lines) in autosomes (right) to selected alleles in the Z chromosomes (left) in a male F1 hybrid F I G U R E B 1 Study of F Z versus F A for low migration rates. Data points in panel (a) denote averages of F Z for migration rates N = 0.125 (circles), 0.25 (squares), 0.5 (triangles), and 0.75 (crosses). Each set of averages is computed from 128 replicate simulations with N = 10 4 , 2 = 10 − 4 , r, r � = 10 − 2 and Δt C = N. Panel (b) describes the fraction of simulation paths that reach a given bin for F Z

A PPE N D I X C WE A K LI N K AG E B E T WE E N G E N E S
As the rate of crossing over between genes is increased from its value for closely linked genes, r � = r, differences between F Z versus F A decrease ( Figure C1).

F I G U R E B 2
Averages for I A versus F A and I Z versus F Z obtained from the simulations in Figure B1 (a) (b) F I G U R E C 1 Study of F Z versus F A as the rate of crossing over between gene sequences is increased. Data points denote averages of F Z versus F A for crossover rates r � = 0.01 (circles), and r � = 0.1 (squares). Averages are computed from 128 replicate simulations with s = 0.01, N = 10 4 , 2 = 10 − 4 , r = 10 − 2 , and Δt C = N