An important goal in evolutionary biology is to understand the genetic changes underlying novel morphological structures. We investigated the origins of a complex wing pattern found among Amazonian Heliconius butterflies. Genome sequence data from 142 individuals across 17 species identified narrow regions associated with two distinct red colour pattern elements, dennis and ray. We hypothesise that these modules in non-coding sequence represent distinct cis-regulatory loci that control expression of the transcription factor optix, which in turn controls red pattern variation across Heliconius. Phylogenetic analysis of the two elements demonstrated that they have distinct evolutionary histories and that novel adaptive morphological variation was created by shuffling these cis-regulatory modules through recombination between divergent lineages. In addition, recombination of modules into different combinations within species further contributes to diversity. Analysis of the timing of diversification in these two regions supports the hypothesis of introgression moving regulatory modules between species, rather than shared ancestral variation. The dennis phenotype introgressed into Heliconius melpomene at about the same time that ray originated in this group, while ray introgressed back into H. elevatus much more recently. We show that shuffling of existing enhancer elements both within and between species provides a mechanism for rapid diversification and generation of novel morphological combinations during adaptive radiation.

Butterflies show an amazing diversity of patterns on their wings. In fact, most of the 18,000 species of butterfly can be distinguished on the basis of their wing pattern. Much of this diversity is thought to arise through novel switches in the genome that turn genes on in new contexts during wing development, thereby producing new patterns. Here we study a set of switches that control the expression of optix, a gene that places red patches onto the wings of Heliconius butterflies. We show that two patterning switches—one that produces red rays on the hindwing and the other a red patch on the base of the forewing—are located adjacent to one another in the genome. These switches have each evolved just once among a group of 16 species but have then been repeatedly shared between species by hybridisation and introgression. Despite the fact that they are now part of a common pattern in the Amazon basin, these two pattern components actually arose in completely different species before being brought together through hybridisation. In addition, recombination among these switches has produced new combinations of patterns within species. Such sharing of genetic variation is one way in which mimicry can evolve, whereby patterns are shared between species to send a common signal to predators. Our work suggests a new mechanism for generating evolutionary novelty, by shuffling these genetic switches among lineages and within species.

Funding: This study was funded by the Biotechnology and Biological Sciences Research Council http://www.bbsrc.ac.uk grant number H01439X/1 to CDJ; the European Research Council MimEvol grant to MJ and the European Research Council Speciation Genetics to CDJ. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Here we focus on the H. melpomene lineage, in which the Amazonian dennis-ray phenotype has evolved recently from a red-banded ancestor [ 31 ]. We carry out a population genomic analysis on H. melpomene and its relatives, H. elevatus and H. timareta, to identify putative regulatory modules associated with distinct red pattern elements. Previously, population genetic evidence has suggested that mimicry among H. melpomene, H. elevatus, and H. timareta has evolved through sharing of the dennis-ray allele by repeated adaptive introgression at the optix locus [ 16 ]. This is especially surprising in the case of H. elevatus, which forms part of the “silvaniform” clade that diverged from H. melpomene around 4 million years ago [ 14 ]. Our analysis here indicates that the origin of the red pattern elements is considerably more complex than has been previously supposed, with the dennis and ray elements of the widespread dennis-ray pattern having distinct evolutionary origins in different clades within the genus.

(A) Examples of the principal red wing patterns: dennis-ray H. e. pseudocupidineus and H. m. malleti, dennis-only H. m. meriana, ray-only H. t. timareta f. contigua, and band H. m. rosina. (B) Association analysis across 96 genomes showing statistical association for the dennis (red dots) and ray phenotypes (orange dots). The horizontal black line represents significance after Bonferroni correction. Boxes show exon positions. (C) Sliding window analysis of fixed differences between specific comparisons to identify dennis- and ray-associated sites (orange and red lines—see also S1 Fig and S1 Table ). (D) Recombination breakpoints allowed separate isolation of regions fully associated with dennis and ray phenotypes. Informative haplotypes are shown (H. elevatus, H. m. malleti, H. t. timareta f. contigua, H. m. meriana, and H. m. rosina, phenotypes shown above). Genotypes are indicated as D/d for dennis present/absent and R/r for ray present/absent. A H. melpomene recombinant hybrid was heterozygous in the dennis region and homozygous for ray-absent, as expected, but was not informative for precise breakpoint delineation because of missing data. See Dryad depository for plot data [ 32 ].

