Rampant Genome-Wide Admixture across the Heliconius Radiation

Abstract How frequent is gene flow between species? The pattern of evolution is typically portrayed as a phylogenetic tree, yet gene flow between good species may be an important mechanism in diversification, spreading adaptive traits and leading to a complex pattern of phylogenetic incongruence. This process has thus far been studied mainly among a few closely related species, or in geographically restricted areas such as islands, but not on the scale of a continental radiation. Using a genomic representation of 40 out of 47 species in the genus, we demonstrate that admixture has played a role throughout the evolution of the charismatic Neotropical butterflies Heliconius. Modeling of phylogenetic networks based on the exome uncovers up to 13 instances of interspecific gene flow. Admixture is detected among the relatives of Heliconius erato, as well as between the ancient lineages leading to modern clades. Interspecific gene flow played a role throughout the evolution of the genus, although the process has been most frequent in the clade of Heliconius melpomene and relatives. We identify Heliconius hecalesia and relatives as putative hybrids, including new evidence for introgression at the loci controlling the mimetic wing patterns. Models accounting for interspecific gene flow yield a more complete picture of the radiation as a network, which will improve our ability to study trait evolution in a realistic comparative framework.


Introduction
Interspecific hybridization and the resulting gene flow across porous species barriers are increasingly recognized as major processes in evolution, detectable across the tree of life (Hilario and Gogarten 1993;Feliner et al. 2017). Interspecific gene flow has since been demonstrated in several animal taxa, both at deep (Chen et al. 2017) and shallow temporal scales (e.g., Fontaine et al. 2015). Although hybridization might often result in the production of deleterious combinations of alleles at different loci, introgression can also enable adaptation by providing novel variation that may be favored by natural selection, as demonstrated in the iconic adaptive radiation of Darwin's finches (Lamichhaney et al. 2015) and in our own lineage (Huerta-S anchez et al. 2014).
The task ahead is to systematically evaluate the prevalence and importance of interspecific gene flow in fueling speciation in adaptive radiations (Schumer et al. 2014;Feliner et al. 2017). Unfortunately, modeling gene flow requires extensive data and poses greater challenges to computational methods than other processes leading to incongruent signals, such as incomplete lineage sorting or gene duplication (Degnan 2018;Wen et al. 2018). For instance, two studies of the swordtail fishes applied different sequencing and analytical strategies, ultimately reaching dissimilar conclusions on the prevalence of hybridization in Xiphophorus (Cui et al. 2013;Kang et al. 2013). The challenge of characterizing introgression in adaptive radiations remains open and requires both taxonomic completeness and sophisticated methodological approaches.
The Neotropical Heliconius butterflies present an excellent opportunity to study the incidence and importance of gene flow in a recent adaptive radiation, due to the natural propensity of Heliconius and the sister genus Eueides to produce hybrids in the wild (Dasmahapatra et al. 2007;Mallet et al. 2007). The loci responsible for their aposematic wing patterns are especially likely to be shared between species, providing a source of genetic variation in a strongly selected trait and thus likely facilitating speciation (Heliconius Genome Consortium 2012;Pardo-Diaz et al. 2012;Enciso-Romero et al. 2017;Jay et al. 2018;Edelman et al. 2019;Massardo et al. 2020) (summarized in supplementary table 1, Supplementary Material online). In the so called melpomene/cydno/silvaniform clade (MCS), species in sympatry can share variation at up to 40% of the genome (Kronforst et al. 2013;Martin et al. 2013).
It remains unknown whether hybridization and introgression documented in the relatively young MCS clade (4.5-3.5 Ma) (Kozak et al. 2015) are a universal characteristic of this genus. Recent studies of the genus Heliconius based on de novo assembly of genomes (Edelman et al. 2019Seixas et al., 2021 and transcriptomes (Zhang et al. 2019) have identified instances of admixture in other groups, but looked at single individuals in less than half of the 47 recognized species. Furthermore, a study of Heliconius hermathena revealed that a phenotype suggestive of introgression is in fact determined by an ancestral allele expressed in multiple species (Massardo et al. 2020). Full understanding of the frequency of interspecific admixture events requires comprehensive analysis across the diverse species and phenotypic races in the genus.
Here, we generate a comprehensive whole-genome resequencing data set of 145 individuals of 40 among the 47 recognized Heliconius species, and six out of 12 Eueides, encompassing nearly an entire radiation at a continental scale ( fig. 1 and supplementary file 1, Supplementary Material online). With this expanded data set, we investigate the prevalence of hybridization, attempt to quantify its extent across the radiation, and compare the processes producing discordance. We demonstrate varied amounts of phylogenetic incongruence (i.e., conflict between gene trees [Degnan 2018] related to heterogenous levels of gene flow among species). We show that a misleadingly well-supported and resolved tree can be recovered despite incongruence, and that previously unknown, complex hybridization events can thus be missed. Although instances of hybridization across the genome, particularly at adaptive loci, are found across the radiation, we demonstrate that they are far more frequent among the relatives of Heliconius melpomene.

Genomic Trees for the Heliconiines
We first constructed a bifurcating autosomal phylogeny. Mapping genome-wide short read data from 145 individuals in 48 species (supplementary file 1 and table 3, Supplementary Material online) to the H. melpomene reference allowed us to recover 6,848 autosomal and 416 sex-linked high quality, orthologous CDS alignments, respectively (of which 6,725 and 406 include at least one outgroup sequence). The mean length of a quality-trimmed autosomal alignment was 1,387 base pairs (supplementary table 4, Supplementary Material online), and the average parametric aLRT support for the estimated maximum likelihood gene trees was 0.675 6 0.060 with all samples included. Filtering for autosomal exome sites with biallelic, nonsingleton SNPs without missing data produced a 122,913 bp supermatrix. This underpins a maximum likelihood tree resolved with full bootstrap support, except for uncertain placement of the Heliconius clysonymus/hortense/ telesiphe clade (bootstrap support 62/100), which is also placed differently in the coalescent analyses ( fig. 1

Genome-Wide Incongruence and Discordance
Although the species tree topologies recovered by various approaches from the autosomal markers are very similar, multiple indices show incongruence among individual gene trees. The Robinson-Foulds pairwise distance is high: 0.745 out of 1.0 for the autosomal phylogenies and 0.699 for the sex-linked loci, indicating that any two gene trees are very likely to contain multiple differing nodes. Among the 56 nodes separating species and major subspecies, less than a half (26) are resolved in an autosomal majority rule consensus tree (supplementary fig. 4, Supplementary Material online). The relative tree certainty (Salichos et al. 2014) on a 0-1 scale is a low 0.322 when using all gene trees, and increases to 0.397 for the 1,000 best gene trees. Many branches with high support in the coalescent trees score low on the internode certainty (IC) measures (IC/ICA; supplementary fig. 4, Supplementary Material online). Brower and FIG. 1.-Bifurcating model of Heliconiini phylogeny obscures the underlying incongruence. The topology was inferred by the multispecies coalescent method ASTRAL-III from 6,725 orthologous autosomal genes. All nodes, except for the CHT clade are supported in 100% bootstrap replicates. In addition, three types of support values are presented: ASTRAL branch support for a quadripartition (on a 0-1 scale); percentage of quartets in individual gene trees containing a specific node; and the Internode Certainty measure of gene tree incongruence (0-1 scale; X indicates the node is not found in the consensus). Colored circles indicate conflict between the topology of the ASTRAL autosomal tree, and the estimates from: mitochondrial data under ML (blue); sex-linked genes (Z chromosome) in ASTRAL (green); concatenated SNPs under ML (red); the alternative coalescent method MP-EST (violet). Branch lengths are in coalescent units, and arbitrarily set to 1.0 for the terminal branches. Branch colors correspond to previously defined clades (Kozak et al. 2015): red-Heliconius melpomene/cydno; orange-silvaniforms (MCS); violet-Heliconius wallacei; green-Heliconius doris; blue-Heliconius sara; crimson-Heliconius clysonymus; scarlet-Heliconius erato; brown-Eueides.
Garz on-Orduña (2018) suggested that incongruence reported previously (Kozak et al. 2015) was an artifact of missing data. Here, we present nucleotide matrices that are nearly complete (>96%). Although modern statistical phylogenetic methods are typically robust to far higher levels of missing data (Wiens and Morrill 2011;Roure et al. 2013), we still find substantial discordance and incongruence, which shows they are not merely artifacts.
Conflicts are further highlighted by the varied quartet support, whereby many of the nodes reported as certain by ASTRAL are only found in a fraction of the gene trees (second set of support values in fig. 1). This discrepancy is especially exacerbated in the low quartet support for the position of species in the MCS clade, as well as at the placement of the small clades of Heliconius aoede, Heliconius wallacei, Heliconius doris, and H. clysonymus. The statistics are not strongly affected by the exact choice of markers. When the 6,367 single exons are used instead of entire genes, the total normalized quartet score for the entire ASTRAL tree decreases only slightly from 0.847 to 0.806, while the quartet support of individual nodes changes by no more than 10 percentage points (supplementary fig. 3, Supplementary Material online).
We found that the Z chromosome gene trees are more

Hybridization between Species Has Been Common throughout the Radiation of Heliconius
For the first time, we test for gene flow across the entire radiation, including nearly all species (N ¼ 27) in clades other than the previously studied MCS (N ¼ 13). Across the radiation, an inference of species networks reveals a pattern of gene sharing within all major clades of the radiation (figs.

Mosaic Genomes in the Heliconius erato Clade
Proportionally fewer admixture events are identified in the network of the SEC clade (three interspecific admixture edges among 21 lineages; fig. 2B) than among other Heliconius (10 edges between 23 lineages). Nonetheless, support for admixture in this lineage is strong as well. In particular, Heliconius hecalesia shares a large portion of its genome with either the ancestor of the (H. telesiphe,(hortense, clysonymus)) clade (c ¼ 0.38, fig. 2B), or just with H. clysonymus (c ¼ 0.16). Furthermore, there is some evidence for an exchange between the CHT and Heliconius sara clades (c ¼ 0.06). The TreeMix AG also uncovers the CHT-H. hecalesia admixture, but places both in different positions in the tree, such that H. hecalesia appears more closely related to the H. sara clade, and admixes with both CHT and H. erato  The admixture during in the evolution of H. hecalesia and the CHT clade is evident in the pattern of variation among rooted triplets of taxa, examined using the D statistic (Durand et al. 2011) (table 1). The results are highly positive and statistically significant for all tests where H. hecalesia is the recipient of admixture from either the H. clysonymus or H. sara clades (table 1). However, consistent with the phylogenetic patterns, there is evidence for stronger gene flow between H. hecalesia and H. clysonymus (D ¼ 0.35; P < 0.0001) or H. hortense (D ¼ 0.38; P < 0.0001) than the very differently patterned H. sara (D ¼ 0.17; P < 0.0001).
Explicit coalescent modeling also favors models where admixture occurs across species boundaries (fig. 4B-E). In general, these models are consistent with the network and the AG. Specifically, the CHT clade is at the nexus of admixture events, exchanging alleles with H. erato and H. hecalesia

Adaptive Introgression at the Wing Pattern Loci
Our nearly exhaustive sampling of Heliconius species provides a uniform framework in which to gauge the amount of introgression across the functionally important loci that modulate adaptive phenotypic variation. Topologies around wing pattern loci differ from the species tree (P < 0.001, SH test) and in many cases, primarily at the optix and cortex loci, show multiple departures (supplementary table 5 The clustering of H. hecalesia and H. clysonymus/H. hortense is the only case where there is strong evidence for adaptive introgression among the 19 species of the SEC clade. Sequences from the three species indeed cluster exactly at the 360,000:380,000 interval of scaffold HE670865 (aLRT > 0.95; P < 0.001, SH test), which aligns to the specific region of the optix locus controlling red patterns in H. erato (Supple et al. 2013) (fig. 5). This indicates that the alleles governing the red pattern in the three species are more similar than expected from the autosomal phylogenies. H. hecalesia/ clysonymus/hortense cluster also at the wntA interval (HE667780: 450,000:490,000) (Martin et al. 2012), to the exclusion of the phenotypically more different H. telesiphe.
Considering the heterogeneity observed at the rest of the genome, the discordance in the regions associated with wing patterning may be a product of ILS and not hybridization. We tested this possibility by comparing Bayesian species tree and species network models. At an interval within the optix locus (i.e., red patterns; table 2), we find strong support for the network model over the simple tree model (Bayes Factor ¼ 242), and 99.8% of the posterior estimates are networks with at least one admixture edge (m ¼ 3.92, r ¼ 1.47). However, the posterior includes 480 topologies, the most common found only in 11.64% of the posterior. This topology (supplementary fig. 12, Supplementary Material online) implies six incidences of gene flow throughout the SEC clade, and places H. hecalesia in a soft polytomy in the CHT clade. In addition, the inferred age of the CHT clade (1.6-1.4 Ma; supplementary fig. 12, Supplementary Material online) is much lower than expected from a relaxed molecular clock estimate (4.5-2.7 Ma) (Kozak et al. 2015). Although this discrepancy could be caused by the use of a strict clock here, all the other split times are consistent between the strict and relaxed clock estimates.
At the Cr interval (cortex: yellow patterns), there is similarly overwhelming support for a network structure over a tree (BF ¼ 411), and 99.4% of the posterior are networks with an average of 3.24 reticulations (r¼1.54). Unexpectedly, the most frequent topology In the MCS clade, in addition to corroborating previous reports of introgression around color pattern regions, we identify several new cases. For instance, at the cortex locus (see table 2), which is responsible for the diverse white and yellow patterns across the genus (Nadeau et al. 2016), H. melpomene and H. timareta alleles cluster with silvaniforms (scaffold HE667780:310,000:330,000; aLRT > 0.95). At the wntA locus (HE668478:450,000:490,000), sequences of H. heurippa cluster with H. cydno, upholding the view that speciation of the former involved a yellow-patterned race of H. cydno (Enciso-Romero et al. 2017), although the rest of the data places H. heurippa unequivocally as sister to H. timareta (figs. 1-3). Most of the variation in the optix region is consistent with the genome-wide lack of resolution in the H. melpomene/cydno/silvaniform clade and confirms known events. The greatest number of discordant branches are among the H. melpomene/cydno clade at 360,000:380,000 (fig. 5), the section controlling both H. melpomene (Wallbank et al. 2016) and H. erato red ray patterns (Supple et al. 2013;Van Belleghem et al. 2017). Intriguingly, alleles from H. hecale clearei cluster with the Heliconius pardalinus/Heliconius elevatus sequences in eight out of 60 windows on the optix scaffold, perhaps related to the complete loss of orange patterning in this uniquely black and white silvaniform ( fig. 5).

In-Depth Sampling Reveals Widespread Admixture
We interrogated an extensive data set of 6,725 autosomal genes sequenced in nearly all species of a continental-scale adaptive radiation to investigate the prevalence of genomewide admixture. We identified up to 13 cases of gene flow between species as a major source of phylogenetic incongruence ( fig. 2), and demonstrated that admixture shaped the evolution of Heliconius throughout their history. Coalescent modeling revealed admixture between deeply diverged lineages, as well as a complex history of gene flow in the SEC clade of the genus. Although Heliconius is recognized as a foremost example of interspecific gene flow, most of the studies (reviewed in supplementary table 1, Supplementary Material online) focused on H. melpomene and relatives, known to hybridize in the wild with notable frequency (Mallet et al. 2007). Recent studies highlighted new cases in other clades within the genus (Edelman et al. 2019;Zhang et al. 2019), but limited taxonomic and geographic representation of Heliconius diversity made it difficult to assess reliably how many species have admixed (Thawornwattana et al. 2021). Here, we include 40/47 species and highlight the importance of admixture in shaping this complex radiation across time. We used the previously investigated clade of H. melpomene, H. cydno and silvaniforms as a test case, where our approach supports other work documenting extensive admixture, including: hybridization during speciation of H. heurippa (Salazar et al. 2008(Salazar et al. , 2010; admixture between H. cydno/timareta and subspecies of H. melpomene (Martin et al. 2013;Nadeau et al. 2013;Enciso-Romero et al. 2017   can detect known events increases our confidence in the detection of additional instances across the radiation. Inclusion of all species in the SEC clade made it possible to pinpoint the extensive admixture between H. hecalesia and 1) the ancestor of the H. clysonymus clade (c ¼ 0.38); 2) H. clysonymus itself (c ¼ 0.16). Adaptive gene flow between the three species is plausible, as H. hecalesia is sympatric with the other two species in parts of its range (Rosser et al. 2012). H. clysonymus Â H. hecalesia and H. hortense Â H. hecalesia hybrids have been found in the wild (Mallet et al. 2007). To a lesser extent, some degree of gene flow is certain between H. clysonymus and H. sara clades, although even with rich data it remains difficult to reconstruct specific events when several recently diverged species are involved, as the exact parameter values in the coalescent models depend on the sampling of taxa ( fig. 4D and E). The problem is especially acute in the reconstruction of introgression histories at the wing pattern loci, where no specific topology is strongly supported, and even top-scoring networks are difficult to interpret given the differences in wing phenotypes of putatively introgressing species ( fig. 5:1-4 and supplementary figs. 12-14, Supplementary Material online). As variation in Heliconius wing patterns appears to be governed by short regulatory elements that differ even between mimics (Concha et al. 2019), detailed investigation will be necessary to identify the specific functional regions within the broader intervals investigated here. Nonetheless, Heliconius are unusual in that introgression of unlinked loci enables rapid evolution of complex patterns, which comprise a patchwork of elements sometimes derived from different sources (Wallbank et al. 2016). Many genomic studies of interspecific gene flow have found introgressions of small genome regions driven by natural selection on beneficial alleles, such as multiple abiotic tolerance factors in Helianthus debilis into Helianthus annuus (Whitney et al. 2010), the hypoxia resistance EPAS1 haplotype (Denisovans ! anatomically modern Tibetans) (Huerta-S anchez et al. 2014), the ALX1 alleles determining diverse beak shapes among Darwin's finches (Geospiza) (Lamichhaney et al. 2015), or the Agouti variant conferring protective coat color (Lepus americanus ! Lepus timidus) (Giska 2019). Only in a few other systems is there evidence for adaptive introgression at multiple loci, including hominins (reviewed by Gokcumen 2020), and the Lonchura finches (Stryjewski and Sorenson 2017). Similar to Lonchura, the evolution of the key adaptive trait in Heliconius (patterning) has involved introgressions at multiple loci and between different combinations of species.

Challenges of Inferring Interspecific Gene Flow
While admixture is rampant, it remains difficult to describe it with precision. Even though all approaches suggest that gene flow occurred, the exact sources and direction are not estimated consistently between methods. The PhyloNet maximum pseudolikelihood (MPL) networks contain 13 reticulation edges ( fig. 2), five of which are also recovered  (Sheppard et al. 1985). However, more recent research has demonstrated that loci that have been defined from intraspecific crosses in different species map to homologous regions of the genome (e.g., see Joron et al. [2006]). Moreover, candidate protein-coding genes have been identified and, in some cases, the intervals containing functional variation have been localized (see Key References). Scaffold numbers refer to the Hmel v1 assembly. HW, hindwing; FW, forewing. a Brown patterns in Heliconius cydno and Heliconius pachinus (Chamberlain et al. 2011). b The Pushmipullyu supergene controlling most of the wing patterning in Heliconius numata (Joron et al. 2011). by the TreeMix AG ( fig. 3): H. hecalesia-CHT clade; H. hecalesia-H. erato; H. melpomene-H. timareta; H. melpomene-H. cydno; H. melpomene clade-silvaniforms. The TreeMix AG does not detect the exchanges in and between the small H. egeria and H. hecuba clades, or some of the events previously documented between species of the silvaniform group and H. melpomene (Jay et al. 2018;Zhang et al. 2016). The discrepancies are expected between two widely different techniques, as the TreeMix AG algorithm assumes that the underlying sequence of events was largely tree-like (Pickrell and Pritchard 2012). The AG approach, based on allele frequencies, was designed with assumptions more appropriate at the level of recently diverged taxa and may be affected by issues of multiple substitution. Similarly, the presented D statistics need to be taken with caution. Although the factors affecting the specificity of D have not been formally determined, it is likely to be affected in clades more distant from the H. melpomene reference genome, as worse read mapping results in lower overall number of sites for comparison (supplementary table 3, Supplementary Material online) and thus possibly an unfavorable signal-to-noise ratio. In comparison, the sensitivity of the D statistic decreases both when the population size is large relative to the divergence time, as is the case for many widespread Heliconius species, and when gene flow was ancient (Zheng and Janke 2018). More surprising is the disparity between two approaches computing over gene trees, MPL networks and PHRAPL. In case of the latter, the ability to evaluate a large number of models with an extensive data set of thousands of gene trees comes at the cost of less accurate parameter estimates, e.g. when compared with Approximate Bayesian Computation approaches (Jackson et al. 2017). Furthermore, the computational burden is reduced by limiting the questions to rooted triplets of taxa and subsampling intraspecific allelic diversity, thus losing many of the benefits of comprehensive sampling.
Other recent studies of introgression among Heliconius encountered similar difficulties. For instance, an analysis of whole genomes of 20 species (Edelman et al. 2019) identified the same key patterns (e.g., uncertain placement of H. hecalesia; gene flow H. hecalesia-CHT; H. melpomene-silvaniforms), but with overall low confidence and without the ability to ascertain if the proposed events involved unobserved lineages. Recent application of full likelihood coalescent modeling produced more robust results (Thawornwattana et al. 2021), but included only six out of 20 species in the SEC clade, making it impossible to infer the exact sources of admixture. The representation of intraspecific variation is also important: in our study we sample only the nonmimetic H. hermathena hermathena and thus cannot replicate the results of Massardo et al. (2020), who discovered introgression at cortex between H. erato and its mimic H. hermathena vereatta. Despite the difficulties in matching sufficient data with robust analytical tools, all approaches used in our and other studies point to H. hecalesia as a product of hybridization.

Neither Concatenation nor Coalescent Trees Adequately Represent Species History
There has been a marked shift over recent years away from phylogenetic methods that involve concatenation of data, and toward approaches that involve coalescent modeling. Methods for inferring a species tree by modeling the incomplete sorting of loci represent an improvement on the assumption that there is a common evolutionary history across all genomic regions (Heled and Drummond 2010;Liu et al. 2010;Mirarab et al. 2014), although the variation in results demonstrates the need for better analytical tools, as well as more complete data. Across the tree of life, from birds (Reddy et al. 2017) and mammals (Chen et al. 2017) to land plants (Zhong et al. 2013) and fungi (Shen et al. 2016), treatment of individual gene trees under MSC methods has yielded substantially different results to simple concatenation. In contrast, our Heliconius trees are consistent with previous work Kozak et al. 2015;Zhang et al. 2016;Edelman et al. 2019), but clarify some uncertainties, including the placement of the H. hecuba and H. egeria groups, relationships in the H. sapho clade, and the position of H. besckei. Nonetheless, similar to other large phylogenetic studies (Brawand et al. 2014;Fontaine et al. 2015), none of the individual gene trees showed exactly the same topology as the autosomal MP-EST species tree, suggesting that the wellsupported bifurcating trees do not fully represent the underlying signal in the genomes within this clade. Network modeling clearly demonstrates that introgression has been important throughout the evolution of the genus, and yet this process could easily be overlooked with many of the modern phylogenomic methods.
The comprehensive analysis of the large butterfly genus shows the important role of adaptive introgression at multiple loci in shaping radiations. The main appeal of studying adaptive radiations is their power for analyzing trait evolution in a comparative framework, and a growing number of studies are looking at several Heliconius characters through this lens (e.g., Briscoe et al. 2013;Sculfort et al. 2020). It is increasingly clear that many key adaptive traits are determined by introgressed sequences, and thus a comparative approach reliant on a single bifurcating species tree would give highly incomplete results (Hahn and Nakhleh 2016;Bastide et al. 2018). To expose the hidden uncertainties of phylogeny and capture the potential of adaptations to be shared between species, future work must utilize approaches that reflect the nonbifurcating reality of evolving genomes.

Sampling
Heliconius can be divided into two deep lineages ( fig. 1 and  supplementary file 1, Supplementary Material online). The first consists of H. erato, H. sara, H. clysonymus and relatives, and throughout the text we refer to this group as SEC. This lineage can be subdivided into three smaller clades of species closely related to H. sara, H. erato, and H. clysonymus. The clade of H. clysonymus, H. hortense, and H. telesiphe turned to be of special interest and we refer to it as CHT. The second major lineage within the genus are species related to H. melpomene, divided into five groups: the H. melpomene/H. cydno group; H. numata and relatives, often called "silvaniforms"; and the clades of H. doris, H. wallacei, and H. aoede. The first two clades are often grouped together and referred to as the H. melpomene/cydno/silvaniform clade (MCS).
We sampled 40 out of 47 Heliconius, as well six of the 12 species in the sister genus Eueides, and the monotypic Dryadula and Agraulis as outgroups (supplementary file 1, Supplementary Material online). Genomes of 11 species were re-sequenced for the first time: Heliconius atthis, Heliconius antiochus, Heliconius egeria, Heliconius leucadia, Heliconius peruvianus, Eueides aliphera, Eueides lampeto, Eueides lineata, Eueides isabella, Eueides vibilia, Agraulis vanillae. Material of sufficient quality could not be obtained for the remaining seven species of Heliconius and six of Eueides. Data for the other 37 species included in the study were published previously (Heliconius Genome Consortium 2012; Briscoe et al. 2013;Kronforst et al. 2013;Martin et al. 2013;Supple et al. 2013;Nadeau et al. 2016;Wallbank et al. 2016;Enciso-Romero et al. 2017;Jay et al. 2018). To enhance coalescent modeling by sampling genetic diversity (Edwards et al. 2016), we included individuals from distant populations and diverse wing pattern races when possible. Our full data set totaled 145 individuals and included multiple individuals of most species.

DNA Sequencing
All sequencing data used in this study, novel and previously published, were generated with the Illumina technology with 100 bp paired-end reads, insert sizes of 250-500 bp and read coverage from 12Â to 110Â. In case of the new samples, DNA was extracted with the DNeasy Blood and Tissue kit (Qiagen) from 30 to 50 lg of thorax tissue homogenized in buffer ATL using the TissueLyser (Qiagen); purified by digesting with RnaseA (Qiagen); and quantified on a Qubit v.1 spectrophotometer (LifeTechnologies). Whole genome libraries with an average insert size of 500 bp were sequenced on a HiSeq 2500 to a mean coverage depth of 50.7Â (range: 33.1-67.8Â).

Exome Alignments and Gene Trees
Protein-coding genes can be effectively treated as discrete markers for multilocus phylogenetics (Edwards et al. 2016). Exonic markers were chosen over noncoding loci because 1) reads from distantly related species map better at the CDS; 2) orthologous sequences can be identified with greater confidence. We minimized paralogy by narrowing the gene set to 1:1:1 orthologs between H. melpomene, Danaus plexippus, and Bombyx mori identified by OrthoMCL (Li et al. 2003;Heliconius Genome Consortium 2012). Alignments of entire protein-coding, single-copy genes were extracted with an inhouse script (gene_fasta_from_reseq.py [Martin 2017]). For this analysis, we ignored all genes within scaffolds linked to the color pattern loci (see table 2), as the exact boundaries of these loci have not been established for most species (Van Belleghem et al. 2017) and linkage to loci involved in adaptive color pattern differences might mislead phylogenetic inference.
We trimmed the alignments with TrimAl v1.2 (Capella-Guti errez et al. 2009), removing any sequences that contained >50% missing data (see the command line in Supplementary Methods). Furthermore, high entropy sections of each alignment were excluded by Block Mapping and Gathering with Entropy (BMGE) (Criscuolo and Gribaldo 2010) with a moderately relaxed PAM100 similarity matrix. Individual ML gene trees were estimated in FastTree v2.1 (Price et al. 2010) with parametric aLRT nodal support (Anisimova and Gascuel 2006). Species tree and network analyses listed below were conducted using rooted gene trees inferred from the 6,725 autosomal and 406 sex-linked CDS genes.

Incongruence in the Data
To assess how much the topologies of gene trees differ from one another, we calculated the Robinson-Foulds distance (RF) (Robinson and Foulds 1981) for all pairs of trees, normalized by dividing the observed distance by the maximum possible RF between the two trees. This statistic was calculated using PAUP* v4 (Swofford 2002) across the entire data set of 145 samples and 6725 genes. As some of the differences are expected to arise from the lack of intraspecific resolution, the calculation was repeated on a thinned data set of 57 high coverage individuals representing all species, with additional individuals included in species with strong geographic structure (supplementary file 1, Supplementary Material online).
We identified highly incongruent nodes by computing 50% majority rule trees and using four information-theory measures (Salichos and Rokas 2013;Salichos et al. 2014) on the reduced data set of 57 samples at 6,725 genes. These analyzes were performed in RAxML v8 (Stamatakis 2014). The information-theory measures included the: 1) IC, which compares the frequency of a bipartition to the frequency of the most common alternative; 2) the IC All (ICA), which considers all alternatives with support ! 5%; 3) the corresponding tree certainty (TC), calculated and normalized over all bipartitions: and 4) the tree certainty All (TCA) (Salichos et al. 2014). All four measures are expressed on a scale from -1 (when the bipartition is not found in any of the gene trees) to 1 (when the bipartition is found in all). Because various patterns of missing data can impact the IC/TC values (Salichos et al. 2014), and many of our gene trees are incompletely resolved, we examined the impact of input data quality on these scores by repeating the procedure with the 1,000 most resolved phylogenies.

