A large deletion at the cortex locus eliminates butterfly wing patterning.

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 ; Cieslak ). 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 ; Giovannini ), melanin-based coat color in domestic mammals (Cieslak ), plumage patterns, and colors in birds (Domyan and Shapiro 2017; Price-Waldman and Stoddard 2021), and chromatophore distribution in squamates and fishes (Guo ). These variants can affect genes involved in the enzymatic production, transport, and deposition of pigments such as the melanin and carotenoid pathways (Mundy ; Courtier-Orgogozo and Martin 2020), or can be caused by genes that affect signaling or cell type differentiation [as with chromatophore distribution in squamates (Kuriyama )]. 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 ), and a wild-type regulatory variant in wall lizards causes changes to ventral scale color (Andrade ). The gene Agouti has repeatedly been mapped in cases of pigment differences, both in cases of domestication [e.g. in the chicken (Yu ), horse (Rieder ), and rabbit (Fontanesi )], and in natural variation [e.g. in warblers (Toews ), humans (Bonilla ), and snowshoe hare (Jones )]. 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, cis-regulatory mutations, and copy number variation (Norris and Whan 2008; Gratten ). 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 ; Makino ; Moyers ). 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 ; Langevin ), and in overo horses, where foals homozygous for a mutated Endothelin Receptor B (EDNRB) develop Lethal White Syndrome (Metallinos ). “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 ). 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 ), 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 ). 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 ).

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 (Jiggins ; Van Belleghem ), 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 out-crossing, 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).

Fig. 1.

Phenotypes of H. melpomene BWPK butterflies. a) Ivory phenotypes in the H. melpomene “Piano Keys” in the UT Austin stock (BWPK). The 3 color states (from left to right “Dark BWPK,” “Pale BWPK,” and “ivory”) depend on the allelic dosage of a codominant mutation. Ivory homozygotes are inviable and only found among offspring from 2 pale heterozygotes. b) Magnified view of the forewing red band region. In the [ivory WT/−] state, abnormal scales have formed in the red region, with red pigment in granules and scales curled, while in [ivory −/−], all scales are yellow or white. c) Magnified view of a central hindwing region that is black in [ivory WT/WT]. The [ivory WT/−] wing has all cover scales as yellow or white while all ground scales remain black. All scales are yellow or white in [ivory −/−]. Complete replacement of melanic scales by yellow-white scales in the abdomen (d) and legs (e) from [ivory −/−] homozygotes. f) Cartoon of image locations for (b–e). Scale bars: b–c = 50 µm; d = 500 µm; e = 100 µm.

The genetics of wing patterning have been extensively studied in Heliconius (McMillan ). 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 (Wallbank ; Zhang ; Lewis ; Morris ;), the signaling ligand WntA (Mazo-Vargas ; Van Belleghem ), and cortex, a cdc20 homolog with a currently unidentified function (Nadeau ; Livraghi ). 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.

Methods

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, RNAse-treated, 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 ) retrieved from LepBase (Challis ) with BWA-MEM using default parameters (Li 2013), and variants called using GATK v4.1 with tools HaplotypeCaller and GenotypeGVCFs using default parameters (McKenna ). 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 ). SNP association was carried out in PLINK v.1.9, with 1,000 permutations (Purcell ). For phylogenetic analyses, we followed the TWISST pipeline(Martin and Van Belleghem 2017); 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 ) 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 ).

CRISPR/Cas9 mutagenesis

We designed sgRNAs against the putative H. erato cortex promoter corresponding to N20NGG sites within the distal promoter/5′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/µl, respectively. Embryonic injections were performed as previously described (Livraghi ) within 1–3 h after egg laying, after which larvae were allowed to develop on a diet of P. biflora until adult emergence.

Table 1.

CRISPR mutagenesis of the cortex distal promoter/5′UTR in H. erato.

Batch	sgRNA target	Injected	Pupae	Adults	Crispants	Survival	Penetrance
1	sgProm1: 5′AGCGGTAACTTTATCGCGAT	512	42	37	3	7.20%	0.6%
2	sgProm2: 5′ATGGGATATGTTTTAATAGT sgProm3: 5′AACCTGTCACATCAGTACAG	695	78	73	1	11.2%	0.1%