Generally, the patterns on butterfly wings are a good system in which to link genetic changes to the developmental processes that generate diversity [ 18 , 19 ]. Wing colour patterns are mosaics of scales, each with a single colour, produced by a combination of pigment and ultrastructure. The relative positions of differently coloured scales are established during larval and pupal wing development [ 20 ]. Wing development is thought to be broadly conserved in insects, with wing developmental genes showing similar expression patterns between flies and butterflies [ 21 – 23 ]. This therefore raises the question: how is this conserved landscape of wing development translated into the diversity of butterfly wing patterns? In Heliconius, pattern diversity is controlled by a surprisingly small number of genomic regions with large effect sizes [ 24 , 25 ]. In particular, genetic mapping and gene expression studies have shown that red elements are associated with expression of the transcription factor optix across all Heliconius species [ 26 , 27 ]. In the absence of fixed coding sequence changes between wing pattern forms, this implies that red pattern variation is controlled by differential regulatory control of optix [ 27 ]. Population genomic studies have identified a region of non-coding sequence downstream of optix that is associated with phenotypic change [ 15 , 28 ]. Previous work suggests that there may be several distinct elements within this region. Occasional hybrid phenotypes possess only the “dennis” patch on the base of the forewing or the “ray” elements on the hindwing and have been hypothesised to be rare recombinants, although this has never been tested genetically [ 29 , 30 ]. Similarly, there are also established forms that exhibit only dennis or ray patterns (H. melpomene meriana and H. timareta timareta f. contigua, respectively, Fig 2 ). This suggests that the broad genomic interval already identified might contain discrete regulatory loci that vary the spatial expression of optix in different wing regions, a hypothesis that we can now test with genetic data.

First row: H. burneyi huebneri, H. aoede auca, and H. xanthocles zamora; second row: H. timareta timareta f. timareta, H. doris doris, and H. demeter ucayalensis; third row: H. melpomene malleti, H. egeria homogena, and H. erato emma; fourth row: H. elevatus pseudocupidineus, Eueides heliconioides eanes, and E. tales calathus; and bottom: Chetone phyleis, a pericopine moth. Stars indicate the three species that are the focus of this study. Butterflies figured are from the Neukirchen Collection, McGuire Centre, Florida. The butterflies are from populations in both Ecuador and Peru.

Here we explore the origins of adaptive novelty among the wing patterns of Heliconius butterflies. These wing patterns are under strong natural selection for mimicry and warning colour, as well as being important mating signals [ 13 ]. The rapid radiation in Heliconius is accompanied by an even more rapid diversification in mimicry patterns as well as convergence among species found in a given locality [ 14 ], both through independent convergent evolution and via introgression of gene regions between races and species [ 15 , 16 ]. Mimetic convergence reaches its peak among red dennis-ray pattern phenotypes in the Amazon ( Fig 1 ), where 11 or more Heliconius species, as well as pierine butterflies and pericopine moths, share the same pattern. In addition to near perfect convergence in wing patterns in a given locality, there is also often striking divergence of patterns between localities as populations adapt to the many different mimicry complexes spread across the Neotropics [ 17 ]. This diversity provides an opportunity to study the genetic and developmental basis of evolutionary novelty.

Much of our understanding of modularity in regulatory evolution comes from Drosophila, in which the loss of trichomes on the larval cuticle [ 5 ], the gain of melanic wing spots [ 8 – 10 ], or changes in abdominal pigmentation [ 3 , 11 ] have been shown to involve evolutionary changes in cis-regulatory elements. These elegant developmental studies demonstrate the underlying logic of regulatory modularity, whereby novel expression domains can arise without disrupting existing function. These studies have also established a paradigm in which small effect mutations alter transcription factor binding sites in these regulatory modules and in combination produce large effect alleles [ 5 ]. Similar conclusions come from recent work in other taxa, including mice and jewel wasps [ 2 , 12 ]. This might seem to imply that the evolution of novel regulatory alleles is relatively gradual, requiring the evolution of many small effect substitutions, but recent adaptive radiations can show extremely rapid rates of morphological change. The role of regulatory modularity therefore remains to be tested in adaptive radiations in which morphological variation evolves very rapidly.

