Genomic architecture and introgression shape a butterfly radiation

Abstract

We used 20 de novo genome assemblies to probe the speciation history and architecture of gene flow in rapidly radiating Heliconius butterflies. Our tests to distinguish incomplete lineage sorting from introgression indicate that gene flow has obscured several ancient phylogenetic relationships in this group over large swathes of the genome. Introgressed loci are underrepresented in low-recombination and gene-rich regions, consistent with the purging of foreign alleles more tightly linked to incompatibility loci. Here, we identify a hitherto unknown inversion that traps a color pattern switch locus. We infer that this inversion was transferred between lineages by introgression and is convergent with a similar rearrangement in another part of the genus. These multiple de novo genome sequences enable improved understanding of the importance of introgression and selective processes in adaptive radiation.

Fig. 1:. Phylogeny and phylogenetic networks of Heliconius show lack of support for bifurcating tree.

A. All nodes resolved in a majority of species trees are shown in this cladogram (heavy black lines), while the poorly resolved silvaniform clade is collapsed as a polytomy (Fig. S20). The 500 colored trees were sampled from 10 kb non-overlapping windows and constructed with maximum likelihood. B, C. High-confidence tree structure (black) and introgression events (red) are shown as solid lines. Dashed red lines indicate weakly supported introgression events. Grey branch ends are cosmetic. The melpomene-silvaniform clade is shown in B, the erato-sara clade in C. Euclidean lengths of solid black lines are proportional to genetic distance along the branches. Scale bars in units of substitutions per site. Breaks at the base in B indicate that the branch leading to H. doris has been shortened for display.

Fig. 2:. Local evolutionary history in the erato-sara clade is heterogeneous across the genome.

A. Each bar represents a chromosome, in terms of the H. erato reference (14). Colored bands represent tree topologies of each 50 kb window; colors correspond to the topologies in B, with black regions showing missing data. B. The eight most common trees are shown. The value in the top left corner is the percentage of all 50 kb windows that recover that topology. C. Each histogram corresponds to the topology of the same color in B, and shows the distribution of the number of consecutive 50 kb windows with that topology. Arrows indicate long blocks in inversions.

Fig. 3:. Chromosomal architecture is strongly correlated with local topology.

Tree 1 is shown in red, and Tree 2 is shown in blue, as in Fig. 2. A. Tree 1 shows a negative relationship with chromosome size, while Tree 2 shows a positive relationship. Lines are linear regressions with chromosome 21 excluded. Numbers along top indicate chromosome number. B. Each chromosome was divided into 10 equally sized bins, and the occupancy of each topology in each bin was calculated as the number of windows that recovered the topology in the bin divided by the number of windows that recovered the topology in the chromosome. C. Windows are binned by recombination rate, and boxes show the fraction of each tree in each bin for each chromosome separately. Numbers above boxes are the number of windows in each bin. D. Boxes show the relationship of tree topology with coding density. Asterisk denotes significance at 5% level (paired t-test, p<0.025). In all boxplots, central line is median, box edges are first and third quartile, and whiskers extend to the largest value no further than 1.5*(inter-quartile range).

Fig. 4:. Parallel evolution of a major inversion at the cortex supergene locus.

A. Map of 1.7 Mb region on chromosome 15. Coordinates are in terms of Hmel 2.5, and ticks are in Mb. Tree topology colors correspond to those in Fig. 2. Genes are shown as black rectangles; cortex is highlighted in yellow. Each line shows the mapping of a single contig. Aligned sections of each contig are shown as thick bars, while unaligned sections are shown as dotted lines. Arrows indicate the strand of the alignment. The H. erato group breakpoints are shown with red vertical lines, while the H. numata breakpoints are shown with green vertical lines. B. Evolutionary hypotheses consistent with the topology observed in this inversion in the context of the previously estimated phylogenetic network. The three species used in the triplet gene tree method – H. erato, H. telesiphe, and H. sara – are shown as black lines, while lineages not included are shown as grey lines. C. Histogram of internal branch lengths (T₂) in windows with the topology H. erato, (H. telesiphe, H. sara). The inferred ILS distribution is shown as a dashed line, and the inferred introgression distribution is shown as a dotted line. The average internal branch length in the inversion is shown as a green vertical line. D. Histogram of normalized D_XY (T₃) between H. telesiphe and H. sara. Mean normalized D_XY in the inversion is shown as a green vertical line