Results

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 wild-type 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 ) (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.

Fig. 2.

Genome-wide association and autozygosity mapping in H. melpomene BWPK. a) Cartoons of the allelic variation present in the H. m. BWPK stock. b) The hindwing veins, showing the composite effects of different alleles depicted in (a). c) GWAS for the Dennis element; Manhattan plot of Wald test P-values, with fixed SNPs in red, and a magnified annotation of the region around the gene optix. d) GWAS for ivory, with fixed SNPs in blue and a magnified annotation of the cortex region, with the point of the arrow at the annotated TSS. e) Comparison of red scales in [ivory WT/WT], [ivory −/−], and cortex crispant wings. f) Some ivory butterflies have a small red dot at the base of the ventral hindwing, the only scales on the entire butterfly that are not white or yellow. They share the atypical phenotype of red scales from cortex crispants (e).

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 (Nadeau ; Livraghi ), and is a hotspot gene of wing pattern evolution in many species, including in industrial melanism in the peppered moth (Van’t Hof , 2019).

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 ).

Ivory butterflies carry a large deletion including the cortex promoter

Many phenotypes involve structural genomic variation, both under artificial selection (Bruders ), and in natural adaptation (Junker ; Mérot ). 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).

Fig. 3.

Ivory is associated with a large deletion in the gene cortex (a) mean coverage plot for the region around cortex, with [ivory WT/WT] in green, [ivory WT/−] in red, and [ivory −/−] in blue. b) Alignments of de novo scaffolds, indicating the precise breakpoints of the ivory deletion. Colors indicate the orientation of sequences relative to the reference scaffold depicted in (a).

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 G0 crispants) have been observed to have high levels of embryonic lethality, with very few crispant clones observed (Livraghi ).

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 (Hanly ), adult ovaries (Pinharanda ), and embryos, adult head and adult abdomen (Dasmahapatra ), and looked at splicing and alignment of the 5′ 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.

Fig. 4.

Knockout of the distal promoter phenocopies cortex protein-coding knockouts. a) Mapping of RNAseq reads at cortex suggests that there are 2 promoters and TSS. Numbers of intron-spanning reads are indicated (see Sashimi plots in Supplementary Figs. 2 and 3). Knockout of the distal promoter in H. erato causes a transformation of black scales to white/yellow in a phenocopy of cortex protein coding knockouts. b) A butterfly with extensive clones but that emerged poorly, indicating pleiotropic effects not observed in the ivory mutants, while c) much smaller clones similar to those reported by Livraghi , with clone positions indicated by pink arrows. d) High magnification images of mutant clones. Scale bars: 100 µm.

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 ; 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 (Martin and Van Belleghem 2017), 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.

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 G0 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 ). 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 ), 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′ 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 , 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 (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 ), and Papilio butterflies (VanKuren ), and also has a role in wing pattern polyphenism in Junonia coenia (van der Burg ). This is another example where the same gene has been the target of 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 ; Ficarrotta ), signaling pathway components like WntA (Martin ), or terminal effector genes like Yellow (Martin ), all of which have molecular functions that support a role in developmental patterning or pigment synthesis (McMillan ). 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 (Nadeau ), 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 ), 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 ). Recurrent deletions of a pitx1 enhancer in sticklebacks, up to 8 kb in length, have caused convergent pelvic reduction (Chan ), 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 ). 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 ), and a 141-kb deletion upstream of agouti/ASIP in Japanese quail (Nadeau ).

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 ; Edelman ). 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 cis-regulatory 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 ), and multiple promoters are also common in human genes (Singer ). 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′UTR which is missing miRNA seed sites on its 5′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′ 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 ). 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′UTR of the neighboring gene parn as well as 2 miRNAs (Surridge ), 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′ 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-to-phenotype relationships (Courtier-Orgogozo ), 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 ). In wild populations of the polymorphic butterfly H. numata, cortex itself is situated in the center of an inversion polymorphism under balancing selection (Maisonneuve ). Further work will be needed to determine the subgene level elements that are within the ivory deletion and mediate homozygous inviability.

Table 2.

Known cases of heterozygote advantage in animal breeding (see Gephebase for references).

Gene	Species	Heterozygote advantage (HA) trait for breeding	Detrimental recessive effect	Nb of HA alleles	Mutation types
ACAN	Cattle	Short stature	Lethal	1	Frameshift (4-bp ins)
				1	New start codon (SNP)
ACAN	Horse	Short stature	Lethal	1	Loss of 6 a.a. (18-bp del)
				1	Frameshift (1-bp del)
ALX1	Cat	Brachycephaly	Frontonasal dysplasia	1	Loss of 4 a.a. (12-bp del)
BBS9	Pig	Fast growth	Lethal (misexpressed BMPER)	1	Large SV (212-kb del)
BMP15	Sheep	Increased female fecundity	Female infertility	5	Missense (SNP)
				2	Nonsense (SNP)
				1	Frameshift (1-bp ins)
				1	Frameshift (17-bp del)
				1	Cis-regulatory (SNP)
