Incomplete recombination suppression fuels extensive haplotype diversity in a butterfly color pattern supergene

Supergenes can evolve when recombination-suppressing mechanisms like inversions promote co-inheritance of alleles at two or more polymorphic loci that affect a complex trait. Theory shows that such genetic architectures can be favoured under balancing se- lection or local adaptation in the face of gene flow, but they can also bring costs associ- ated with reduced opportunities for recombination. These costs may in turn be offset by rare ‘gene flux’ between inverted and ancestral haplotypes, with a range of possible out- comes. We aimed to shed light on these processes by investigating the BC supergene, which underlies three distinct wing colour morphs in Danaus chrysippus, a butterfly known as the African monarch, African queen and plain tiger. Using whole-genome resequencing data from 174 individuals, we first confirm the effects of BC on wing col- our pattern: background coloration is associated with SNPs in the promoter region of yellow, within an inversionted part of the supergene, while forewing tip pattern is most likely associated with a copy-number-variable part of the same supergene. We then show that haplotype diversity within the supergene is surprisingly extensive: there are at least six divergent haplotype groups that experience suppressed recombination with re- spect to each other. Despite high divergence between these haplotype groups, we identify an unexpectedly large number of natural recombinant haplotypes. These evid- ently arose through crossovers between adjacent inversion ‘modules’ as well as through double crossovers within inversions. Furthermore, we show that at least one of the es- tablished haplotype groups probably arose through recombination between two pre-ex- isting ones. Moreover, on at least two occasions, double crossovers within an inversion have led to the transfer of alleles for dark colouration in the promoter of yellow onto a different haplotype background. Overall, our findings paint a picture of dynamic evolu- tion of supergene haplotypes, fuelled by incomplete recombination suppression. 2 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 2 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint

Genomic architectures that promote co-inheritance of certain combinations of alleles at multiple loci can be beneficial if they maintain adaptive allele combinations in the face of genetic mixing [1]. Such genetic architectures are referred to as supergenes, because they facilitate the inheritance of complex phenotypes in a simple Mendelian fashion, (re- viewed in [2]). Supergenes are often associated with chromosomal inversions, which have the unique feature of suppressing recombination between haplotypes with distinct inversion orientations while allowing free recombination between haplotypes with the same orientation. This effectively divides the population into two subpopulations over a defined portion of the genome [3]. Inversion supergenes have been found to underpin trait polymorphisms under balancing selection, such as alternative life-history and repro- ductive strategies [4–8] and mimetic coloration [9]. Inversions can also be favoured dur- ing local adaptation in the face of gene flow between two distinct environments if they maintain locally adapted combinations that differentiate ecotypes [10–12]. Here we in- clude these cases of locally adapted inversions under the umbrella term ‘supergenes’, provided they still act to maintain allele combinations that would otherwise be broken down through gene flow and recombination. Genomic studies of species that exhibit local adaptation frequently uncover inversions underpinning differences between eco- types [13–16], and ‘bottom-up’ analyses have identified signatures consistent with loc- ally adapted inversions even when the specific traits they may be associated with are unknown [17–19]. This suggests that supergene architectures may be ubiquitous in nature, highlighting a need to better understand their evolutionary dynamics. Recent theoretical and empirical work suggests that supergenes may evolve over time in their structure, composition, and effects on phenotype and/or fitness. First, supergene structure can evolve through recurrent chromosomal rearrangements. Studies across multiple systems show that supergenes often comprise several rearrangements, imply- ing stepwise growth in the region of suppressed recombination (e.g. in Heliconius but- terflies [20], Danaus butterflies [21], and Solenopsis fire ants [6]). In addition to allowing the incorporation of additional co-adapted alleles at other loci, subsequent rearrange- ments of the supergene region can cause recombination suppression between more than two distinct haplotypes [21]. Second, even without physical expansion, the reper- toire of traits affected by an inversion supergene may expand over time as additional al- leles become established at other loci within the region of recombination suppression [11,17]. The fitness consequences of a supergene could also change if selective pres- sures change over time or across space, potentially limiting the value of a supergene in a changing environment or during dispersal into a new environment [22]. Even in a stable environment, inversion supergenes may be subject to accumulation of increased mutational load compared to the rest of the genome due to their reduced opportunities for recombination and reduced effective population size (Ne) (due to the effective subdi- vision of the population in that part of the genome) [23–25]. Some theoretical models 3 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 3 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint suggest that this may lead to heterokaryotype advantage through sheltering of recess- ive deleterious mutations (associative overdominance) [23–25], or failure of locally ad- apted inversions to reach high frequency [26]. Another process that contributes to the evolution of inversion supergenes is rare recom- bination in heterokaryotypes, which results in some degree of gene flux between haplo- types. While single crossovers within inversions are usually strongly suppressed due to the production of unbalanced chromosomes [27], large inversions may allow for double- crossovers, which result in the exchange of genetic material between arrangements, ef- fectively creating a mosaic haplotype. In addition, gene flux of short fragments can oc- cur through non-crossover gene conversion, which has been shown experimentally in Drosophila to occur at least as often within inversion heterokaryotypes as in syntenic re- gions [28–30]. Population genomic analyses in a range of different species support the existence of considerable gene flux within inversions [31–34], and suggest that a drift- flux equilibrium can be reached [3]. Gene flux provides a means by which supergenes may avoid some of the costs of recombination suppression described above, potentially facilitating their long-term persistence [6,23,34]. In some supergenes, sequence diver- gence suggests persistence of polymorphism for millions of years, even through mul- tiple speciation events [21,34]. However, the role of gene flux in supergene temporal dy- namics remains under-explored. A compelling example of recombination creating a third distinct haplotype is seen in the supergene that controls reproductive morphology and behaviour in the ruff [5,34]. On the other hand, a mimicry supergene in Papilio butter- flies shows phylogenetic relationships consistent with wholesale ‘allelic turnover’, in which new haplotypes arise (possibly via recombination) and replace ancestral ones [35]. Case studies such as these provide valuable insights into the range of processes contributing to supergene evolution. Danaus chrysippus, a butterfly known as the African monarch, African queen, and plain tiger, presents an opportunity to investigate the contributions of the above processes to the evolution of a large, ancient supergene. Like other milkweed butterflies, D. chrysip- pus has bright warning patterns that advertise its toxicity. Across its range, it is divided into several parapatric morphs with distinct warning patterns that meet in a broad hybrid zone in eastern central Africa (Fig. 1A) [36,37]. Phenotypic differences in two forewing traits are controlled by the ‘BC supergene’ on chromosome 15, which links at least two colour patterning loci, the ‘B’ locus which affects whether background colour is dark (dominant B allele) or pale (recessive b allele); and the ‘C’ locus which determines whether the apical black tip and white band is present (recessive c allele), or absent (dominant C allele) [21,38]. Previous genomic comparisons of a limited number of popu- lations revealed the existence of three divergent haplotype groups (also called ‘alleles’) at the BC supergene, which were named according to the morphs in which they are found: ‘orientis’ (Bc, Southern Africa), ‘klugii’ (bC, East Africa) and ‘chrysippus’ (bc, 4 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 4 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint West Africa, North Africa, Mediterranean, Asia) [38]. Note that while a distinct morph ‘al- cippus’ is found in West Africa, it only differs in its hindwing phenotype, controlled by a locus on a different chromosome, but shares its forewing phenotype and corresponding BC supergene haplotype with the chrysippus morph. Chromosome-scale assemblies of each of the three haplotypes revealed that, rather than comprising a single large inver- sion, the BC supergene is underpinned by a complex, modular rearrangement including several inversions and a large copy-number variable (CNV) region [21]. This complex architecture helps to explain how recombination is suppressed between more than two divergent haplotype groups. Phylogenetic analysis showed that the rearrangements probably occurred in a stepwise manner, beginning several million years ago, and that the supergene has persisted through multiple speciation events [21]. Intriguingly, an ad- ditional rearrangement has recently occurred and risen to high frequency in a small re- gion of East Africa: a fusion of chr15 to the female-specific W chromosome [38,39]. Be- cause crossing-over is limited to males in Lepidoptera, this fusion represents an addi- tional mechanism of recombination suppression, effectively creating a new female-spe- cific sub-group of the chrysippus haplotype [38]. There is also evidence for rare recom- bination within the supergene based on phenotypes in crossing experiments, and a pu- tative recombinant haplotype assembly from a hybrid-zone individual [21,40]. Taken to- gether, the above findings suggest that, despite its age, the BC supergene continues to evolve, raising questions about the full extent of diversity at this locus and the role of re- combination in shaping this diversity. Here we investigate haplotype diversity and evolution of the BC supergene using gen- omic data from 174 D. chrysippus individuals representing 14 regions across Africa and Southern Europe. We first confirmed that the supergene links loci controlling two distinct wing colour pattern traits. We then describe haplotype diversity and identify additional divergent haplotype groups beyond those previously described. Comparable levels of genetic diversity in each haplotype group and across most of the discrete structural modules of the supergene argues against ongoing allelic turnover, and instead suggest a long-term polymorphism, probably driven by local adaptation. Perhaps most in- triguingly, we find abundant evidence of recombination and gene flux between haplo- type groups, with crossovers having occurred both between adjacent inversions and within individual inversions, leading to exchange of functional colour pattern alleles. These results support a role for somewhat rare recombination in the diversification and long-term persistence of supergene haplotypes. 5 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 5 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint

Widespread sampling confirms that the BC supergene is associated with forewing col- our variation and is the main axis of genetic variation We analysed resequenced genomes of 174 individual butterflies spanning the three core D. chrysippus forewing morphs, ‘chrysippus’ (note that for our purposes, this also includes form ‘alcippus’, which differs only in its hindwing), ‘orientis’, and ‘klugii’, as well as intermediate morphs (Fig. 1A; S2 Table). Our sampling includes areas of mono- morphism for each morph: western Africa, South Africa and southeastern Kenya, re- spectively. We also sampled the hybrid zone (sampled in Rwanda and central Kenya, where all three forewing morphs and intermediates are found), and North Africa and the Mediterranean (where two of the three morphs, and intermediates are found). Twelve sequenced individuals were bred in captivity from parents of mixed or unknown origins and were not assigned to any geographic group. Pairwise FST between three mono- morphic locations in South Africa (TSW population, orientis morph), Nigeria (NGA popu- lation, chrysippus morph) and Eastern Kenya (WAT population, klugii morph), is close to zero genome-wide, with the exception of a few large peaks, including the largest on chromosome 15 (chr15), as described previously (Fig. 1B; S2 Table)[38]. Genome-wide principal components analysis (PCA) for all autosomes excluding chr15 further supports minimal genetic structure outside of the BC supergene, with the exception that samples from the small island of St. Helena and north of the Saraha cluster separately from the rest of sub-saharan Africa (Figure A in S1 Text). 6 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 6 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint Figure 1. Broad geographic distribution of Danaus chrysippus morphs and our continent-wide sampling of individuals spans wing pattern variation that is under- pinned largely by variation on chromosome 15. A. The geographic distribution of the three main D. chrysippus forewing morphs; chrysippus, orientis, and klugii as well as putative heterozygotes (data from manual curation of GBIF record and other scientific collections, following the approach in [37]; left) and sampled locations of D. chrysippus individuals used for genomic analysis in the present study (right). Morphs and their corresponding inferred genotypes at the B and C loci are shown below. Note that for our purposes the forewing morph ‘chrysippus’ includes the west-African morph ‘alcippus’, which has the same forewing phenotype but differs in its hindwing phenotype. B. Genome-wide patterns of FST 7 176 177 178 179 180 181 182 183 184 7 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint between morphs (smoothed across discrete 50kb windows) highlighting significant differentiation along the length of the chr15 supergene region. C-D. Genome-wide associations between SNP variation and C.

colouration and D. forewing band presence/absence (both N=172) reflecting significant asso- ciations between this region and wing-pattern variation. Significant associations, as identified via permuta- tion tests (p < 0.05), are represented by black points. Genome-wide association (GWAS) analysis confirms that the BC supergene region on chr15 is associated with two forewing traits: the pale/dark background wing colour (also known as the B locus; Fig. 1C), and the presence/absence of the apical black tip and white band (also known as the C locus; Fig. 1D). In addition to a large peak of associ- ation in the supergene, there are several peaks of association with pale/dark variation on other chromosomes, suggesting that other loci also contribute to this trait. A repeat of the GWAS using only samples from the hybrid zone (NYA, NRB, MPL populations, N=87) recapitulates the same general result (Figure B in S1 Text), confirming that the observed associations are not driven by correlated selection on other traits that follow a similar geographic distribution. The cluster of SNPs most strongly associated with back- ground wing colour fall between two genes: yellow, which is known to be involved in melanin synthesis and coloration in other insects [41], and the achaete-scute complex protein T3-like, which has been shown to be involved in wing scale development across Lepidoptera [42]. The SNPs most strongly associated with the presence or absence of the forewing black tip were found to fall within the copy-number-variable (CNV) region of the supergene (Figure C in S1 Text). These findings suggest that the B and C loci are distinct, and that the supergene is favoured because it maintains locally adapted com- binations of alleles at these loci, and probably others among the ~150 genes it encom- passes, in the face of gene flow. Extensive genetic variation at the BC supergene We previously described three divergent haplotype groups (sometimes called ‘alleles’) of the BC supergene, corresponding to the three forewing morphs, maintained by sev- eral inversions, intra-chromosomal translocations, and copy-number differences that suppress recombination [21,38]. We therefore expected that diploid genotypes across the BC supergene would form six distinct genetic clusters corresponding to the three homozygous and three heterozygous states. However, as we describe below, several lines of evidence suggest that there are additional divergent haplotypes, and possible recombinants. First, both a distance-based neighborNet network and principal compon- ents analysis for the complete supergene (excluding the CNV region) show more com- plex genetic structuring than we expected (Fig. 2A, Figure A in S1 Text). Three clusters representing homozygous genotypes for the three previously-described haplotype groups are identifiable (Fig. 2A), and include individuals previously matched to these genotypes [21]. One large cluster of individuals is intermediate between the known chrysippus and klugii homozygotes (Fig. 2A) and includes known heterozygotes 8 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 8 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint between chrysippus and klugii from crosses and suspected wild heterozygotes based on wing patterns. However, numerous other individuals are dispersed across the net- work at varying distances between the three homozygous clusters. We therefore hypo- thesised that in addition to heterozygotes among the three known haplotype groups, our sampling may include additional divergent haplotype groups and/or recombinants among the three common haplotypes. Figure 2. Genetic clustering and ancestry painting suggest additional divergent haplotype groups, and recombinants. A. NeighborNet network for unphased diploid genotypes across the BC supergene (excluding the CNV region). Each tip represents one diploid individual, and the network is constructed based on average pairwise genetic distances considering both haplotypes in each individual. We therefore expect ‘heterozygous’ individuals carrying two distinct haplotypes to be at inter- mediate positions in the network. The single individual represented by a black dot is the Danaus melanip- pus outgroup. B. Map of localities for sequenced individuals. Pie charts in panels A and B represent in- ferred ancestry components for each diploid individual from Admixture analysis with k=4 source popula- tions (See Figures D-H in S1 Text for plots with other values of k). Note that for highly sampled localities, an arbitrary subset of sequenced individuals is shown in panel B. C. Ancestry painting across the central 9 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 9 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint portion of chr15 including the BC supergene for 27 representative individuals, including homozygotes, heterozygotes and putative recombinants. The 27 selected individuals are indicated in panel A. See Fig- ures I and J in S1 Text for ancestry painting for all individuals. The CNV region is excluded from the plot for convenience as it cannot be reliably genotyped (represented as a white gap). Given the complex relationships observed, we attempted to naively identify distinct hap- lotype groups using Admixture [43] analysis, in which each individual is modelled in terms of its ancestry components from k source populations. Because recombination should occur freely among haplotypes with the same structural arrangement, and should only be suppressed among those with different arrangements, the supergene re- gion should comprise a set of semi-isolated sub-populations. We predicted that “hetero- zygous” individuals (i.e. carrying two divergent haplotypes) will appear as admixed, with approximately 50% ancestry contributions from two source populations, while homozy- gotes should be assigned 100% ancestry from a single source. All models with k=3 (matching our a priori hypothesis) and above captured the three previously-described haplotype groups corresponding to the klugii, chrysippus and orientis morphs. These are represented in clusters of homozygous (100% ancestry) individuals from eastern, western and southern Africa, respectively, while most hybrid-zone individuals appear heterozygous, with 50% ancestry proportions (Figures E-H in S1 Text). However, as de- scribed below, the models with k=4 to k=6, which all have lower cross-validation error than k=3 (Figure D in S1 Text), provide compelling evidence for the existence of addi- tional divergent haplotype groups that had not been previously described. We have chosen to present the ancestry assignments according to the model with k=4 sources in Fig. 2 (see pie charts in Fig. 2A and 2B), because our sampling included ho- mozygotes for all four putative haplotype groups. The fourth source population is in- ferred to contribute ~100% ancestry to two South African individuals, and ~50% ances- try to several others. We therefore infer that a fourth arrangement of the supergene oc- curs in southeast Africa. This haplotype group is cryptic in the sense that it is associated with a forewing phenotype indistinguishable from that of orientis. Hereafter, we refer to this as the ‘karamu’ haplotype group. In the model with k=5 sources (Figure G in S1 Text), the fifth cluster captures the previ- ously described neo-W fusion - a lineage of chr15 haplotypes that arose recently in East Africa when a copy of chr15 from the chrysippus haplotype group fused to the female- limited W sex chromosome [38,39]. These females are each assigned ~50% ancestry from this source, consistent with W being haploid in females. We previously showed that the neo-W-chr15 fusion occurred recently and spread to high frequency in southern Kenya over the past 2,200 years [38]. The high sequence similarity and complete lack of recombination of the neo-W explains how it is identifiable as a distinct genetic cluster despite its recent separation from the standard chrysippus haplotype group. 10 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 10 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint In the model with k=6 sources (Figure H in S1 Text), several individuals from North Africa and the Mediterranean are assigned to a sixth distinct cluster, suggesting yet an- other divergent haplotype group. We hypothesise that, similar to the karamu haplotype from southeast Africa, another distinct arrangement of chr15 likely occurs north of the Sahara. However, subsequent analyses (described below) suggest that none of our se- quenced individuals are homozygous for this haplotype, limiting our ability to further in- vestigate its status. Ancestry painting confirms additional haplotype groups, as well as recombinants To investigate fine-scale haplotype ancestry and recombination across the BC su- pergene region, we applied two ancestry painting methods. These aim to assign ances- try for phased-inferred haplotypes based on pre-defined source populations, using either a hidden markov model [44] or defined windows of SNPs (see Methods for de- tails). We used populations WAT, TSA and NGA as reference sets representing the klu- gii, orientis and chrysippus haplotype groups, due to their monomorphic clustering (Fig. 2B). In addition, we used the two individuals from the SPA population showing 100% karamu ancestry to represent this fourth haplotype group (Fig. 2B). Both ancestry paint- ing methods agree strongly with our inferences from the Admixture analysis. All four haplotype groups occur as intact haplotypes in most individuals, with large numbers of both homozygotes (individuals 1-11 in Fig. 2C) and heterozygotes (individuals 12-18 in Fig. 2C). See Figures I and J in S1 Text for all individuals. Heterozygotes are particu- larly common in the hybrid-zone (populations NYA, MPL and NRB in Rwanda and cent- ral Kenya, Figures I and J in S1 Text). Ancestry painting also identifies individuals with mosaic ancestry consistent with recom- bination between the divergent haplotypes. Most of these putative recombinant haplo- types appear to be simple chimeras that result from a single crossover located between two adjacent inversions (see for example individuals 19-21 in Fig. 2C). However, there are also instances of more fine-scale mosaic ancestry consistent with recombination within inversions (see for example individuals 22-25 in Fig. 2C). While gene conversion can lead to the transfer of small haplotype tracts within inversions (tract lengths vary across taxa but estimates suggest they are typically < 4kbp in length; [45]), this process cannot explain the scale of mosaic ancestry tracts we observe (tens to hundreds of kb). Instead, these patterns are consistent with occasional double crossovers, which allow genetic exchange within the inversions without the formation of unbalanced gametes (see Discussion). Some of these recombinant haplotype patterns appear multiple times in unrelated individuals, implying that they may be increasing in frequency in the popula- tion (see for example individuals 23-24 in Fig. 2C). Unsurprisingly, most recombinant haplotypes are found in the hybrid zone (Figures I and J in S1 Text). However, we were surprised to find two distinct recombinant haplotypes present in the genomes from the 11 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 11 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint small island of St. Helena (STH population), where only one of the six sequenced indi- viduals did not carry a recombinant haplotype (Figures I and J in S1 Text). Ancestry painting failed to assign ancestry to one of the haplotypes in some individuals from North Africa and the Mediterranean (see for example individuals 26-27 in Fig. 2C). These are individuals that were assigned to the 6th cluster in the admixture analysis with k=6 sources (Figure H in S1 Text). This supports the existence of an additional di- vergent haplotype of the BC supergene that occurs north of the Sahara. Other members of this population carry a chrysippus-like haplotype, with some evidence for recombina- tion with orientis (seen in individuals 25-27 in Fig. 2C). Because no individual is homo- zygous for the additional divergent haplotype, we were unable to include it as a source for ancestry painting, and it was not included in our subsequent analyses. In summary, our findings suggest that there are (at least) six divergent haplotype broups of the BC supergene that experience suppressed recombination with respect to each other. Three were previously assembled and found to have distinct arrangements (chrysippus, klugii and orientis); one results from fusion of chr15 to the W chromosome (neo-W-chrysippus); and two are newly identified here based on genetic clustering, but probably also represent distinct structural arrangements (karamu and the unnamed hap- lotype from North Africa and the Mediterranean) (Table 1). 12 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 12 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint Table 1. Six divergent haplotype groups of the BC supergene Haplotype Name Geographic re- gion(s)