One of the major impediments to evolutionary innovation is the constraint on genetic change imposed by existing function [ 1 ]. Mutations that confer advantageous phenotypic effects in a novel trait will often result in negative pleiotropic effects in other traits influenced by the same gene. Several mechanisms have been proposed by which evolution can circumvent such constraints, resulting in phenotypic diversification. In particular, the modularity of cis-regulatory elements [ 2 – 6 ] means that novel modules can encode new expression domains and functions without disrupting existing expression patterns [ 6 , 7 ]. This modularity underlying gene regulation has led to the assertion that much of morphological diversity has arisen through regulatory evolution [ 6 ].

Results and Discussion

We took advantage of natural phenotypic variants in which the two red elements, dennis and ray, occur separately to identify putative functional regulatory regions controlling red pattern within the H. melpomene clade. Genomic analysis of 96 individuals from the melpomene-timareta clade revealed two distinct regions that showed strong association with the dennis and ray pattern elements, respectively. Our analysis included a race of H. melpomene, H. m. meriana, from the Guiana shield, which possesses the forewing dennis patch but not ray, as well as H. t. timareta f. contigua from Ecuador, which possesses ray but not dennis, plus a recombinant individual from an H. melpomene hybrid zone in Ecuador with dennis but not ray (Fig 2A). Across all 96 individuals, there were significant genotype-by-phenotype associations across all genome regions surveyed. This “background” signal of genotype-by-phenotype association is likely due to the presence of genetically divergent species in our dataset that are to some degree confounded with phenotype. Nonetheless, our analysis identified a peak of genotype-by-phenotype association spanning roughly 50 kb and located from 60–110 kb 3ʹ of the optix gene, similar to what has been observed previously (Fig 2B) [28]. This region also corresponds closely to that identified recently in the mimetic species H. erato [15], implying convergence in the regulatory architecture controlling wing pattern mimicry at a finer scale than has been previously demonstrated [28,31].

Furthermore, within this region in our data, distinct adjacent peaks of association were observed for the dennis and ray elements. Focusing specifically on fixed single nucleotide polymorphism (SNP) differences between alternative red phenotypes revealed two distinct peaks of association (Fig 2C). One, approximately 10 kb in length, contained SNPs perfectly associated with the red dennis patch. The other adjacent region was broader, roughly 25 kb, and contained SNPs perfectly associated with red hindwing rays.

We next used broader taxonomic sampling to further refine these intervals and identify exact sequence haplotypes associated with each of the two phenotypic elements (Fig 2D). To identify recombination breakpoints around dennis and ray haplotypes, we generated a high-quality sequence alignment by de novo assembly of each individual genome and then identified contigs across the associated region using the Basic Local Alignment Search Tool (BLAST). For the dennis region, alignment was assisted by a sequenced fosmid clone from H. m. aglaope (dennis phenotype) to complement the reference genome (derived from a non-dennis butterfly). The final alignment included the 96 melpomene-timareta individuals used for association analysis and a further 46 individuals that included species with no red (H. cydno) and species from the more distantly related silvaniform clade including H. elevatus, which has the dennis-ray pattern. The distal end of the dennis region, relative to optix, was delineated by a rapid loss of phenotype-associated variants across all species sampled, whilst the proximal end was determined by a single fixed recombination event in the race H. m. meriana (dennis but no ray phenotype), generating a region of ~7 kb fully associated with dennis. For ray, a breakpoint in the ray-only H. t. timareta f. contigua defined the distal end, whilst a recombination in H. m. meriana defined the proximal end, resulting in a larger ~37 kb region (Fig 2D). Each haplotype group was characterised by diagnostic SNPs as well as a fixed architecture of indel variation (Fig 3). These analyses therefore support the hypothesis derived from phenotypic evidence, that dennis and ray phenotypes are controlled by adjacent distinct genetic elements. In combination with previous work showing differential expression of optix across a wide diversity of Heliconius species and races, this provides clear genetic evidence for modularity in the cis-regulatory control of optix.

PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 3. Schematic overview of the dennis region alignment architecture. Each group of alleles is characterised by a complex structure of insertion and deletion variation. Horizontal bars represent aligned allele sequences for the outgroups (A), silvaniforms (B), H. elevatus (C), dennis morph melpomene/timareta (D), non-dennis morph melpomene/timareta (E), and the cydno/heurippa/pachinus clade (F). There are three general forms of the region: an “outgroup” allele (X) with many indels and lacking large sections of sequence; a silvaniform/dennis allele (Y), which possesses indel 1 and lacks indel 2; and a non-dennis/cydno allele (Z), which lacks indel 1 and possesses indel 2. Dotted lines and blue colouration indicate the point where dennis-morph melpomene/timareta alleles become silvaniform like (left) and where the dennis-only morph, H. m. meriana, recombines with the other non-rayed melpomene (right). Fixed SNPs in perfect association with the dennis phenotype are indicated by yellow asterisks and form two clusters. See Dryad depository for alignment [32]. https://doi.org/10.1371/journal.pbio.1002353.g003