Species Trees
Naïve "total evidence" phylogenies were estimated from concatenated exonic SNPs in RAxML v8 with 100 bootstrap replicates, GTRþC model and ascertainment bias correction (Stamatakis 2014). The history of the matriline was approximated from the whole-mitochondrial alignment with partitions determined by PartitionFinder v1.1 (Lanfear et al. 2012). In addition, we used two MSC approaches to estimate the species tree under the assumption of Incomplete Lineage Sorting. MP-EST v1.4 maximizes a pseudo-likelihood function over the distribution of taxon triples extracted from gene tree topologies, and provides a measure of incongruence based on the proportion of triples shared by the gene trees, similar to the RF distance (Liu et al. 2010). An MSC phylogeny was inferred in MP-EST from the 6,725 gene trees, evaluating support by re-estimating the tree 100 times with random samples of 500 gene trees. Since MP-EST may be misled by errors in gene tree reconstruction (Mirarab and Warnow 2015), we compared the results with ASTRAL-III, a fast quartet method that accounts for polytomies and low support values in the input . ASTRAL quantifies discordance by computing how many of the gene trees contain the quartets making up the species tree (Sayyari and Mirarab 2016).
Recombination between exons scattered across a long genomic interval may lead to conflicting signals, which could be obscured by performing phylogenetic inference at the level of a complete coding gene, a practice criticized as "concatalescence" (Gatesy and Springer 2013). To account for this possibility, we repeated the ASTRAL-III species tree inference using individual autosomal exons longer than 500 base pairs. By restricting analysis to these 6,367 longer exons, we ensured sufficient information content, while eliminating the possibility that recombination between exons of one gene was responsible for conflicting signals.

