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.

Keywords: adaptive introgression; admixture; phylogenomics; radiation.

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.

Fig. 2.

Evidence of introgression is found across the entire Heliconius radiation. Networks inferred under maximum pseudolikelihood (MPL) based on 6,725 autosomal ML gene trees distinguish between introgression and incomplete lineage sorting, revealing several admixture events. Numbers on the edges indicate the inheritance probabilities (Wen et al. 2018), which correspond to the proportion of the data supporting the grouping. (A) The analysis of Heliconius melpomene and cognates reveals previously undetected introgressions closer to the root of the tree. Known events in the MCS clade are recapitulated, demonstrating sensitivity of the approach. (B) Fewer admixtures occurred in the evolution of the Heliconius erato/sara clade, but Heliconius hecalesia may be a recent hybrid.

Fig. 3.

The extent of interspecific gene flow varies across the tree. TreeMix inference of splits and mixture from autosomal SNPs. Migration edges (1–8) are inferred on a phylogenetic tree built from allele frequencies under a Gaussian genetic drift approximation. Colors of the edges correspond to the proportion of the genome exchanged.

Fig. 4.

Ambiguous genomic composition of Heliconius hecalesia. (A) Although usually recovered as the sister species of Heliconius hermathena and Heliconius erato in bifurcating phylogenies, H. hecalesia shares variation with the clades of Heliconius sara and Heliconius clysonymus. Principal Component Analysis of variation in the autosomal SNPs within the H. erato/sapho clade. First two PCs account for over half of the variation. (B) Models of divergence history in the H. erato/Heliconius telesiphe group inferred by PHRAPL, where orange edges indicate gene flow. Models including up to five terms for coalescence, gene flow and variation in population sizes were fitted to four sets of three taxa, extracted from the 6725 gene trees. The bar in panel 4E shows time in units of 4N (diploid population size; see supplementary table S6, Supplementary Material online).

Fig. 5.

Discordance at the red patterning locus. Branches in positions varying from the species tree are labeled red. The ML tree was estimated for a 20 kb region on the optix scaffold (HE670865:360,000:380,000), containing the intervals most strongly associated with wing patterns in Heliconius erato and Heliconius melpomene. Intraspecific relations for species not discussed in text are collapsed. Outgroups and parametric support values <0.9 not shown. 1. Heliconius telesiphe sotericus, 2. Heliconius hortense, 3. Heliconius clysonymus hygiana, 4. Heliconius hecalesia formosus; 5. Heliconius hecale felix, 6. Heliconius hecale clearei; 7. Heliconius timareta timareta, 8. Heliconius melpomene malleti, 9. Heliconius elevatus; 10. Heliconius melpomene melpomene, 11. Heliconius heurippa.