We have previously hypothesised that the dennis-ray mimicry pattern introgressed as a single genomic block between H. melpomene and H. timareta, as well as more distantly between H. melpomene and H. elevatus [16]. Our new data suggest a much more complex history than previously recognised, with dennis and ray having quite distinct origins. As expected, a maximum likelihood (ML) phylogeny shows that the ray alleles fall within the H. melpomene clade, indicating an origin derived from an ancestral H. melpomene phenotype. In contrast, however, alleles producing the dennis phenotype originated within the silvaniform clade, which diverged from H. melpomene around 4 million years ago (Fig 4) [14]. Members of this clade have mottled orange/red, black, and yellow “tiger” patterns and are mostly co-mimics of butterflies in the tribe Ithomiini, whereas the melpomene-cydno clade are all co-mimics of other Heliconius species. Nonetheless, the silvaniforms commonly have orange patches on the base of the forewing, which in some cases are remarkably similar to the dennis patch of H. melpomene. In particular, the form H. hecale metellus has a dennis-like phenotype (Fig 4A), which suggests a plausible ancestral phenotype that might have provided the source of the dennis allele in H. melpomene.

PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 4. Modularity of dennis and ray and phylogenetic relationships among alleles at both loci. (A) The phylogeny of species used in this study (left, from [14]) and their respective combinations of dennis and ray haplotypes (right). Note that although genetic data for Heliconius elevatus roraima are not analysed here or included in our alignment, this taxon is shown for completeness. Representative butterfly phenotypes are shown with their respective subspecies or form names. H. e. pseudocupidineus is abbreviated to pseudoc. (B) ML trees with bootstraps (generated using de novo assembled genomes) show that the dennis alleles cluster with the silvaniform clade (left), whilst the ray alleles cluster with the melpomene/timareta clade (right). Terminal nodes are coloured by silvaniform clade (green, except H. elevatus), dennis/ray (red, H. melpomene/timareta; burgundy, H. elevatus), melpomene/timareta non-dennis/ray forms (blue), and cydno/pachinus/heurippa clade (purple). Outgroup species are in black (see also S2 Fig for sample labels). H. m. meriana and H. t. timareta f. contigua have dennis-only and ray-only patterns, respectively, and cluster with their expected phenotypes. All trees were rooted to H. aoede. Bootstraps are given as percentages of 1,000 iterations. https://doi.org/10.1371/journal.pbio.1002353.g004

Various scenarios might explain this complex history. Sharing of variation between species can be explained by either retention of ancestral polymorphism or introgression through hybridisation. We can directly test these alternative scenarios using dated trees inferred from our alignments. In order to provide comparable trees, we used the divergence date between the silvaniform and melpomene-cydno clades derived from a recent species tree for the Heliconinii [14] to calibrate the dennis and ray region phylogenies. These trees support introgression and rule out ancestral polymorphism because the dates of coalescence of the H. melpomene and H. elevatus dennis and ray alleles are significantly more recent than the divergence of these two species. These species last shared a common ancestor at around 3.96 Ma (95% highest posterior density [HPD] interval 3.18–4.81 Ma) [14]. In contrast, the dennis allele shared between H. elevatus and H. melpomene/timareta diverged around 1.95 Ma (2.79–1.25 Ma HPD). The divergence of the ray allele is even more recent and shared a common ancestor between H. elevatus and H. melpomene/timareta around 0.66 Ma (0.93–0.43 Ma HPD). The recent origin of these alleles is also supported by low levels of genetic diversity within these clades, with the average pairwise sequence divergence among the dennis alleles only 1.5%, including those from H. melpomene, H. timareta, and the more distantly related H. elevatus (Table 1, top). This is less than that found among the same individuals in flanking sequence (2.4%) and comparable to that among the red-forewing-banded “postman” group for the same locus (1.6%), which includes only more closely related melpomene and timareta individuals. The ray alleles also show only 1.1% average pairwise sequence divergence at the ray locus, similarly less than in the postman group at the same locus. Although sequence diversity is likely to be reduced in these regions because of functional constraint, it seems likely that such constraint is similar across different clades in the phylogeny, so the relatively low levels of diversity within the ray and dennis clade support their recent origin.