Admixture Networks
Following the identification of problematic nodes based on gene tree statistics, information criteria and conflicting nodes in species trees, we applied two distinct network approaches (Hahn and Nakhleh 2016). First, we determined the admixture graph (AG) in TreeMix v1.13 by identifying the pairs of taxa sharing more than the expected proportion of allelic variation (Pickrell and Pritchard 2012). Individuals were again assigned to taxa, distinguishing major clades within well-represented species as separate lineages. Relations between taxa were inferred from allele frequency data computed in PLINK v1.9 (Chang et al. 2015), based on the matrix of autosomal SNPs. As TreeMix assumes individual SNPs to be represented across samples and independent, the original matrix was filtered to remove sites with <95% complete data. We identified and pruned sites that could be linked within species, using the pairwise linkage disequilibrium estimator in PLINK v1.9 with default settings.
Second, we modeled both hybridization and incomplete lineage sorting under the MSC network framework implemented in PhyloNet v3.5. Networks were computed under the MPL criterion from the 6725 rooted autosomal phylogenies, considering only nodes with support >0.8 (-b 0.8) and starting with the MP-EST species tree. To search for ancient admixtures between deeper branches of the tree, we conducted an analysis with one species from each of the seven Heliconius clades (H. melpomene, H. numata, H. doris, H. wallacei, H. erato, H. telesiphe, H. aoede) and Eueides, as a full run with all the samples was not feasible. The optimal network was determined by calculating the Bayesian information criterion from the maximum likelihood and the number of lineages, admixture edges and gene trees in each model (supplementary file 3, Supplementary Material online) (Yu et al. 2014).

