Genomic islands of divergence in hybridizing Heliconius butterflies identified by large-scale targeted sequencing

Abstract

Heliconius butterflies represent a recent radiation of species, in which wing pattern divergence has been implicated in speciation. Several loci that control wing pattern phenotypes have been mapped and two were identified through sequencing. These same gene regions play a role in adaptation across the whole Heliconius radiation. Previous studies of population genetic patterns at these regions have sequenced small amplicons. Here, we use targeted next-generation sequence capture to survey patterns of divergence across these entire regions in divergent geographical races and species of Heliconius. This technique was successful both within and between species for obtaining high coverage of almost all coding regions and sufficient coverage of non-coding regions to perform population genetic analyses. We find major peaks of elevated population differentiation between races across hybrid zones, which indicate regions under strong divergent selection. These 'islands' of divergence appear to be more extensive between closely related species, but there is less clear evidence for such islands between more distantly related species at two further points along the 'speciation continuum'. We also sequence fosmid clones across these regions in different Heliconius melpomene races. We find no major structural rearrangements but many relatively large (greater than 1 kb) insertion/deletion events (including gain/loss of transposable elements) that are variable between races.

Figure 1.

Sampling sites, phylogenetic relationships and colour pattern loci. (a) Sampling sites of butterflies used for creating fosmid libraries and for target enrichment sequencing. (b) The phylogenetic relationships of the Heliconius species referred to in the text. (c) Chromosome maps showing the relative positions of the colour pattern loci that are referred to (modified from Baxter et al. [26] and Joron et al. [27]). The wing colour pattern elements controlled by the H. melpomene loci are also indicated in (b).

Figure 2.

Genetic differentiation (F_ST) between H. m. aglaope and H. m. amaryllis across the colour pattern regions (HmYb region (a), HmB/D region (b)) and three unlinked BACs (c). F_ST is uncorrected for sample size and calculated as a 10 kb moving average at 100 bp increments. The threshold in the colour pattern regions indicates the upper 95% CI from 10 000 bootstrap resampling replicates of 1000 bp (the minimum number of sites with data in each 10 kb window) of the unlinked BACs. Peaks of F_ST > 0.3 are shaded in grey. Coding regions (black) and introns (grey stripes) are shown at the bottom of the colour pattern regions; annotations of the unlinked BACs were performed using RNAseq data and automated gene prediction (pipeline to be published in the forthcoming genome paper).

Figure 3.

Genetic differentiation (F_ST) across the colour pattern regions (a,b) and three unlinked BACs (c) at three levels of divergence: within-species—H. m. amaryllis to H. m. aglaope (grey); between closely related species—H. timareta to H. m. aglaope (orange); and between more distantly related species—H. numata to H. m. aglaope (black). The thresholds are the upper 95% CI from 10 000 bootstrap resampling replicates of 1000 bp of the unlinked BACs. Regions showing peaks of F_ST between H. m. agaope and H. m. amaryllis are highlighted in pink. Coding regions (dark green) and introns (pale green stripes) are shown.