Our dated trees can also be used to infer the relative timing of introgression events. Here the data indicate that the ray allele originated within H. melpomene at around the same time as dennis introgressed into the H. melpomene clade from an ancestor of H. elevatus, sometime around 1.85 Ma (2.53–1.25 Ma HPD; Fig 5 and S3 Fig). This suggests that the characteristic mimetic dennis-ray phenotype first came together within H. melpomene at that time. In contrast, H. elevatus did not acquire the ray allele until about a million years later and perhaps persisted during this time as part of the Guiana Shield dennis-only mimicry ring. The dennis and ray alleles of H. timareta are each nested within those of H. melpomene and are more recent than the divergence of these species, implying introgression from H. melpomene into H. timareta and consistent with previous analyses [33]. Nonetheless, dennis and ray events also differ in timing, as most of the H. timareta ray alleles diverge from their H. melpomene relatives around 1 Ma, but the dennis alleles diverged only around 0.45 Ma (Fig 5 and S3 Fig). H. timareta ray alleles are polyphyletic with respect to H. melpomene, also supporting multiple introgression events and recombination between the regions. Further sampling will be needed to resolve more clearly the timing and number of introgression events between these species.

PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 5. Hypothesis for the origins and introgression of the dennis and ray regions inferred from dated trees. Key events indicated on the x-axis are (A) introgression of dennis from H. elevatus into H. melpomene, inferred from the coalescence time of alleles sampled from these two species; (B) the origin of ray alleles within H. melpomene, inferred from the coalescence of dennis and non-dennis alleles within H. melpomene; (C) introgression of ray into H. timareta; (D) introgression of ray into H. elevatus, and (E) introgression of dennis into H. timareta. Dating of phylogenies was carried out using BEAST, and full dated trees with support limits are provided in S3 Fig. The species tree is derived from [14]. https://doi.org/10.1371/journal.pbio.1002353.g005

In addition to recombination between lineages, there is also shuffling of alternate alleles at these regulatory modules within species. Across most of their range, H. melpomene and H. timareta have either postman or dennis-ray haplotypes across the entire studied region. However, the race H. m. meriana has dennis alleles but shows recombination across the adjacent ray locus that removes the ray phenotype. Similarly, a single recombination event in H. t. timareta f. contigua produced a phenotype with ray but not dennis (Fig 2D; Fig 4). H. elevatus, like H. melpomene, also has a dennis-only race found in the Guiana Shield, which likely represents another case of enhancer shuffling within species, although we have not sampled this species here (Fig 4). Hence, although alleles within these regulatory modules are now highly divergent and presumably arose through accumulation of a number of mutations of small effect, novel phenotypes could have arisen rapidly through recombination between modules both within and between species.

We have demonstrated several aspects of the genetic architecture of wing pattern that have contributed to evolutionary innovation in the Heliconius radiation. First, distinct genetic elements are associated with different patches of red on the butterfly wing. This supports the hypothesis of regulatory modularity, which should facilitate evolutionary innovation. Second, the origin of the dennis-ray phenotype in H. melpomene involved a combination of evolutionary tinkering of existing patterns and introgression between species. Finally, we show that diversity within three lineages (H. elevatus, H. melpomene, and H. timareta) has been generated by shuffling of these distinct regulatory modules among populations and species. Within all three lineages, some populations possess one or other of these elements, providing further flexibility in pattern evolution. Our data imply that recombination between lineages can generate novel phenotypic combinations and demonstrate how modularity in the cis-regulatory control of key genes can drive the rapid evolution of novel morphologies. Although the evolution of novel regulatory modules may involve many mutational steps [5], these can subsequently be exchanged between lineages and shuffled into new combinations enabling rapid adaptive evolution. Recent studies showing that adaptation can proceed via gene flow of preadapted genetic modules between nearby populations or species suggest that similar mechanisms may be important in other radiations. In sticklebacks, adaptation to freshwater involves movement of alleles through the marine landscape [34]. Mosquitoes, Darwin’s finches, and even humans also show evidence for introgression of alleles between species that facilitate adaptation [35–37]. The extent to which recombination between regulatory alleles can contribute to morphological novelty in these other groups of organisms remains to be seen.