Gene Flow in the H. erato Lineage
Admixture has been studied only in a few species of the SEC clade (Edelman et al. 2019;Massardo et al. 2020). Here, we focused on the events involving H. hecalesia, H. clysonymus and cognates. The extent of genome-wide similarity between species clusters was illustrated with PCAs of variation in the matrix of autosomal SNPs, calculated for the SEC clade in the R package adegenet (Jombart and Ahmed 2011). To test for gene flow we calculated the D statistic (Durand et al. 2011), derived from the allelic configurations of two taxa, P 1 and P 2 , their relative P 3 and an outgroup. Under the null hypothesis, no admixture occurred between P 3 and either P 1 or P 2 . Under the alternative, P 3 exchanged alleles with either P 1 or P 2 , resulting in the excess of the corresponding allelic configuration. Significance of the result can be assessed by block jackknifing (Martin et al. 2015). Specific tests were conducted for gene flow between H. hecalesia (P 2 ) and H. clysonymus, H. hortense, H. telesiphe or species from the H. sara clade (P 3 ). The sister species P 1 was either the allopatric H. erato from French Guiana, or the parapatric H. erato from Amazonia, and the outgroup was always H. melpomene.
To address the hypothesis of hybrid origin of H. hecalesia or the H. clysonymus/H. hortense pair explicitly, we modeled various scenarios of divergence with gene flow under the maximum likelihood criterion as implemented in PHRAPL (Jackson et al. 2017). PHRAPL is a maximum likelihood method to assess complex scenarios including lineage coalescence, gene flow and population growth. PHRAPL uses simulations to compare all possible scenarios under a range of parameter values sampled from a predefined grid, using gene tree topology as the summary statistic (Jackson et al. 2017). We fitted models of speciation history to four distinct triplets of taxa: 1) H. erato (Western and Andean populations), H. hecalesia and H. clysonymus (including the sister species H. hortense); 2) H. telesiphe, H. clysonymus and H. charithonia (including the sister species H. peruvianus); 3) H. erato, H. hecalesia and H. telesiphe; 4) H. telesiphe, H. charithonia and H. sara. The four sets of taxa were chosen based on the earlier evidence for introgression from either MPL networks and TreeMix (1, 2) or the D statistics (1-4). To diminish the computational burden, we subsampled up to three tips per species from the 6,725 gene trees with 30 replicates per tree (Jackson et al. 2017). For each triplet, we evaluated a total of 48 models, including parameters for coalescence, symmetric migration, and change in population size between lineages, estimated at the default values for a PHRAPL grid search. The optimal model was selected by weighted Akaike's information criterion (wAIC).

