A large deletion at the cortex locus eliminates butterfly wing patterning

Abstract As the genetic basis of natural and domesticated variation has been described in recent years, a number of hotspot genes have been repeatedly identified as the targets of selection, Heliconius butterflies display a spectacular diversity of pattern variants in the wild and the genetic basis of these patterns has been well-described. Here, we sought to identify the mechanism behind an unusual pattern variant that is instead found in captivity, the ivory mutant, in which all scales on both the wings and body become white or yellow. Using a combination of autozygosity mapping and coverage analysis from 37 captive individuals, we identify a 78-kb deletion at the cortex wing patterning locus, a gene which has been associated with wing pattern evolution in H. melpomene and 10 divergent lepidopteran species. This deletion is undetected among 458 wild Heliconius genomes samples, and its dosage explains both homozygous and heterozygous ivory phenotypes found in captivity. The deletion spans a large 5′ region of the cortex gene that includes a facultative 5′UTR exon detected in larval wing disk transcriptomes. CRISPR mutagenesis of this exon replicates the wing phenotypes from coding knock-outs of cortex, consistent with a functional role of ivory-deleted elements in establishing scale color fate. Population demographics reveal that the stock giving rise to the ivory mutant has a mixed origin from across the wild range of H. melpomene, and supports a scenario where the ivory mutation occurred after the introduction of cortex haplotypes from Ecuador. Homozygotes for the ivory deletion are inviable while heterozygotes are the targets of artificial selection, joining 40 other examples of allelic variants that provide heterozygous advantage in animal populations under artificial selection by fanciers and breeders. Finally, our results highlight the promise of autozygosity and association mapping for identifying the genetic basis of aberrant mutations in captive insect populations.


