Ancestral gene flow at the base of Heliconius

(A) Posterior probabilities of species trees for major Heliconius clades inferred from BPP analysis of 200-locus blocks (each spanning 1–3 Mb) across the genome under the multispecies coalescent (MSC) model with no gene flow. The y-axis shows the posterior probability of species trees and ranges from 0 to 1. Colors correspond to the nine most common MAP trees, summarized by lumping species in the melpomene-silvaniform clade (Bes, Num, Mel) into a single tip (BNM); see Figure S1 for full trees. Proportions of coding and noncoding blocks with each tree as a MAP tree are shown in parentheses. (B, C) Two scenarios of early divergence of Heliconius: (B) erato-early versus (C) aoede-early. Each tree is a MAP tree from a block having the MAP tree supporting each scenario, with estimates of branch lengths (posterior means). (D, E) Introgression histories estimated under an MSC model with introgression (MSci) corresponding to the two scenarios based on coding loci in chromosome 1 (3,517 coding loci; see Figures S7–S8). (F) Proportions of trees of scenarios 1 or 2 in each chromosome in order of increasing number of loci (used as a proxy for chromosome size). Chromosome 21 is placed at the right end. Number of blocks is shown on top of each bar. Tal: Eueides tales, Mel: H. melpomene, Bes: H. besckei, Num: H. numata, Bur: H. burneyi, Dor: H. doris, Aoe: H. aoede, Era: H. erato, Sar: H. sara.

Major introgression events in the melpomene-silvaniform clade

(A) Blockwise estimates of species trees of the melpomene-silvaniform clade inferred from 200-locus blocks across the genome under the MSC model with no gene flow using BPP (see Table S10 for data, Tables S11–S12 for full results). MAP trees are labelled in decreasing order of frequency across blocks. Proportions of coding and noncoding blocks with each tree as a MAP tree are shown in parentheses. (B) The MSci model with introgression events that can explain the three major groups of trees in (A). Branch lengths are based on posterior means of divergence/introgression times estimated from 5,341 noncoding loci in chromosome 1 (Table S15). Each internal node is given a label, which is used to refer to a population above the node, e.g. the population between nodes r and bn is referred to as branch bn. Each horizontal arrow represents a unidirectional introgression event, e.g. the arrow from tcmeph1 to n3 represents tcmeph→Num introgression at time τtcmeph1 = τn3 with probability φtcmeph1→n3. (C) Continuous migration (IM) model for the pardalinus-hecale clade, allowing bidirectional gene flow among three species. (D) IM model for the cydno-melpomene clade. For C and D, branch lengths are based on estimates from noncoding loci in chromosome 1 (Figure S14; Tables S17–19). Arrow sizes are proportional to estimated migration rate (M = Nm). Bes: H. besckei, Num: H. numata, Par: H. pardalinus, Ele: H. elevatus, Hec: H. hecale, Tim: H. timareta, Cyd: H. cydno, Mel: H. melpomene.

Introgression history of the chromosome 15 inversion region (15b)

(A) Blockwise estimates of species trees of the inversion (15b) and its flanking regions (15a and 15c) on chromosome 15 in the melpomene-silvaniform clade inferred from 200-locus blocks and 100-locus blocks under the MSC model without gene flow using BPP (Table S20). Tree legends are grouped by whether H. numata clusters with H. ismenius (blue), or H. numata with P1 inversion (Num11) clusters with H. pardalinus (red), or other relationships (green). (B–E) Possible scenarios of the origin and introgression route of the P1 inversion. Ism: H. ismenius, Num00: H. numata homozygous uninverted 15b (H. n. laura and H. n. silvana), Num11: H. numata homozygous for inversion P1 (H. n. bicoloratus), Eth: H. ethilla. See Figure 2 legend for codes for other species.

Revised phylogeny and introgression history of Heliconius

Phylogeny of the erato-sara clade was previously estimated using a similar approach [26]. Arrows represent introgression events; dotted lines indicate events that are uncertain. Introgression of chromosomal inversions are region-specific (13b and 15b). Grey shading represents a period of continuous gene flow. Triangle represents the origin of pollen feeding near the base of Heliconius. Stars indicate two alternatives for where the P1 inversion on chromosome 15 originated; green branches indicate lineages with the inversion (polymorphic in Num, fixed in Par). Two rows at the bottom show status of inversions 13b and 15b, with + indicating inversion, – indicating standard orientation, and +/– inversion polymorphism. Him: H. himera, Sia: H. hecalesia, Tel: H. telesiphe, Dem: H. demeter; see Figures 1–3 legends for other species codes.