Introgression at the Color Pattern Loci
The history of the loci associated with aposematic wing phenotypes is distinct due to strong selection (Mö st et al. 2020), and thus the seven scaffolds in the Hmel v1 assembly containing these loci (table 2) were treated separately. Each scaffold alignment was partitioned into windows of 20 kb, sliding by 10 kb, discarding windows with <1000 polymorphic sites [script sliPhy3.py (Martin 2017)]. The topology for every window was reconstructed and tested for significant differences from the MP-EST species tree using the SH test (Shimodaira and Hasegawa 1999) in RAxML. In order to understand precisely how the Hmel v1 reference corresponds to the specific red control loci of H. erato (Supple et al. 2013), the optix (B/D) scaffold HE670865 was aligned against the H. erato B/D BACs (Papa et al. 2008) with mLAGAN (Brudno et al. 2003).
After finding discordant topologies at the optix, cortex, and WntA loci, we tested if the variation can be explained by incomplete lineage sorting alone, or whether gene flow in the SEC clade needs to be taken into account. We focused on specific intervals around the key genes, where departures from the species tree were found for SEC: HE670865:310,000:460,000, containing the optix protein CDS and several regulatory elements (Supple et al. 2013;Van Belleghem et al. 2017); HE667780:570000:750000, including cortex (Nadeau et al. 2016;Van Belleghem et al. 2017) and the novel inversion identified in some species of the SEC clade (Edelman et al. 2019); HE668478:450,000:500,000, containing WntA (table 2).
We analyzed all three genomic intervals separately under coalescent models in BEAST2. In each case, the interval was divided into windows of 10 kb treated as separate partitions, and the alignments were reduced to the relevant species: H. hecalesia; the CHT clade; H. sara and H. charithonia; H. erato, H. himera and H. hermathena. First, the coalescent tree model, where incongruence is purely due to incomplete lineage sorting, was fitted in StarBEAST2 (Ogilvie et al. 2017). Each partition was assigned the HKYþC substitution model with four rate categories, and a speciation-only tree model was selected for compatibility with the species network analysis. Based on the age of the Eueides-Heliconius split (Chazot et al. 2019), we used a strict molecular clock with an estimated exome-wide rate of 0.003 substitutions per site per MY. Two independent chains of 200 million cycles with 50% burnin were executed for each analysis; possible bias in priors was assessed by executing empty prior runs; and the convergence and effective sample sizes of the numeric parameters were visualized in Tracer (Rambaut et al. 2014). Second, the coalescent network was estimated in BEAST2 (Zhang et al. 2018): a model where incongruence can be a result of either ILS or admixture between species. Additional priors included: species net diversification rate (exponential with mean 1.0, corresponding to doubling every 1 My) and turnover rate (beta distribution: a ¼ 3.0; b ¼ 1.0, corresponding to low probability of admixture). The relative fit of the tree and network models to the same data was estimated by computing Bayes Factors from marginal likelihoods determined by Path Sampling with 20 steps and chains of 10 million cycles