Introduction
Captive populations can harbor variation that is not observed in the wild. Aside from domestication for agricultural or commercial yield, animals and plants have often been selectively bred for "beautiful," "interesting," or "unnatural" esthetic traits, especially color and pattern variants (Driscoll et al. 2009;Cieslak et al. 2011). In recent decades, the genetic basis of artificially selected variation has begun to be mapped in many organisms used in agriculture, floriculture and the pet trade, including the genetic mapping of a large number of genes responsible for flower color variation (Park et al. 2007;Giovannini et al. 2021), melanin-based coat color in domestic mammals (Cieslak et al. 2011), plumage patterns, and colors in birds (Domyan and Shapiro 2017;Price-Waldman and Stoddard 2021), and chromatophore distribution in squamates and fishes (Guo et al. 2021). These variants can affect genes involved in the enzymatic production, transport, and deposition of pigments such as the melanin and carotenoid pathways (Mundy et al. 2016;, or can be caused by genes that affect signaling or cell type differentiation [as with chromatophore distribution in squamates (Kuriyama et al. 2020)]. The extensive set of artificial and domesticated genetic variants has been used repeatedly both for studies of genotype-to-phenotype relationship, and for understanding the genetics of disease states in humans (Andersson 2016).
Notably, several cases have been identified where a gene underlying an artificially selected variant is also responsible for natural variation in other populations or species. This is the case with BCO2, where an artificially selected protein coding variant in cattle affects milk color (Berry et al. 2009), and a wild-type regulatory variant in wall lizards causes changes to ventral scale color (Andrade et al. 2019). The gene Agouti has repeatedly been mapped in cases of pigment differences, both in cases of domestication [e.g. in the chicken (Yu et al. 2019), horse (Rieder et al. 2001), and rabbit (Fontanesi et al. 2010)], and in natural variation [e.g. in warblers (Toews et al. 2016), humans (Bonilla et al. 2005), and snowshoe hare (Jones et al. 2018)]. This reuse of hotspot genes in wild and captive populations can even be observed at the within-species level; multiple separate variants in Agouti have been linked to pigment variation in sheep-both in wild and captive populations, and caused by protein coding mutations, cisregulatory mutations, and copy number variation (Norris and Whan 2008;Gratten et al. 2010). As such, studies of domesticated variation can inform our understanding of natural variation too, both on macro-and microevolutionary scales.
As the primary selective force in captive populations is the eye of the selector, one may often observe variation that is not seen in nature because it is not optimally fit (Andersson and Georges 2004;Marsden et al. 2016;Makino et al. 2018;Moyers et al. 2018). Fanciers and breeders have purposefully selected for variation that causes deleterious effects in combination with color and pattern differences, applying a regime of selection that favors the phenotype of interest over fitness. This applies especially to distinctive coloration, which can be straightforward to observe and maintain in a captive stock. Examples include coat color variants for Merle dogs, where homozygotes for a retrotransposon insertion in the SILV gene have an increased risk of deafness and blindness (Strain et al. 2009;Langevin et al. 2018), and in overo horses, where foals homozygous for a mutated Endothelin Receptor B (EDNRB) develop Lethal White Syndrome (Metallinos et al. 1998). "Lemon frost" geckos have been selected for their unique color, but exhibit an increased risk of iridophoroma, the formation of tumors from iridophores (analogous to melanoma) (Guo et al. 2021). Such reduced fitness has also been identified in floriculture, where white petunias with mutations to the AN1 gene have deficient vacuole acidification and a weakened seed coat (Spelt et al. 2002), though an artificially selected mutation to the same gene causing a similar effect in the morning glory does not appear to carry the same deleterious fitness effects (Park et al. 2007). Whether they are the direct targets of selective breeding, or collateral effects of inbreeding and bottlenecks, deleterious mutations that would be purged by natural selection in the wild are common in captive populations (Bosse et al. 2019).
The majority of work mapping artificially selected or domesticated variation in animals has taken place in vertebrates, with the notable exceptions of the economically important silkworm (Futahashi and Osanai-Futahashi 2021). Heliconius butterflies, a model system for the genetic study of color pattern adaptations in the wild Van Belleghem et al. 2021), have also been maintained in captivity by butterfly breeders continuously since at least the mid-20th century, including at notable tourist attractions like Butterfly World (BW) in Florida. There, they have been selectively bred for hundreds of generations with intentional, artificial selection upon pattern variations that deviate from phenotypes normally observed in the wild. Over an extended period of inbreeding, selection, and occasional outcrossing, novel color variations have appeared in this captive population, including the "Piano Keys" (PK) pattern, and a deleterious mutation that occurred within the PK genetic stock, here dubbed ivory (Fig. 1).
The genetics of wing patterning have been extensively studied in Heliconius (McMillan et al. 2020). Linkage and association mapping have pinpointed regulatory regions around a small toolkit of genes as being responsible for much of the pattern variation seen in wild populations, including the transcription factor optix Zhang et al. 2017;Lewis et al. 2019;Morris et al. 2020;), the signaling ligand WntA Van Belleghem et al. 2017), and cortex, a cdc20 homolog with a currently unidentified function Livraghi et al. 2021). We used whole-genome resequencing, association mapping and autozygosity mapping to determine that the ivory mutation is caused by a large deletion at the patterning gene and evolutionary hotspot cortex, which likely occurred de novo in the captive population.

Stock history
The captive population of Heliconius melpomene was initially collected by JRG Turner and others from around Central and South America, and was later transferred from his genetic research stocks at the University of Leeds to the former London Butterfly House (LBH) maintained by Clive Farrell from 1981 as a multirace hybrid population. They were then acquired from Tom Fox of LBH by R. Boender at the MetaScience Butterfly Farm in Florida in about 1985. These hybrid stocks formed a core of Boender's inhouse hybrid Heliconius display stock when BW opened in 1988. Only one known introduction of wild-caught H. melpomene occurred to this stock, of H. melpomene cythera collected in the early 1980s by R. Boender and T. Emmel around Tinalandia, Ecuador. After this introduction, Boender began selection for a novel pattern that is termed "Piano Key" in this paper. The exact mix of races that were transferred from LBH to BW are not known, but contain pattern alleles that are characteristic of H. melpomene races of Suriname and the Guianas (LEG, personal observation).
Since the Piano Key phenotype was first noticed by R. Boender over 30 years ago, the Piano Key stock population has been maintained in separate compartments for around 400 generations at BW and has long been fixed for the PK phenotype (RB, personal communication). Thus, we referred to it as H. melpomene BWPK throughout. Within the last decade, novel mutant butterflies appeared with fewer melanic scales, initially termed "pale PK." These were subsequently separated to form a selected subpopulation at BW. Soon after, virtually pure white butterflies incapable of flight began appearing in this selected stock, herein referred to as ivory. These were subjected to intentional artificial selection. In 2015, "Pale PK" were sent to UT Austin and maintained by LEG in climate-controlled greenhouses (see Supplementary Table 1 for phenotypes and collection dates) to investigate the genetics of this mutation.