col- our (B locus gen- otype) Forewing tip and band (C locus genotype) Conventional name of morph Origin chrysippus West Africa, North Africa & Mediterranean, Asia Pale (b/b) Present (c/c) chrysippus, alcippus Ancient klugii East Africa Pale (b/b) Absent (C/-) klugii Ancient orientis Southern Africa Dark (B/-) Present (c/c) orientis Ancient karamu Southeast Africa Dark (B/-) Present (c/c) orientis Recombination: klugii + orientis neo-W-chrysippus East Africa (Hybrid Zone) Pale (b/b) Present (c/c) N/A (never homozygous) W-fusion: chrysippus Unnamed North Africa & Mediterranean Unknown (no homozygotes) Unknown (no homozygotes) Unknown (no homozygotes) Ancient Diversity and divergence patterns are consistent with long-term polymorphism One possible explanation for the coexistence of multiple haplotypes associated with the same wing phenotype in the same population (e.g. karamu and orientis haplotypes in the SPA population) is that we are witnessing an ongoing ‘allelic turnover’ event [35], in which a new haplotype is replacing an older one. If a new haplotype arose through a novel structural rearrangement, it is expected to have undergone a bottleneck of N=1 at its conception [3], which would result in dramatically reduced diversity immediately fol- lowing its origin. Older haplotypes should have had the opportunity to recover diversity following their initial bottleneck, but are nevertheless still expected to have lower di- versity than the genomic background level because inversions effectively divide the species into sub-populations with lower effective population size across the inverted re- gion (even for the ancestral orientation) [3]. We therefore expected that all haplotype groups may show reduced within-group diversity within the inverted regions, but we hy- pothesised that a more pronounced reduction may be observed in those with derived in- version orientations: orientis and chrysippus [21] and possibly the newly identified karamu haplotype. Consistent with expectations, diversity is lower within the inverted re- gions of chr15 compared to collinear regions in all four haplotype groups (Fig. 3b). The 13 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 13 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint extent of diversity reduction is notably similar in the four groups, with the exception that orientis has lower diversity in inverted region 1.1, which is known to be uniquely inverted in orientis. No haplotype group has diversity approaching zero, as would be expected under rapid, recent allelic turnover, so these patterns are more consistent with a long- term polymorphism at the BC supergene. Figure 3. The supergene region has reduced genetic diversity in all haplotype groups and elevated genetic differentiation between them, but no evidence for a very recent origin of any haplotype. A. A graphical representation of the structural variation across the chrysippus, orientis, and klugii haplotypes described in [21]. The supergene evolved through rearrangements of four main regions (indicated with different colours, regions with an inverted orientation relative to the ancestral state are indicated with *). Regions 1.1, 1.2, 2, and 4 remain in single-copy and can be reliably aligned and genotyped. B-E. Population summary statistics computed in non-overlapping 50kb windows for representative populations of chrysippus (NGA population), orientis (TSW), klugii (WAT) and karamu (two homozygous individuals from the SPA), plotted across chr15 (excluding the CNV region). Note that the reference genome is from a klugii haplotype, hence regions 1 and 2 are immedi- 14 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 14 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint ately adjacent. Statistics plotted are: B. Nucleotide diversity (π) C. absolute pairwise divergence (dXY), D. net pairwise divergence (da), and E. pairwise differentiation (FST). Genetic divergence and differentiation between haplotype groups are also consistent with relatively ancient origins. Elevated FST and net divergence (da) are restricted to the rearranged parts of the chromosome (Fig. 3), consistent with recombination suppres- sion being restricted to the supergene (note that regions 1.1 and 4 are only inverted in orientis and are collinear in klugii and chrysippus). Similar levels of differentiation are seen between these three previously-described haplotype groups across most of the supergene (excluding region 1.1), which may indicate that the arrangements all oc- curred around the same time, although this pattern would also be expected under muta- tion-drift-flux equilibrium [3]. However, divergence and differentiation are both notably lower between klugii and the newly-described karamu haplotype group across most of the supergene, consistent with either a more recent divergence between these haplo- types, or a higher rate of gene flux between them. This is intriguing, because these hap- lotypes are associated with distinct wing phenotypes. This raises the possibility that gene flux between haplotype groups has allowed the exchange of functional alleles. Recombination between supergene haplotypes has exchanged wing colour alleles We set out to test whether the evolution of the karamu haplotype may have arisen through recombination between haplotype groups. Specifically, we applied topology weighting [46] using TWISST2 to ask whether there are tracts in which karamu is more closely related to orientis than to klugii. This revealed that karamu clusters broadly with klugii throughout the length of the supergene, but there are multiple narrow tracts at which it clusters more closely with orientis, particularly in inverted region 2 (Fig. 4A). An- cestry painting using Loter [44] confirms that one such narrow tract is a 20kb region containing the gene yellow and part of its promoter, including nine of the ten SNPs most strongly associated with wing background colouration according to our GWAS (Fig. 4B). This therefore suggests that incorporation of the orientis-like allele at yellow (i.e. the B allele at the B locus) through recombination causes the karamu haplotype to produce a dark wing colour phenotype matching that of orientis, despite otherwise being more klu- gii-like in its overall ancestry. We hypothesise that a similar allelic exchange at the C locus causes karamu individuals to share the forewing black tip phenotype of orientis, but we are unable to test this hypothesis, as this trait maps to the CNV region in which genotyping and ancestry assignment are unreliable. 15 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 15 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint Figure 4. Two independent cases of recombination within the BC supergene with phenotypic consequences. A. Topology weightings across chr15 showing how the karamu haplo- type is related to the klugii and orientis haplotypes. Upper panel shows three possible rooted genealogical topologies. Second panel shows weights for each topology along the chromosome, smoothed with a 20kb span. Arrows above the plot indicate the locations of inversions. Third panel shows unsmoothed topology weightings across a 1.5 Mb region corresponding to Inversion 2. B. Ancestry painting from Loter [44] across a 100 kb region within Inversion 2 showing ancestry tracts for two homozygous karamu individuals 16 400 401 402 403 404 405 406 16 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint compared to two representative individuals homozygous for the orientis and klugii haplotypes. Coding re- gions are indicated below the plot, with the candidate gene for background coloration yellow indicated. Green triangles represent the top 10 SNPs for background colour in our GWAS (Fig. 1C). There is evid- ence for recombination throughout the supergene region, and specifically in the vicinity of yellow, consist- ent with the hypothesis that orientis ancestry at this locus (i.e. the B allele) is associated with darker col- ouration in karamu individuals. C and D. A second example of recombination in the promoter of yellow. Plots are as described for panels A and B, except showing relationships between two individuals from North Africa & Mediterranean with a chrysippus-like haplotype (according to Admixture and ancestry painting), but dark background colouration. Again, orientis ancestry in the promoter region of yellow sug- gests that recombination allowed the transfer of the B allele into a different genetic background, causing darker wing colouration. We identified a second probable case of recombination affecting background coloura- tion. Several individuals from North Africa and the Mediterranean are homozygous for chrysippus-like haplotypes according to chromosome painting (see for example indi- vidual 25 in Fig. 2C), but have dark background colouration similar to that seen in ori- entis. We again tested for fine-scale mosaic ancestry in two of these individuals using topology weighting and ancestry painting. A very similar pattern to the karamu case is observed: these individuals cluster strongly with chrysippus throughout the supergene, but carry narrow orientis-like tracts (Fig. 4C), including a 10kb tract in the promoter re- gion of yellow encompassing all ten SNPs most strongly associated with dark back- ground colouration (Fig. 4D). This represents a distinct event in which orientis-like al- leles (i.e. the B allele at the B locus) have been incorporated into a different genetic

