Major patterns in the introgression history of Heliconius butterflies

(a) Sample information. (b) Coordinates of chromosome inversion regions considered in this study. (c) Number of coding and noncoding loci by chromosomal region in four versions of the 'etales-9spp' dataset. (d) Proportions of posterior species tree estimates (maximum a posteriori [MAP] trees) from bpp multispecies coalescent (MSC) analysis without gene flow of all four versions of the 'etales-9spp' dataset in blocks of 200 loci, with BNM clade (H. besckei, H. numata, H. melpomene) collapsed into a single tip, summarized into major regions: auto (all autosomes excluding inversion regions), the Z chromosome (chr 21, excluding the inversion region), and individual inversion regions (2b, 6b, 6c, 13b, and 21b). n = number of blocks. Tree indices correspond to those in Figure 1—figure supplement 2. A full list of MAP trees for each chromosomal region is provided in (e). Full results without BNM lumping are provided in (f and g). (e) Proportions of MAP trees from bpp MSC analysis without gene flow of all four versions of the 'etales-9spp' dataset in blocks of 200 loci, with BNM clade collapsed into a single tip, summarized by chromosomal region. Genome-wide summaries are given in (d). (f) Full results from bpp MSC analysis without gene flow of the 'etales-9spp' dataset, summarized into major regions of the genome. See also legend to (d). Trees are in a decreasing order of the total frequencies across all four versions of the dataset. Tree indices correspond to those in Figure 1—figure supplement 1. (g) Full results from bpp MSC analysis without gene flow of the 'etales-9spp' dataset, summarized by chromosomal region. See legend to (f). (h) Bayes factor B = M_1E/M_1D for comparing the two MSC-I models in Figure 1D (M_1E; scenario 1) and Figure 1E (M_1D; scenario 2), calculated for coding and noncoding loci for each chromosomal region using thermodynamic integration with 32 Gaussian quadrature points. M_1E (scenario 2) is preferred if B>100, M_1D (scenario 1) is preferred if B < 0.01, and the test is not significant (n.s.) otherwise (i.e. 100 < B < 0.01). The last three columns are adjusted values of Bayes factor to account for numerical instability of the thermodynamic integration estimates (see ‘Methods’ and Appendix 1—figure 4). (i) Posterior means and 95% HPD intervals (in parentheses) of parameters from the bpp MSC-I model in Figure 1D (scenario 1) using noncoding and coding loci for each chromosomal region of 'etales-8spp' dataset (see a). Multiple posterior peaks, if present, are reported, for example there are two peaks, denoted chr1 and chr1-alt, from the coding loci in chromosome 1. Plots of these estimates with 95% HPD intervals are provided in Figure 1—figure supplement 3. Plots of species tree with estimated divergence times, introgression times, and introgression probabilities are provided in Figure 1—figure supplement 5. (j) Posterior means and 95% HPD intervals (in parentheses) of parameters from the bpp MSC-I model in Figure 1E (scenario 2) using noncoding and coding loci for each chromosomal region. Multiple posterior peaks, if present, are reported, for example, there are two peaks, denoted chr1 and chr1-alt, from the coding loci in chromosome 1. Plots of these estimates with 95% HPD intervals are provided in Figure 1—figure supplement 3. Plots of the species tree with estimated divergence times, introgression times, and introgression probabilities are in Figure 1—figure supplement 6. (k) Maximum likelihood estimates (MLEs) of parameters under the isolation-with-migration (IM) model obtained from 3s for each of the 28 pairs of species in the 'etales-9spp' dataset (H. erato reference, minDP12), using Eueides tales (Tal) as an outgroup. We used the likelihood ratio test (LRT) statistic to test if the model with gene flow (M2) was preferred to the model without gene flow (M0) using the p-value threshold of 1%. The ‘auto’ dataset used all autosomal loci excluding the inversion regions (see c for the number of loci). Plots of migration rate and divergence time estimates are provided in Appendix 1—figure 1. Plots of all parameter estimates with confidence intervals are provided in Figure 1—figure supplement 7. (l) Number of coding and noncoding loci in each chromosomal region in the ‘hmelv25-res' dataset used in bpp analysis and the ‘hmelv25-all' dataset use in 3s analysis. (m) Proportions of MAP trees, with minimum, median, and maximum posterior probabilities shown in parentheses, from bpp MSC analysis without gene flow of the 'hmelv25-res' dataset in blocks of 200 loci, summarized into four major regions: auto (all autosomes excluding 2b and 15b), 2b and 15b inversion regions, and chromosome 21 (Z chromosome). Trees are in a decreasing order of combined frequency across all chromosomes. Tree indices correspond to those in Figure 2—figure supplement 1. (n) Proportions of MAP trees, with minimum, median, and maximum posterior probabilities shown in parentheses from bpp MSC analysis without gene flow in blocks of 200 loci, summarized by chromosomal region. See legend to (m). (o) MLEs of parameters under the MSC-M model obtained from 3s for each of the 55 pairs of species in the 'hmelv25-all' dataset using H. aoede as an outgroup. We used the LRT statistic to test if the model with gene flow (M2) was preferred to the model without gene flow (M0) using the p-value threshold of 1%. See Figure 2—figure supplement 2A for plots of these estimates. The ‘auto’ setting used all autosomal loci excluding the inversion regions. See l for the number of loci. (p) MLEs of parameters under the MSC-M model obtained from 3s for each of the 55 pairs of species in the 'hmelv25-all' dataset using H. erato as an outgroup. See legend to (o). See Figure 2—figure supplement 2 for plots of these estimates. (q) Posterior means and 95% HPD intervals (in parentheses) of parameters obtained from the bpp MSC-I model of Figure 2B using noncoding loci from each non-15b chromosomal region, and under the MSC-I model m1 of Figure 3—figure supplement 3 for the chromosome 15b inversion region. Plots of the estimates are provided in Figure 2—figure supplement 4. 'n/a' means the parameter does not exist in the model. (r) Posterior means and 95% HPD intervals (in parentheses) of parameters from bpp MSC-I models are provided in Figure 2B and Figure 3—figure supplement 3 (model m1) using coding loci. See legend to (q). (s) Posterior means and 95% HPD intervals (in parentheses) of parameters from bpp MSC-M models for the pardalinus-hecale clade: ((Par, Ele), Hec), with six continuous migration rates (M = Nm) among the three species, two species divergence times (τ) and five population sizes (θ). Plot of the estimates is provided in Figure 2—figure supplement 6A. (t) Posterior means and 95% HPD intervals (in parentheses) of parameters obtained from bpp MSC-M models for the pardalinus-elevatus clade: (Par, Ele), with two continuous migration rates (M = Nm). All populations were assumed to have the same population size (θ). Plot of the estimates is provided in Figure 2—figure supplement 6B. (u) Posterior means and 95% HPD intervals (in parentheses) of parameters from bpp MSC-M models for the cydno-melpomene clade: ((Tim, Cyd), Mel), with six continuous migration rates (M = Nm) among the three species, two species divergence times (τ), and five population sizes (θ). Plot of the estimates is provided in Figure 2—figure supplement 6C. (v) Proportions of MAP trees, with minimum, median, and maximum posterior probabilities shown in parentheses, from bpp MSC analysis without gene flow of the ‘silv-chr15' dataset in blocks of 200 loci or 100 loci, summarized by three regions: the inversion region (15b) and flanking regions (15a and 15c). Trees are in decreasing order of combined frequency across all chromosomal regions and from both options of block size. (w) Posterior means and 95% HPD intervals (in parentheses) of parameters obtained from five bpp MSC-I models in Figure 3—figure supplement 3.