Imaging
Pinned specimens were imaged with a Nikon D5300 camera mounted with a Micro-NIKKOR 105 mm f/2.8 lens. Scale phenotypes were imaged with a Keyence VHX-5000 microscope mounted with a VH-Z100R lens.
Short-read DNA sequencing DNA was extracted from thoracic tissues of 37 H. m. BWPK individuals using the Qiagen Dneasy Blood & Tissue Kit, RNAsetreated, and used to prepare a multiplexed sequencing library with the TruSeq PCR-free DNA protocol. Samples were sequenced on an Illumina NovaSeq S1 150 bp PE run, yielding 13Â mean coverage per individual. Sequencing reads are accessible via NCBI SRA under the project accession number PRJNA610063.

Genomic analyses
Samples were aligned to Hmel2.5 (Pinharanda et al. 2019) retrieved from LepBase (Challis et al. 2016) with BWA-MEM using default parameters (Li 2013), and variants called using GATK v4.1 with tools HaplotypeCaller and GenotypeGVCFs using default parameters (McKenna et al. 2010). Variant sites were accepted if they were biallelic and the quality (QUAL) value was !30. SNPs were phased with Beagle 4.1 (Browning and Browning 2007). Phased SNP variants were used to perform principal component analysis (PCA) using the Eigensoft module SmartPCA (Price et al. 2006). SNP association was carried out in PLINK v.1.9, with 1,000 permutations (Purcell et al. 2007). For phylogenetic analyses, we followed the TWISST pipeline ; briefly, data were phased in Beagle 4.1 with a window size of 10,000 and overlap of 1,000. Trees were built from windows of 50 SNPs, and then analyzed with TWISST using the 5 groups "East," "West," "Ecuador," "Atlantic," and "BWPK" (Supplementary Table 3).

De novo assembly for indel breakpoints
To find precise indel breakpoints, a subset of samples were de novo assembled with Velvet (Zerbino and Birney 2008) using default parameters. The resulting contigs were searched with BLAST (Camacho et al. 2009) for the genomic region including the deletion. Scaffolds that bridged the indel breakpoints were selected, and aligned to the Hmel2.5 reference genome with MAFFT (Katoh et al. 2002).

CRISPR/Cas9 mutagenesis
We designed sgRNAs against the putative H. erato cortex promoter corresponding to N20NGG sites within the distal promoter/5 0 UTR. In order to mutagenize the locus, we designed 3 sgRNAs (Table 1) using the "find CRISPR sites" algorithm within the Geneious software. Guide specificity and off-target effects were assessed by scoring against the H. erato reference genome. sgRNAs displaying low off-target scores were then synthesized commercially by Synthego, and mixed with Cas9 at a concentration of 500:500 ng/ ml, respectively. Embryonic injections were performed as previously described (Livraghi et al. 2021) within 1-3 h after egg laying, after which larvae were allowed to develop on a diet of P. biflora until adult emergence.