through recombination within an inversion, leading to an altered phenotype. Historial gene flux also occurred between the major supergene haplotypes Based on the extensive evidence for recombination observed, we hypothesised that even the three originally described haplotype groups (orientis, klugii and chrysippus) may have been historically shaped by gene flux. To explore this possibility, we first ran another topology weighting analysis to examine fine-scale relationships between these three groups. This generally agrees with the supergene structure: in each rearrange- ment, the predominant genealogy clusters groups according to whether they carry the rearranged or ancestral arrangement (Figure L in S1 Text). However, there is hetero- geneity in these relationships, and occasionally, complete switches to a different rela- tionship within each rearranged region (Figure L in S1 Text). These patterns are con- sistent with historical gene flux through double-crossover events. However, considering the high absolute divergence observed between these three haplotype groups (Fig. 3), there is no compelling evidence for very recent exchange of large haplotype blocks through double-crossovers. Finally, we attempted to quantify gene flux between the three major haplotype groups by fitting isolation-with-migration (IM) models, which are conventionally applied to model 17 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 17 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint species divergence with gene flow. Here haplotype groups are modelled as panmictic populations that became isolated at some point in the past, and gene flux between the distinct haplotype groups is modelled as ‘migration’ [33]. IM models were fitted using gIMble [47], with separate models for each rearranged region and each pair of haplo- type groups, because not all regions are inverted between all pairs. For example, re- gions 1.1 and 4 are inverted in orientis only, so we expect approximately free recombin- ation between klugii and chrysippus in these regions. Of the eight pairwise tests in which the regions are inverted, three had best-fitting models without gene flux, while the remaining five had estimated effective migration rates (me) ranging from 2.42x10-11 to 3.09x10-6 (this equates to a range of 3.02x10-5 to 4.89 effective migrants per generation (Me) assuming Me = me4Ne (recipient); Fig. 5; Table 2). Higher rates are generally seen between klugii and orientis, but there is no apparent trend with the size of the re- arranged region. Of the four tests that involved collinear regions, two had the highest estimated effective migration rates, and two had low or no inferred gene flux. This is probably because these latter two collinear tracts (specifically region 1.2 between chrysippus and orientis, and region 2 between klugii and chrysippus) are nevertheless physically translocated with respect to each other, probably impairing correct meiotic pairing. Taking these results together, we conclude that there is considerable evidence for history of gene flux in at least some of the rearranged regions of the supergene, but our estimated rates of gene flux should be interpreted with caution given the small size of the regions analysed and resulting sampling noise and reduced inference power (we note that one model failed to optimise), and the modelling restriction of unidirectional gene flux. 18 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 18 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint Figure 5. gIMble analysis highlights gene flow between specific supergene mod- ule alleles. A. A graphical representation of structural variation across the chrysippus, orientis, and klu- gii haplotypes replotted with data from [21] (regions with an inverted orientation relative to the ancestral state are indicated with *). B. Coalescent modelling between chrysippus and klugii individuals (top), chrysippus and orientis individuals (middle), and klugii and orientis individuals (bottom) highlighting which of the three model types, divergence only, divergence with gene flow from population A to B, and diver- gence with gene flow from B to A had the lowest composite likelihood for each region and comparison. Box borders indicate whether regions being compared are collinear (dotted line) or inverted (solid line) between morphs. * indicates the one case where the model failed to optimise. Table 2. Summary of gIMble results highlighting the best fitting models for each region and haplotype pair combination. Populations tested are listed in the order A-B. DIV = and divergence with no migration; IM AB = isolation with unidirectional migration from A to B; IM BA = from B to A. Divergence times are in generations. 19 469 470 471 472 473 474 475 476 477 478 479 480 481 19 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint

Structural variants that suppress recombination have been found to contribute to local adaptation and the maintenance of complex phenotypic polymorphisms. However, to understand how these variants persist and change over time, it is necessary to go bey- ond simply associating genotype and phenotype, and to investigate how sequences evolve within these genomic regions. Here we show that the BC supergene of Danaus chrysippus, which underpins wing colour pattern variation, displays surprising haplotype diversity, with a greater number of divergent haplotype groups than there are described phenotypic morphs. Counterintuitively, this diversity is partly fuelled by recombination. Our findings suggest that supergene evolution is a highly dynamic process, raising a number of questions that should be considered in future structural variant research. A paradox of our study is that we report evidence for both historical and recent recom- bination in a genomic region that was first identified because of its role in suppressing crossing over. In reality, the realised rate of recombination in the BC supergene must be low, otherwise we would not have been able to describe the highly diverged haplotype groups, nor to identify recombinants between them. Indeed, although a considerable fraction of our sequenced individuals carry haplotypes showing evidence for recent re- combination, many of the recombination breakpoints are shared among individuals, im- plying that the observed recombinant haplotypes stem from a smaller number of recom- bination events. The modular structure of the supergene means that single crossovers between adjacent inversions can generate modular ancestry mosaics. However, the more complex mosaic ancestries of some haplotypes imply double crossovers within in- versions. It is probably rare for two crossovers to occur within a single inversion of just a few megabases, due to crossover interference, especially considering that chromosome map lengths of butterflies in the same family tend to average around 50 cM [48,49], im- plying an average of one crossover per meiosis per bivalent. However, we note that once a double-crossover event has occurred, transferring and effectively un-inverting genetic material between inversion haplotypes, a more complex fine-scale ancestry mo- saic can emerge through unrestricted single crossovers in subsequent generations. This process, repeated over time, is expected to lead to more gene flux near the centre of in- versions than near the breakpoints, giving rise to a characteristic “suspension bridge” pattern of divergence between inversion haplotypes, as is seen for example in Droso- phila melanogaster inversion 3R [32]. The lack of such a pattern in any of the D. chrysippus BC inversions requires further investigation. In addition to double crossov- ers, there is probably a considerable contribution to gene flux from non-crossover gene conversion, as has also been described in Drosophila [28]. However, as noted, the res- ulting flux is insufficient to cause homogenisation of allele frequencies as seen in the re- mainder of the genome, where FST ~ 0. The similar levels of divergence observed between three of the haplotype groups across the inversions suggests that a point of 20 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint mutation-drift-flux equilibrium has been reached [3]. Our findings do not provide direct evidence of the potential costs or benefits of recom- bination within the supergene region. The simplest models of supergene evolution imply that recombinant haplotypes have reduced fitness due to broken associations between co-adapted or locally adapted alleles. One piece of indirect evidence in support of this idea is the observed excess of recombinant haplotypes on the small island population of St. Helena, which we believe to be recently bottlenecked, and therefore subject to greater drift and less efficient purifying selection than mainland populations. On the other hand, it is also likely that certain recombinant haplotypes are occasionally more fit than the prevailing haplotypes, either because they provide new allelic combinations as- sociated with a fit phenotype [5,22,34], or because they allow purging of deleterious mutations [23]. Our finding that several haplotypes with the same complex ancestry mo- saics appear in multiple individuals may imply a process of allelic turnover, in which an existing ‘allele’ (or haplotype group) is in the process of being replaced by a fitter recom- binant one. Indeed, the complex pattern of relationships among even the established haplotype groups suggests that recombination and allelic turnover may have occurred repeatedly throughout the evolutionary history of the BC supergene. As a result, it is likely that none of the established haplotype groups we describe here closely resemble the ‘original’ haplotypes captured when the rearrangements first occurred. A process of dynamic turnover has been suggested in a different supergene system in Papilio butter- flies [35], and may prove with further work to be the norm for supergenes, as it is for sex chromosomes in some taxa [50]. A puzzling feature of the BC supergene is the apparent long-term persistence of mul- tiple haplotype groups with (apparently) the same phenotypic effects. Specifically, the newly described karamu haplotype found in eastern parts of South Africa is associated with a colour pattern indistinguishable from that of the orientis haplotype found through- out South Africa. Although we initially hypothesised that this may represent a case of ongoing allelic turnover (in which the karamu haplotype may eventually replace the ori- entis haplotype), this is not supported by the similarly high level of diversity seen in both haplotype groups, implying that neither has experienced a recent or ongoing selective sweep. An alternative explanation is that both of these haplotypes may be adapted to different geographic regions by influencing traits other than colour pattern. It is certainly possible that the BC supergene region, which contains approximately 150 protein-cod- ing genes, contributes to other important ecological traits. Simulations have shown that such adaptations may continue to accumulate after inversions have become established [11]. It is notable that the karamu haplotype appears to have arisen through recombina- tion between the southern orientis and eastern klugii haplotypes, and appears to be most common in a geographic region intermediate between these two (though our sampling is too sparse to be sure). It is therefore plausible that this haplotype confers an 21 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 21 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint adaptive advantage by combining an optimal warning pattern for southern Africa with other environmental adaptations suited to east Africa. A similar scenario may explain the existence of a dark-coloured morph in North Africa andn the Mediterranean with a haplotype that is highly similar to that of the pale-coloured chrysippus morph, but carries orientis-like variants at the promoter of yellow that result in dark colouration. In this case, it may be that recombination allowed for allelic replacement at one locus while maintaining the broader locally-adapted haplotype. Further investigation into the eco- logy and fitness of different genotypes and phenotypes is required to test these hypo- theses. Another factor that may shape haplotype diversity is vulnerability of inversions to the ac- cumulation of deleterious mutations. Both theoretical and empirical studies have shown that inversions can accumulate excess deleterious load as a result of their lower effect- ive population size compared to the rest of the genome and their reduced opportunities for recombination [3,23,24]. This may drive turnover if new or recombinant haplotypes have lower mutational load. It has also been suggested that load may promote poly- morphism through associative overdominance, in which heterokaryotypes carrying dif- ferent sets of deleterious recessive alleles have increased fitness relative to homokaryo- types due to masking of the deleterious recessives [23]. While examples of heterokaryo- type advantage are known [24], it is difficult to prove that associative overdominance is the cause, especially given that the conditions under which it is likely to evolve are highly restrictive [25]. In D. chrysippus, there is no compelling evidence for heterozygote advantage: the three most common haplotype groups each occur in large regions of monomorphism, and broad clines suggest that polymorphism in the hybrid zone reflects a balance between selection (local adaptation) and extensive dispersal [37] rather than heterokaryotype advantage. Further investigation is needed to specifically rule out het- erokaryotype advantage, especially for the newly identified karamu haplotype and the other divergent haplotype found north of the Sahara, which are only known from poly- morphic locations. It is worth noting that the recent fusion of chr15 (chrysippus haplo- type) to the W chromosome in East Africa [38] represents the origin of yet another diver- gent haplotype that is only present in the heterozygous state (butterfly females are ZW and do not undergo recombination), and therefore would be sheltered from the effects of deleterious recessives. Finally, it is highly probable that the unusual physical structure of chromosome 15 is at least part of the explanation for the excessive diversity at the BC supergene. Most not- ably, the chromosome carries a large copy-number-variable region which began to grow around 7.5 MYA, and now comprises over a third of the chromosome in the klugii haplo- type [21]. We hypothesise that this region, which is also rich in transposable elements, creates fertile ground for the emergence of new divergent haplotypes through two mechanisms of recombination suppression [21]: It could directly prevent crossing over 22 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 22 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint by impeding proper meiotic chromosome pairing in individuals with different comple- ments of copies, and it could indirectly promote recombination suppression by propagat- ing new inversions through ectopic recombination. Whether the physical characteristics of this chromosome also facilitated its recent fusion to the W chromosome in East Africa remains to be determined. The tendency of this one chromosome to be especially prone to rearrangement echoes a similar pattern seen in the avian chromosome 4, which has undergone both ancient rearrangements (possibly associated with an ancient su- pergene [51]), and more recent fusions and further rearrangement in the formation of neo-sex chromosomes [52–54]. An interesting feature that is shared by the D. Chrysip- pus BC supergene, the Heliconius numata P supergene [24] and the fireant social su- pergene [6], is the presence of multiple adjacent inversions that appear to share break- points. This modular organisation means that crossovers that occur between two inver- sions are not suppressed, allowing the distinct modules of the supergene to evolve somewhat independently. Modular arrangements may simply be a side-effect of ‘reuse’ of inversion breakpoints (e.g. in repeat-rich regions), but it is plausible that supergenes with this structure may be more likely to persist long-term if it avoids some of the costs of complete recombination suppression. In conclusion, we have shown that a supergene with an apparently simple phenotypic association shows unexpected diversity at the haplotype level. Our findings bring to light the nuanced relationship between supergenes and recombination, in which incomplete recombination suppression can fuel haplotype diversification, and possibly support long- term persistence. Further work is needed to understand the selective and ecological drivers underlying the observed diversity. Our study adds to a growing body of work re- vealing the indirect effects of structural variation on the evolution of surrounding se- quences [24,33,55,56]. Further work across a diverse range of taxa, along with popula- tion genetic modelling will help to uncover the generality of these phenomena. 23 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 23 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint

Data collection, whole genome sequencing and genotyping We analysed whole genome sequence data (Illumina, 150 bp paired-end) from 174 but- terflies, including 158 wild-caught butterflies representing each of the main colour morphs across Africa and the hybrid zone (Fig. 1A, S2 Table) and 12 captively reared individuals. Samples collection was conducted with the permission of landowners and under appropriate permits where relevant: NACOSTI/P15/3290/3607, NACOSTI/ P15/2403/3602, National Commission for Science and Technology, Kenya; MINEDUC/S&T/459/2017, Ministry of Education, Rwanda; MPB5667 Mpumalanga Parks and Tourism Agency, South Africa; FAUNA0615202, Northern Cape Department of Environment and Nature Conservation, South Africa; EMDEP006/17, Environment & Natural Resources Directorate, St Helena Government; Permit 604 (2017), Council of Fuerteventura, Spain. A subset of individuals had been sequenced previously [21,38,57]. For those newly se- quenced in this study, the same protocols were followed. Briefly, genomic DNA was ex- tracted from ethanol-preserved tissue using the Qiagen DNEasy Blood and Tissue Kit (Qiagen, Redwood City, CA). Illumina library preparation and sequencing (150bp, paired-end) was performed by Novogene (Cambridge, UK), using a Novaseq 6000 in- strument. Reads were mapped to the Dchry2.2 reference genome (ENA accession GCA_916720795.1; [57]) using BWA mem v0.7.17 [58] with default parameters, and PCR duplicates were removed using PicardTools v2.11.1 (https://broadinstitute.git- hub.io/picard/). Indel realignment was then carried out using GATK v.3.8 [59]. To identify any substantial variation in read depth across the dataset, mean read depth was calculated using Mosdepth [60]. One individual was found to have substantially higher read depth (SM18W01) and as a result the bam file for this individual was randomly downsampled to 30% of its original size using SAMtools ([61]; samtools view -s 0.30) to eliminate biases in genotyping. Genotyping was carried out using BCFtools [62], with genotypes on the Z chromosome (chr1) called by defining ploidy based on known sex (i.e. females as haploid and males as diploid). Genotypes were then filtered using BCFtools to remove any positions called fewer than twice and to leave only genotypes with an individual depth of >= 7 and with genotype quality of >=30. Additional filtering re- moved SNPs with >10% missing data (corresponding to a haploid count of 313 across autosomes and 224 on the Z chromosome). Population structure, diversity and divergence Genomic variation across our dataset, spanning different regions of the genome, was assessed by producing PCAs from different sets of SNPs using PLINK [63]. Two SNP sets were used, representing (i) the autosomal regions of the genome excluding chr15 24 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 24 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint as well as the Z sex chromosome chr1 (14,548,985 SNPs), and (ii) rearranged regions 1, 2, and 4 of the BC supergene on chr15, described in [21] (173,998 SNPs). Genome- wide nucleotide diversity (π), FST, and dXY were computed in non-overlapping 50kb win- dows across the genome using the script popgenWindows.py (github.com/si- monhmartin/genomics_general release 0.2) and net divergence (da) was calculated from this output (calculated as d a=d XY−(π P 1+ π P 2)/2where P1 and P2 represent Popula- tion 1 and Population 2, respectively). Genome-wide association The association between SNP variation and phenotypic variation was analysed using a genome-wide association analysis implemented in PLINK across all chromosomes. Phenotypes of scanned wings from 172 individuals were scored following [37] based on whether the background colour was pale (1), intermediate (2), or dark (3), and whether the forewing band was absent (1), partial (2), or present (3) (wings were not available to phenotype for two individuals resulting in their omission from the GWAS). Empirical p- values were computed using the default adaptive permutation method implemented in PLINK. To confirm that geographic sampling structure alone does not underpin the pat- terns of association observed across the genome we re-ran the GWAS using only the 87 individuals from the hybrid zone for which we had phenotypes (populations NYA, NRB, and MPL). Visualising variation in the BC supergene using NeighborNet To visualise genetic relationships and clustering across the BC supergene, we gener- ated a phylogenetic network based on pairwise genetic distances among diploid indi- viduals using the Neighbor-Net algorithm [64], implemented in SplitsTree [65]. Pairwise distances were computed using the script distMat.py (github.com/simonhmartin/genom- ics_general release 0.2), which computes the average pairwise distance between the four pairs of sequence haplotypes for each pair of diploid individuals. This analysis was computed using modules 1, 2 and 4 of the supergene region. Admixture We aimed to naively determine the number of distinct haplotype groups of the su- pergene by modelling each individual in terms of its ancestry from a fixed number of source populations. To this end, we ran Admixture v1.3.0 [43] on a concatenated data- set of loci from across the supergene modules 1, 2, and 4 (i.e. excluding the CNV re- gion). Admixture was run specifying a range of clusters from k=3 to k=12 (Figure D in S1 Text) and 20-fold cross-validation error (--cv=20) was computed. Phasing Phasing was carried out to assist with ancestry painting using two consecutive tools, 25 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 25 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint Whatshap v1.1 [66], a read-based phasing software, and Shapeit4 v4.2.2 [67] a statist- ical phasing software using default parameters. Ancestry painting and additional phase correction To visualise genetic clustering of phased haplotypes, we used two approaches for an- cestry painting along the chromosome. Both approaches require defined reference indi- viduals, which were selected based on the Admixture analysis above. Reference indi- viduals used for each analysis are described in the Results. The first approach was Loter [44], which uses a hidden Markov model to assign ances- try and identify breakpoints along the genome. First, phased genotypes were filtered us- ing filterGenotypes.py (github.com/simonhmartin/genomics_general release 0.2) to re- move any sites at which more than 60% of individuals are heterozygous, as such sites are likely to be caused by mis-mapping and could be prone to break tracts of ancestry. Haploid genotypes were output as counts of the minor allele (0 or 1), and Loter was run using the python API in a wrapper script (loter_wrapper.py), with options rate_vote=0, nb_bagging=20. Adjacent SNPs with identical ancestry scores were concatenated into tracts to reduce the output file size. A second approach for ancestry painting that we call ‘distPaint’ was applied to confirm the Loter results, while providing a coarser visualisation at the whole-chromosome scale. This approach uses a sliding-window along the genome and assigns ancestry for each haplotype based on genetic distance (dXY) from the reference sets of individuals. A window size of 200 SNPs was used. Ancestry was assigned if this difference was signi- ficantly lower for one reference population according to a Wilcoxon Rank Sum test (p<=0.01), otherwise ancestry was set to undefined. This was implemented using a cus- tom python script distPaint.py (github.com/simonhmartin/genomics_general). Initial visualisation of the ancestry painting from both approaches revealed obvious cases of phasing errors (Figures K in S1 Text). Such errors could lead to false-positive inference of recombination events. For example, say an individual is phased into two haplotypes, that are then painted with ancestries “c-c-c-c-k” and “k-k-k-k-c” (where c and k represent tracts of chromosome assigned to reference populations chrysippus and klugii, respectively). This could either represent a case of an individual carrying two recombinant haplotypes in which both recombination events happened to occur at the same genomic region, or it could be a case of a single phasing error in an individual that in fact caries two common haplotypes (c-c-c-c-c and k-k-k-k-k). The latter is of course far more likely. We therefore applied a final heuristic phase correction approach to the inferred ancestries which attempts to maximise similarity among haplotypes at the whole-chromosome scale, thereby minimising inferred recombination events (Figure K in S1 Text). For each diploid individual, this heuristic approach iterates over each an- 26 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 26 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint cestry-painted tract in the chromosome (e.g. from Loter) and switches the phase if this would cause an increase in similarity of the two resulting haplotypes to one or more other haplotypes in the full dataset (considering all ancestry blocks up to and including the focal one). Because we are interested in identifying common shared haplotypes, we consider the average similarity of the top three most similar haplotypes in the rest of the dataset. The process is performed left-to-right and then right-to-left for each individual, and then repeated again starting from the first individual, up to a total of 20 iterations for Loter outputs and 100 iterations for distPaint.py outputs (the latter have fewer intervals, allowing for more rapid computation). This algorithm is implemented in the script phasepaint.py (https://github.com/simonhmartin/phasepaint). Topology weighting To complement the ancestry painting described above, we analysed changes in rela- tionships among supergene haplotypes across the chromosome (as an indicator of pos- sible recombination), using topology weighting [46], implemented in TWISST2 (https:// github.com/simonhmartin/twisst2). This approach infers genealogies and where they change along the genome, and outputs a ‘weighting’ for each possible topology for the relationships between defined groups of individuals in each genomic region. The input variants file was further filtered to remove sites with more than 75% heterozygotes (which are likely genotyping errors that could interfere with tree inference) and sites with missing data for more than 50% of individuals. An outgroup was needed for polarisation, for which we used a single Danaus melanippus individual (S2 Table), following the alignment and genotyping procedures described above. Modelling gene flux using gIMble To detect and quantify gene flux between supergene haplotypes, we fitted isolation- with-migration models for individuals homozygous for the chrysippus, klugii, orientis al- leles using the gIMble framework [47]. Under the IM model, gene-flux between su- pergene haplotypes is modelled as migration (gene flow) between panmictic popula- tions. Six klugii individuals (WAT population), seven orientis individuals (TSW popula- tion), and six chrysippus (NGA population) were included in the analysis (S2 Table). Since gIMble carries out pairwise tests this resulted in six pairwise tests for each of the 4 regions tested. A VCF containing only data from these 19 individuals was first prepro- cessed using gimble preprocess. Next, intergenic bed files were produced by using bedtools intersect and bedtools subtract to extract and then exclude genic regions using the reference annotation file. gimble parse was then run for each combination of su- pergene region (i.e. regions 1.1, 1.2, 2, and 4, Fig. 3A) and pairwise comparison, fol- lowed by gimble blocks, gimble windows, gimble info, and gimble tally. gimble optimize was then run for each combination of region and pairwise test specifying each of three models, DIV - a divergence only model, IM_AB an isolation with migration model where gene flow occurs from population A into B, and IM_BA an isolation with migration model 27 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 27 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint where gene flow occurs from population B into A, resulting in a total of 36 models. gimble query was then run to extract model summary information, including composite likelihood, for each model allowing us to determine the best fitting of the three models for each of the 12 region and pair combinations (full results are reported in https://git- hub.com/RishiDeKayne/Danaus_WGS/blob/main/FINAL_GIMBLE_RESULTS_3morph- s.xlsx). For detailed commands and model parameters see https://github.com/ RishiDeKayne/Danaus_WGS/blob/main/Danaus_WGS_commands.txt.

We thank the following people for assistance with sampling, breeding and acquisition of permissions: Rwanda - Constantin Sibomana and Jody Garbe; Kenya - Piera Ireri, Ivy Ng’Iru, Godfrey Etelej, David Smith, Anna Orteu, Jenny York, Owen McMillan, Valarie McMillan; South Africa - Jeremy Dobson; Ghana - Oskar Brattstrom; Spain (Canary Is- lands) - Yeray Monasterio, David Smith; Tunisia - Roger Vila ; Italy: Richard ffrench- Constant; St. Helena: David Pryce. We also thank Brian Charlesworth, Deborah Char- lesworth, David Smith, Richard ffrench-Constant, Chay Graham, Thomas Decroly, Alex- ander Mackintosh, Domink Laetsch and Konrad Lohse for valuable input on this work. Data availability: All raw sequencing reads are available via the European Nucleotide Archive (ENA) - Sample accession numbers are provided in S2 Table. Commands and scripts for bioin- formatics analyses can be found at https://github.com/RishiDeKayne/Danaus_WGS/. Funding: This work was supported by a The Royal Society (grants URF\R1\180682, RGF\EA\ 181071, URF\R\231034 to S.H.M), the Swiss National Science Foundation (grants P2BEP3_195567, P500PB_211005 to R.D.K) and the National Geographic Society (grant WW-138R-17 to I.J.G). 28 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 28 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint

1. Charlesworth B, Charlesworth D. Selection of new inversions in multi-locus genetic sys- tems. Genet Res . 1973;21: 167–183. 2. Thompson MJ, Jiggins CD. Supergenes and their role in evolution. Heredity . 2014;113: 1– 8. 3. Charlesworth B. The effects of inversion polymorphisms on patterns of neutral genetic di- versity. Genetics. 2023;224. doi:10.1093/genetics/iyad116 4. Mérot C, Llaurens V, Normandeau E, Bernatchez L, Wellenreuther M. Balancing selection via life-history trade-offs maintains an inversion polymorphism in a seaweed fly. Nat Com- mun. 2020;11: 670. 5. Küpper C, Stocks M, Risse JE, Dos Remedios N, Farrell LL, McRae SB, et al. A supergene determines highly divergent male reproductive morphs in the ruff. Nat Genet. 2016;48: 79– 83. 6. Yan Z, Martin SH, Gotzek D, Arsenault SV, Duchen P, Helleu Q, et al. Evolution of a su- pergene that regulates a trans-species social polymorphism. Nat Ecol Evol. 2020;4: 240– 249. 7. Li J, Cocker JM, Wright J, Webster MA, McMullan M, Dyer S, et al. Genetic architecture and evolution of the S locus supergene in Primula vulgaris. Nat Plants. 2016;2: 16188. 8. Kim K-W, Bennison C, Hemmings N, Brookes L, Hurley LL, Griffith SC, et al. A sex-linked supergene controls sperm morphology and swimming speed in a songbird. Nat Ecol Evol. 2017;1: 1168–1176. 9. Komata S, Kajitani R, Itoh T, Fujiwara H. Genomic architecture and functional unit of mim- icry supergene in female limited Batesian mimic Papilio butterflies. Philos Trans R Soc Lond B Biol Sci. 2022;377: 20210198. 10. Kirkpatrick M, Barton N. Chromosome inversions, local adaptation and speciation. Genet- ics. 2006;173: 419–434. 11. Schaal SM, Haller BC, Lotterhos KE. Inversion invasions: when the genetic basis of local adaptation is concentrated within inversions in the face of gene flow. Philos Trans R Soc Lond B Biol Sci. 2022;377: 20210200. 12. Berdan EL, Barton NH, Butlin R, Charlesworth B, Faria R, Fragata I, et al. How chromo- somal inversions reorient the evolutionary process. J Evol Biol. 2023;36: 1761–1782. 13. Lowry DB, Willis JH. A widespread chromosomal inversion polymorphism contributes to a major life-history transition, local adaptation, and reproductive isolation. PLoS Biol. 2010;8. doi:10.1371/journal.pbio.1000500 14. Hager ER, Harringmeyer OS, Wooldridge TB, Theingi S, Gable JT, McFadden S, et al. A chromosomal inversion contributes to divergence in multiple traits between deer mouse ecotypes. Science. 2022;377: 399–405. 29 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 29 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint 15. Huang K, Andrew RL, Owens GL, Ostevik KL, Rieseberg LH. Multiple chromosomal inver- sions contribute to adaptive divergence of a dune sunflower ecotype. Mol Ecol. 2020;29: 2535–2549. 16. Kapun M, Fabian DK, Goudet J, Flatt T. Genomic Evidence for Adaptive Inversion Clines in Drosophila melanogaster. Mol Biol Evol. 2016;33: 1317–1336. 17. Todesco M, Owens GL, Bercovich N, Légaré J-S, Soudi S, Burge DO, et al. Massive haplo- types underlie ecotypic differentiation in sunflowers. Nature. 2020;584: 602–607. 18. Harringmeyer OS, Hoekstra HE. Chromosomal inversion polymorphisms shape the gen- omic landscape of deer mice. Nat Ecol Evol. 2022;6: 1965–1979. 19. Morales HE, Faria R, Johannesson K, Larsson T, Panova M, Westram AM, et al. Genomic architecture of parallel ecological divergence: Beyond a single environmental contrast. Sci Adv. 2019;5: eaav9963. 20. Joron M, Frezal L, Jones RT, Chamberlain NL, Lee SF, Haag CR, et al. Chromosomal re- arrangements maintain a polymorphic supergene controlling butterfly mimicry. Nature. 2011;477: 203–206. 21. Kim K-W, De-Kayne R, Gordon IJ, Omufwoko KS, Martins DJ, Ffrench-Constant R, et al. Stepwise evolution of a butterfly supergene via duplication and inversion. Philos Trans R Soc Lond B Biol Sci. 2022;377: 20210207. 22. Roesti M, Gilbert KJ, Samuk K. Chromosomal inversions can limit adaptation to new envir- onments. Mol Ecol. 2022;31: 4435–4439. 23. Berdan EL, Blanckaert A, Butlin RK, Bank C. Deleterious mutation accumulation and the long-term fate of chromosomal inversions. PLoS Genet. 2021;17: e1009411. 24. Jay P, Chouteau M, Whibley A, Bastide H, Parrinello H, Llaurens V, et al. Mutation load at a mimicry supergene sheds new light on the evolution of inversion polymorphisms. Nat Genet. 2021;53: 288–293. 25. Charlesworth B. The fitness consequences of genetic divergence between polymorphic gene arrangements. Genetics. 2024;226. doi:10.1093/genetics/iyad218 26. Jay P, Aubier TG, Joron M. The interplay of local adaptation and gene flow may lead to the formation of supergenes. Mol Ecol. 2024; e17297. 27. Sturtevant AH, Beadle GW. The Relations of Inversions in the X Chromosome of Droso- phila Melanogaster to Crossing over and Disjunction. Genetics. 1936;21: 554–604. 28. Chovnick A. Gene conversion and transfer of genetic information within the inverted region of inversion heterozygotes. Genetics. 1973;75: 123–131. 29. Crown KN, Miller DE, Sekelsky J, Hawley RS. Local Inversion Heterozygosity Alters Re- combination throughout the Genome. Curr Biol. 2018;28: 2984–2990.e3. 30. Korunes KL, Noor MAF. Pervasive gene conversion in chromosomal inversion heterozy- gotes. Mol Ecol. 2019;28: 1302–1315. 30 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 30 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint 31. Corbett-Detig RB, Hartl DL. Population genomics of inversion polymorphisms in Drosophila melanogaster. PLoS Genet. 2012;8: e1003056. 32. Kapun M, Mitchell ED, Kawecki TJ, Schmidt P, Flatt T. An Ancestral Balanced Inversion Polymorphism Confers Global Adaptation. Mol Biol Evol. 2023;40. doi:10.1093/molbev/ msad118 33. Lundberg M, Mackintosh A, Petri A, Bensch S. Inversions maintain differences between mi- gratory phenotypes of a songbird. Nat Commun. 2023;14: 452. 34. Hill J, Enbody ED, Bi H, Lamichhaney S, Lei W, Chen J, et al. Low Mutation Load in a Su- pergene Underpinning Alternative Male Mating Strategies in Ruff (Calidris pugnax). Mol Biol Evol. 2023;40. doi:10.1093/molbev/msad224 35. Palmer DH, Kronforst MR. A shared genetic basis of mimicry across swallowtail butterflies points to ancestral co-option of doublesex. Nat Commun. 2020;11: 6. 36. Smith DAS, Owen DF, Gordon IJ, Lowis NK. The butterfly Danaus chrysippus (L.) in East Africa: polymorphism and morph-ratio clines within a complex, extensive and dynamic hy- brid zone. Zool J Linn Soc. 1997;120: 51–78. 37. Liu W, Smith DAS, Raina G, Stanforth R, Ng’Iru I, Ireri P, et al. Global biogeography of warning coloration in the butterfly Danaus chrysippus. Biol Lett. 2022;18: 20210639. 38. Martin SH, Singh KS, Gordon IJ, Omufwoko KS, Collins S, Warren IA, et al. Whole-chromo- some hitchhiking driven by a male-killing endosymbiont. PLoS Biol. 2020;18: e3000610. 39. Smith DAS, Gordon IJ, Traut W, Herren J, Collins S, Martins DJ, et al. A neo-W chromo- some in a tropical butterfly links colour pattern, male-killing, and speciation. Proc Biol Sci. 2016;283. doi:10.1098/rspb.2016.0821 40. Smith DA. Evidence for autosomal meiotic drive in the butterfly Danaus chrysippus L. Heredity . 1976;36: 139–142. 41. Wittkopp PJ, Vaccaro K, Carroll SB. Evolution of yellow gene regulation and pigmentation in Drosophila. Curr Biol. 2002;12: 1547–1556. 42. Galant R, Skeath JB, Paddock S, Lewis DL, Carroll SB. Expression pattern of a butterfly achaete-scute homolog reveals the homology of butterfly wing scales and insect sensory bristles. Curr Biol. 1998;8: 807–813. 43. Alexander DH, Lange K. Enhancements to the ADMIXTURE algorithm for individual ances- try estimation. BMC Bioinformatics. 2011;12: 246. 44. Dias-Alves T, Mairal J, Blum MGB. Loter: A Software Package to Infer Local Ancestry for a Wide Range of Species. Mol Biol Evol. 2018;35: 2318–2326. 45. Chen J-M, Cooper DN, Chuzhanova N, Férec C, Patrinos GP. Gene conversion: mechan- isms, evolution and human disease. Nat Rev Genet. 2007;8: 762–775. 46. Martin SH, Van Belleghem SM. Exploring Evolutionary Relationships Across the Genome Using Topology Weighting. Genetics. 2017;206: 429–438. 31 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 31 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint 47. Laetsch DR, Bisschop G, Martin SH, Aeschbacher S, Setter D, Lohse K. Demographically explicit scans for barriers to gene flow using gIMble. PLoS Genet. 2023;19: e1010999. 48. Davey JW, Chouteau M, Barker SL, Maroja L, Baxter SW, Simpson F, et al. Major Improve- ments to the Heliconius melpomene Genome Assembly Used to Confirm 10 Chromosome Fusion Events in 6 Million Years of Butterfly Evolution. G3 . 2016;6: 695–708. 49. Shipilina D, Näsvall K, Höök L, Vila R, Talavera G, Backström N. Linkage mapping and genome annotation give novel insights into gene family expansions and regional recombin- ation rate variation in the painted lady (Vanessa cardui) butterfly. Genomics. 2022;114: 110481. 50. Vicoso B. Molecular and evolutionary dynamics of animal sex-chromosome turnover. Nat Ecol Evol. 2019;3: 1632–1641. 51. Mirarab S, Rivas-González I, Feng S, Stiller J, Fang Q, Mai U, et al. A region of suppressed recombination misleads neoavian phylogenomics. Proc Natl Acad Sci U S A. 2024;121: e2319506121. 52. Sigeman H, Strandh M, Proux-Wéra E, Kutschera VE, Ponnikas S, Zhang H, et al. Avian Neo-Sex Chromosomes Reveal Dynamics of Recombination Suppression and W Degener- ation. Mol Biol Evol. 2021;38: 5275–5291. 53. Sigeman H, Zhang H, Ali Abed S, Hansson B. A novel neo sex chromosome in Sylvietta ‐ brachyura (Macrosphenidae) adds to the extraordinary avian sex chromosome diversity among Sylvioidea songbirds. J Evol Biol. 2022;35: 1797–1805. 54. Dierickx EG, Sin SYW, van Veelen HPJ, Brooke M de L, Liu Y, Edwards SV, et al. Genetic diversity, demographic history and neo-sex chromosomes in the Critically Endangered Raso lark. Proc Biol Sci. 2020;287: 20192613. 55. Huang K, Ostevik KL, Elphinstone C, Todesco M, Bercovich N, Owens GL, et al. Mutation Load in Sunflower Inversions Is Negatively Correlated with Inversion Heterozygosity. Mol Biol Evol. 2022;39. doi:10.1093/molbev/msac101 56. Jeong H, Baran NM, Sun D, Chatterjee P, Layman TS, Balakrishnan CN, et al. Dynamic molecular evolution of a supergene with suppressed recombination in white-throated spar- rows. Elife. 2022;11. doi:10.7554/eLife.79387 57. Singh KS, De-Kayne R, Omufwoko KS, Martins DJ, Bass C, ffrench-Constant R, et al. Gen- ome assembly of Danaus chrysippus and comparison with the Monarch Danaus plexippus. G3 . 2021. doi:10.1093/g3journal/jkab449 58. Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics. 2010;26: 589–595. 59. Poplin R, Ruano-Rubio V, DePristo MA, Fennell TJ, Carneiro MO, Auwera GAV der, et al. Scaling accurate genetic variant discovery to tens of thousands of samples. bioRxiv. 2017; 201178. 60. Pedersen BS, Quinlan AR. Mosdepth: quick coverage calculation for genomes and ex- omes. Bioinformatics. 2018;34: 867–868. 32 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 32 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint 61. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Align- ment/Map format and SAMtools. Bioinformatics. 2009;25: 2078–2079. 62. Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. Gigascience. 2021;10. doi:10.1093/gigascience/giab008 63. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81: 559–575. 64. Bryant D, Moulton V. Neighbor-net: an agglomerative method for the construction of phylo- genetic networks. Mol Biol Evol. 2004;21: 255–265. 65. Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006;23: 254–267. 66. Martin M, Patterson M, Garg S, Fischer SO, Pisanti N, Klau GW, et al. WhatsHap: fast and accurate read-based phasing. bioRxiv. 2016. p. 085050. doi:10.1101/085050 67. Delaneau O, Zagury J-F, Robinson MR, Marchini JL, Dermitzakis ET. Accurate, scalable and integrative haplotype estimation. Nat Commun. 2019;10: 5436. 33 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 33 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 26, 2024. ; https://doi.org/10.1101/2024.07.26.605145doi: bioRxiv preprint