cortex b	Heliconius	White patterns	Flightless	1	Large SV (78-kb del)
EDNRB	Horse	White patterns	Lethal	1	Missense (SNP)
FGF3/4/19 cluster	Dog	Dorsal hair ridge	Dermoid sinuses	1	Large SV (133-kb duplication)
FGFR3	Sheep	Increased bone length	Chondrodysplasia	1	Missense (SNP)
FOXI3	Dog	Hairlessness	Ectodermal dysplasia	1	Frameshift (7-bp duplication)
GDF9	Sheep	Increased female fecundity	Female infertility	5	Missense (SNP)
HMGA2	Rabbit	Short stature	Lethal	1	Large SV (12.1-kb del)
HOXC10	Chicken	Head crest	Cerebral Hernia—likely epistatic	1	Large SV (197-bp duplication)
IHH	Chicken	Short legs	Lethal	1	Large SV (11.9-kb del)
KIT	Camel	White patterns	Lethal	1	Frameshift (1-bp del)
KIT	Dog	White patterns	Lethal	1	Frameshift (1-bp ins)
KIT	Donkey	White patterns	Lethal—not proven	1 1	Missense (SNP) Splice Site mutation (SNP)
KIT	Fox	White patterns	Lethal	1 1	Missense (SNP) Splice Site mutation (SNP)
KIT	Horse	White patterns	Lethal	3	Frameshift (1-bp, 4-bp del)
				2	Large SV (1.3-kb, 1.9-kb del)
				2	Nonsense (SNP)
				1	In-Frame Del (54-bp)
				1	Splice Site mutation (SNP)
KIT	Pig	White patterns	Lethal—not proven	1	Splice Site mutation (SNP)
MITF	Buffalo	White patterns	Microphtalmia (inferred from mice)	1 1	Nonsense (SNP) Splice Site mutation (SNP)
MITF	Dog	White patterns	Deafness (risk increasea)	1	Cis-regulatory (195-bp ins)
MITF	Horse	White patterns	Microphtalmia (inferred from mice)	2	Large SV (8.7-kb, 63-kb del)
				1	Frameshift (4-bp del)
				1	Missense (SNP)
MITF	Quail	White patterns	Small size and slow growth	1	Frameshift (2-bp del)
Mlana	Pigeon	White patterns	Eye defects	4	CNV including 3 other genes
MNR2 locus	Chicken	Head comb	Male infertility	1	Large SV (7.4-Mb inversion)
MRC2	Cattle	High muscle	Musculoskeletal defects	1	Frameshift (2-bp del)
				1	Missense (SNP)
MSTN c ,d	Dog	Racing performance	Nonconform appearance	1	Frameshift (2-bp del)
PAX3	Horse	White patterns	Lethal	2	Missense (SNP)
PMEL17	Dog	White patterns	Auditory and ocular defects	6	Splice Site mutation (TE insertions)
PMEL17	Horse	White patterns	Ocular defects	1	Missense (SNP)
RNASEH2Ba locus	Cattle	Milk yield (causal gene unknown)	Lethal (due to RNASEH2Ba KO)	1	Large SV (660-kb deletion)
RYR1	Pig	Meat yield and quality	Hyperthermia	1	Missense (SNP)
PRLR/SPEF2	Pig	Increased female fecundity	Sperm defects (due to SPEF2 KO)	1	Cis-regulatory (TE insertion)
SPINT1	Gecko	White patterns	Metastases	1	Unknown
STX17	Horse	White patterns	Melanoma	1	Cis-regulatory (4.6-kb duplication)
TBXT (Brachyury)	Cat	Short or no tail	Lethal	4	Frameshift (1-bp del; 14-bp ins)
TBXT (Brachyury)	Dog	Short or no tail	Lethal	1	Missense (SNP)
TRPM1	Horse	White patterns	Ocular defects	1	Cis-regulatory (TE insertion)
TRPV4	Cat	Short ears	Chondrodysplasia	1	Missense (SNP)

The corresponding entries and references can be retrieved using a search for “@HeterozygoteAdvantage” or for the gene and species names at www.gephebase.org.

This study.

Recommendations to outcross are recent in the affected stock.

MSTN loss-of-function alleles exist in other bred species. Homozygotes sometimes requires assisted birthing, but are viable and artificially selected.

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 ; Grewe ). This could allow the mapping of other cases like Hindsight in J. coenia, or pseudozorro in Parnassius apollo (Weatherbee ; 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.

Click here for additional data file.

Click here for additional data file.