Heliconius melpomene BWPK-artificially selected wing patterns in an insectary stock
The stock maintained at BW in Florida and at UT Austin is here named H. melpomene BWPK (see Methods-Stock history section). The insectary-bred, artificially selected line is polymorphic for several wing pattern elements found in wild populations, including the red pattern elements (Fig. 1). Additionally, H. m. BWPK includes aberrant pattern features not observed in natural populations. First, in the pattern dubbed "Piano Keys" (Fig. 1a), white or yellow hindwing marginal elements extend along the veins toward the discal cell. While distal yellow hindwing pattern elements are observed in some wild pattern forms of Heliconius, including H. m. cythera from Ecuador and H. cydno, the extent of these marginal elements in H. m. BWPK is much greater than any of these, occurring over most of the area of the hindwing except for the intervein regions that would otherwise be taken by the red hindwing rays.
Occasional butterflies with more extensive regions of yellow or white scales were identified (Fig. 1a, center). Other than the wing pattern differences, both forms of butterflies exhibit typical feeding, mating, egg-laying, and flight behaviors to other insectary-reared Heliconius (e.g. Supplementary Movie 1). However, the offspring of a mating of 2 pale H. m. BWPK will include completely white-winged and -bodied butterflies (Fig. 1c). We dub these butterflies "ivory" in reference to their emergence in the "Piano Keys" stock, and inferred that the pale morph of H. m. BWPK represents a heterozygous state for a mutation that causes the ivory phenotype in the homozygous state.
The ivory butterflies have melanin pigments in their body cuticle and eyes, indicating that their capacity to synthesize and deposit melanin has not been perturbed. However, all scales on the wings and body are white or yellow, with no black or red scales, with the exception of a 1-mm patch of red at the base of the hindwing in some individuals. This is in stark contrast to wildtype or BWPK butterflies which have black scales covering most of the body. Unlike the other BWPK butterflies, ivory homozygote butterflies do not fly and have not been observed to successfully mate or lay eggs. The wings appear to be of a normal strength and sturdiness, but flight is weak and will not occur when prompted by dropping (Supplementary Movie 1).
In order to investigate the genetic basis of these wing patterns, we sequenced 37 butterflies from the UT Austin H. m. BWPK stock to an average coverage of 13Â, including 11 Dark BWPKs, 9 Pale BWPKs and 10 ivory butterflies (Fig. 1a).

Autozygosity mapping identifies associated SNPs
We used a combination of GWAS and patterns of SNP autozygosity to determine the locations of regions associated with Mendelian pattern variants (Fig. 2). As proof of principle, we mapped the Dennis pattern element which was previously described as a cis-regulatory element of the gene optix. GWAS for presence vs absence of Dennis gave a narrow association peak centered near optix, in a region previously identified as associated with the Dennis pattern element (Morris et al. 2020) ( Supplementary Fig. 1). We then examined the inheritance patterns of SNPs, filtering for SNPs where butterflies with no Dennis element were homozygous for the reference allele, and individuals with the Dennis element were either heterozygous or homozygous for an alternate allele (6/29 butterflies had the Dennis pattern element, Supplementary Table 1). In total, 3,430 SNPs matched this pattern of inheritance, with 3,425 in a cluster centered on the gene optix, indicating that autozygosity mapping is appropriate for mapping pattern elements in this data set.
GWAS for ivory gave a broad association peak on chromosome 15 (Fig. 2d). We expected dark morph BWPK butterflies to carry the reference allele, pale morphs to be heterozygous, and ivory morphs to be homozygous for a nonreference allele. We expected that at the causative locus, this allelic segregation pattern should be fixed in every individual. This was the case at just 131 SNPs in the whole genome, all within a 9-kb interval on chromosome 15, within the broad association peak (Fig. 2d, blue marks). These SNPs sit within the first intron of the gene cortex. Of particular note, cortex is a switch gene necessary for differentiation into black type II and red type III scales Livraghi et al. 2021), and is a hotspot gene of wing pattern evolution in many species, including in industrial melanism in the peppered moth ( Van't Hof et al. 2016. We noted that in [ivory WT/À] butterflies, scales in red regions exhibited aberrant morphologies (Fig. 2, e and f). Granular clumps of red pigments could be observed in the body of some scales, whereas wild-type scales are solidly pigmented with no granularity. Additionally, we observed scales that were curled at the edges rather than flat. Similar scale phenotypes occur in the wings of cortex crispant butterflies (Fig. 2, e and f), thus supporting a role for cortex loss-of-function in the ivory color phenotypes (Livraghi et al. 2021).

Ivory butterflies carry a large deletion including the cortex promoter
Many phenotypes involve structural genomic variation, both under artificial selection (Bruders et al. 2020), and in natural adaptation (Junker et al. 2020;. As such, we checked read depth across this region in all our sequenced samples. Immediately adjacent to the block of fixed SNPs identified by autozygosity mapping, we found a large region where read depth in [ivory À/À] butterflies dropped to zero, while in [ivory WT/À] it dropped by half. This indicated that a large deletion spans the first exon and promoter of cortex (Fig. 3a).
To determine the precise breakpoints of this deletion, we de novo assembled short-read data for all [ivory WT/WT] and [ivory   À/À] individuals, and BLASTed the resulting contigs against the cortex region. This allowed us to recover individual contigs from some individuals that bridged the 2 ends of the deletion, giving a precise breakpoint at positions Hmel215003o:1494902-1573177, with a length of 78,275 bp relative to the reference genome (Fig. 3b).
The ivory deletion includes 1 of the 2 cortex transcription start sites The ivory deletion contains the only annotated transcription start site (TSS) and promoter for cortex. While butterflies carrying this mutation can complete development, even with effects in scale development across the whole organism and the inability to fly, coding knockouts of cortex (including clonal G 0 crispants) have been observed to have high levels of embryonic lethality, with very few crispant clones observed (Livraghi et al. 2021).
The ability of ivory mutants to survive in spite of the loss of the promoter, as well as the presence of a long first intron in cortex, led us to hypothesize that there may be a second, alternate TSS. In order to determine if this is the case, we mapped RNAseq from H. melpomene for a variety of samples including larval and pupal wings , adult ovaries (Pinharanda et al. 2019), and embryos, adult head and adult abdomen (Dasmahapatra et al. 2012), and looked at splicing and alignment of the 5 0 end of the transcript. No expression was detected in adult head or abdomen, but in embryos, ovaries, and pupal wings, transcription initiates at position Hmel215003o:1424680 in the third annotated exon, which is adjacent to the coding sequence (the "proximal promoter") ( Fig. 4a; Supplementary Figs. 2  and 3). Expression in larval wings initiates at the annotated TSS (the "distal promoter"), and the exon with the proximal promoter is skipped. In pupal wings, transcription again initiated at the proximal promoter. This indicates that cortex has 2 alternative TSSs, a proximal promoter used in multiple tissues and a distal promoter specific to the larval wing. Only the distal promoter is included in the ivory deletion.

Mutagenesis of the distal promoter phenocopies cortex null effects
We reasoned that if the loss of the distal promoter causes the ivory phenotype, knocking out the promoter with CRISPR/Cas9 should phenocopy the effects of the ivory deletion. We generated G 0 mosaic knockouts (crispants) of the promoter in the related species Heliconius erato (Fig. 4, b-d). Resulting butterflies had extensive yellow/white clones on the wings, however, several failed to emerge from the pupa, and a low survival rate and penetrance were observed (Table 1), suggesting that the few wing crispants obtained are rare mosaic escapers of a loss-of-function experiment with lethal effects (Livraghi et al. 2018). Thus, while the pleiotropic effects of deleting the TSS were reminiscent of the previously reported cortex protein coding knockouts, which also exhibited low survival and penetrance (Livraghi et al. 2021), this experiment did not fully recapitulate the ability of the ivory deletion allele to generate large wing clones. This difference in pleiotropy may be due to the use of H. erato in the CRISPR assay, and further experiments are needed to test the functionality of deleted elements within a H. melpomene stock. However, we conclude from these experiments that the disruption of 5 0 distal elements of cortex are necessary for normal color scale patterning, consistent with a causal role of the 78-kb deletion in underlying the ivory phenotypes.

Demographic origins of the ivory allele
The H. m. BWPK stock was kept in insectaries for 30 years (ca. 400 generations) before producing the ivory mutation. The precise history and ancestry of the stock is not recorded, but we do know that the original stock was generated with a mix of butterflies from multiple locations after the addition of butterflies from the vicinity of Tinalandia, Ecuador in the 1980s. The first pale PK (i.e. heterozygotes for ivory) appeared in the PK culture at BW after 2013 and were shared with LEG for study in 2014-2015. We sought to determine if the ivory deletion or associated variation originated from this introduction event from Tinalandia.
Whole genome PCA clustering of a geographic spread of samples of H. melpomene, as well as sister taxa H. cydno and H. timareta, replicates the geographic and species clustering previously described by Martin et al. (2019), with the addition of a distinct "Atlantic" group within H. melpomene (Supplementary Table 1). Consistent with a large contribution from Ecuadorian ancestry, H. m. BWPK cluster together near H. melpomene from the west of the Andes. Similarly, PCA of just the cortex locus, using many more individuals generated by selective sequencing by Moest et al. (2020) (Supplementary Table 4), also places BWPK closest to western H. melpomene.
Given the mixed ancestry of H. m. BWPK, we examined heterogeneity in genome-wide ancestry to determine if variation around the ivory locus in particular was of Ecuadorian origin. We built phylogenetic trees in windows of 50 SNPs across the genome and using topology weighting to determine the closest neighbor of H. m. BWPK at all genomic positions. Of the 15 possible topologies, the most frequent topology places BWPK as an outgroup to the 4 other H. melpomene taxa. The second most common topology groups BWPK with H. melpomene from western Ecuador (H. m. vulcanus and H. m. plesseni), and the third most common places BWPK with H. melpomene from the East of the Andes. In contrast, sliding windows across the cortex locus (including the selective sequencing data) show the most common tree topology here places BWPK sister to the western Ecuador taxon. This result is supportive of H. m. BWPK having an admixed genome, but with the allele at cortex having originated in Ecuador. Importantly, the ivory deletion was not detected in any of 458 wild individuals, indicating it was likely a de novo mutation within the BWPK stock, though it is possible it occurs at very low frequencies in the wild.

Discussion
Using whole-genome sequencing of a set of related individuals, we were able to map the ivory mutation to a large deletion that removed the distal promoter of the wing patterning gene cortex. Cortex was first identified as a switch between melanic and nonmelanic patterns in both Heliconius and the peppered moth Biston betularia, and has later been mapped as a switch between melanic and nonmelanic scale types in a number of other species, including 3 other geometrid moths (Van't Hof et al. 2019), and Papilio butterflies (VanKuren et al. 2019), and also has a role in wing pattern polyphenism in Junonia coenia (van der Burg et al. 2020). This is another example where the same gene has been the target of Fig. 5. Phylogenetic clustering and weighting imply a mixed origin for H. melpomene BWPK. Clustering by PCA from both whole-genome data (a) and from the cortex locus using additional selective sequencing data (b) show similar geographic clustering to that previously reported by Martin et al. (2019); c) with the addition of separation between H. melpomene from East of the Andes and H. melpomene from French Guiana, Suriname, and Brazil (here termed "Atlantic"). Heliconius m. BWPK form a distinct cluster, closest to the West cluster of H. melpomene. Topology weighting was performed using TWISST , with 5 defined groups (West of Andes, Ecuador, East of Andes, Atlantic, and BWPK) (c). d) Genome-wide, the most common topologies were 14 (BWPK as outgroup), 11 (BWPK clustered with East of Andes), and 1 (BWPK clustered with Ecuador). In contrast to this in (e), TWISST results at the cortex locus showed that topology 1 was most common, supporting an Ecuadorian origin for the cortex variants in H. m. BWPK.
natural selection in a wild population and of artificial selection in a captive population.
The list of genes that have been identified as the targets of selection in butterfly wing pattern variation includes transcription factors like Optix and Bab (Reed et al. 2011;Ficarrotta et al. 2021), signaling pathway components like WntA , or terminal effector genes like Yellow , all of which have molecular functions that support a role in developmental patterning or pigment synthesis (McMillan et al. 2020). Cortex, on the other hand, has no clear functional association with color pattern or cell differentiation, and yet has been mapped in a very broad phylogenetic spread of species, suggesting its function may be ancestral in the Lepidoptera. The gene does not appear to have a role in Drosophila wing development or patterning , limiting our ability to make inferences about its molecular function and interactions. In identifying this large deletion, this study provides novel insight into the genetics of this hotspot locus, which will aid future characterization of its developmental function as a switch gene for color patterning.
The loss of all black and red scales in the ivory mutant supports a model in which white/yellow scales are the "default" state during scale cell differentiation, with scale cell precursors failing to differentiate into red or black cell types. This model was initially proposed based on observation and ultrastructure and pigmentation of scales (Gilbert et al. 1987), as well as on the analysis of pattern homologies of the Heliconiini (Nijhout and Wray 1988). It is likely, given the repeated involvement of this locus in wing pattern evolution, that this function of cortex will prove to be conserved throughout Lepidoptera.

Structural variation and mutations
Large deletions have repeatedly been observed in cases of both natural variation and domestication, with 29 deletions affecting regulatory regions which are larger than 1 kb listed on GepheBase (Courtier-Orgogozo et al. 2020). Recurrent deletions of a pitx1 enhancer in sticklebacks, up to 8 kb in length, have caused convergent pelvic reduction (Chan et al. 2010), and the deletion of a 60.7-kb regulatory region at the AR gene in humans, which removes enhancers present in chimpanzee and conserved in other mammals, leads to loss of penile spines (McLean et al. 2011). Similarly large deletions have also been confirmed in the genetic basis of domestication phenotypes, including a 44-kb deletion, also at pitx1, causing feathered feet in chickens (Domyan et al. 2016), and a 141-kb deletion upstream of agouti/ASIP in Japanese quail (Nadeau et al. 2008).
The large deletion that causes ivory is a form of structural variation (SV). Other structural variants affecting this locus have been found in wild populations of Heliconius; multiple inversions have formed a wing pattern supergene in Heliconius numata, and an independent inversion, similar in size and position, has been found in multiple species of the erato clade (Joron et al. 2011;Edelman et al. 2019). It is possible that this genomic region will prove to be prone to a higher-than-average level of SV, or, conversely, that the region is more tolerant to SV than other parts of the genome.

Regulatory consequences of the ivory mutation
Large deletions have the potential to remove a large section of cisregulatory sequence, as well as transcribed sequences. The ivory deletion causes the loss of 1 of 2 promoters. As many as 40% of developmentally expressed genes in Drosophila have 2 promoters which can cause distinct regulatory programs (Bhardwaj et al. 2019), and multiple promoters are also common in human genes (Singer et al. 2008). This both increases the complexity of gene regulatory interactions and increases the number of transcript isoforms per gene-for example, an alternate promoter in the Drosophila gene Zfh1 creates an isoform that has a shorter 5 0 UTR which is missing miRNA seed sites on its 5 0 UTR, permitting differential degradation of the mRNA by miR-8 (Boukhatmi and Bray 2018). We determined that the ivory deletion contains 1 annotated promoter of cortex, but that in some contexts during development, another TSS is used, much closer to the translation start site. Alternate promoter usage is likely to be common and widespread in animals.
We expect that the use of 2 separate, context-specific promoters at cortex will create a hierarchical regulatory organization, where CREs used in different tissues or at different times can loop to different promoters. Additionally, transcripts from each of the 2 promoters contain a different arrangement of 5 0 noncoding exons, giving them different UTRs, possibly leading to differential translational regulation or degradation of the mRNA.
The ivory deletion does not include protein-coding sequence, which led us to hypothesize that by removing just 1 promoter, cortex would still be expressed in tissues that use the alternate promoter, and that therefore ivory mutants would bypass the pleiotropic, highly lethal effects observed in CRISPR experiments that target coding exons of cortex (Livraghi et al. 2021). However, crispants with deletions in the distal promoter in H. erato did not appear to be free of these pleiotropic effects, and exhibited high lethality and comparatively small clone sizes, dissimilar to the H. melpomene [ivory À/ À] mutants. This may be due to inherent differences in the regulatory function of the distal TSS between these 2 species, or to the existence of other functional sequences within the deletion which might be necessary for generating the pale phenotypes. Specifically, the 78-kb deletion removes part of the 3 0 UTR of the neighboring gene parn as well as 2 miRNAs (Surridge et al. 2011), and it is also possible that a large deletion affects other cortex cis-regulatory elements or affect 3D chromatin interactions on a large portion of the chromosome. Additionally, the ivory butterflies carry a set of fixed SNPs adjacent to the deletion, which could contribute to the phenotype. In summary, the complete loss of black and red scales caused by the ivory deletion and by CRISPR-Cas9 mutagenesis in 2 Heliconius species indicates that either the distal promoter or another functional sequence within the 5 0 noncoding region of cortex are necessary for black and red scale development.

Heterozygote advantage of deleterious mutations in captivity
The cortex mutant phenotype has been artificially selected by a breeder in the heterozygous state, when unusually widespread yellow and white patterns emerged in the stock. The heterozygote carriers of the ivory mutation ([ivory WT/À]) do not appear to have any reduction in fitness or vigor. In contrast, homozygous ivory butterflies of either sex are unable to fly (Supplementary Movie 1), limiting their ability to feed and reproduce-we have not observed any successful matings involving an [ivory À/À] butterfly. With a homozygous fitness of zero, those recessive effects make the ivory mutant allele functionally equivalent to a recessive lethal allele with a heterozygous phenotype under artificial selection such as in PAX3 overo horses (28). This was previously formalized by Hedrick with models for lethal and "near-lethal" mutations with heterozygous advantage, and shown to be only possible under conditions of intense selective breeding (Hedrick 2015). To further probe the empirical evidence for this phenomenon, we compiled cases of heterozygote advantage under artificial selection from the literature using a database of gene-tophenotype relationships (Courtier-Orgogozo et al. 2020), and found 39 analogous gene-to-trait relationships spanning domesticated mammals and birds ( Table 2), 15 of which are recessive- lethal. Large SVs such as the ivory deletion accounted for 12 out of 83 derived alleles in this data set, suggesting that the recessive deleterious effects of macromutations occasionally provide heterozygous states of interest to artificial selection. Of note, 17 out the 40 cases of captive heterozygote advantage involve selection for depigmentation traits, highlighting the trend among breeders and fanciers to select for conspicuous variants (Table 2). Finally, cases of heterozygous advantage also occur in the wild, such as in ruff birds where a large inversion haplotype provides male coloration and behavioral traits that are maintained by sexual selection in spite of being recessive lethal (Kü pper et al. 2016). In wild populations of the polymorphic butterfly H. numata, cortex itself is situated in the center of an inversion polymorphism under balancing selection (Maisonneuve et al. 2021). Further work will be needed to determine the subgene level elements that are within the ivory deletion and mediate homozygous inviability.

Conclusion
By using autozygosity mapping and association on a comparatively small pool of individuals, we were able to identify a structural mutation involved with butterfly wing pattern. This approach may prove fruitful in other studies of butterfly wing patterning or in the identification of de novo mutants, especially if combined with recent developments in whole-genome sequencing from dried museum specimens (Cong et al. 2021;Grewe et al. 2021). This could allow the mapping of other cases like Hindsight in J. coenia, or pseudozorro in Parnassius apollo (Weatherbee et al. 1999;Pierrat and Descimon 2011). The identification and characterization of spontaneous mutants has provided very valuable insights into the genetics of development in model organisms. Further whole-genome sequencing of deleterious mutations in butterflies will help to identify parts of the genome that are functionally required for normal development, which will assist in a more complete understanding of their evolution and development.

Data availability
Whole-genome sequences are accessible on the NCBI SRA under the Bioproject accession number PRJNA610063.
Supplemental material is available at G3 online.