Selection and gene flow define polygenic barriers between incipient butterfly species

Characterizing the genetic architecture of species boundaries remains a difficult task. Hybridizing species provide a powerful system to identify the factors that shape genomic variation and, ultimately, identify the regions of the genome that maintain species boundaries. Unfortunately, complex histories of isolation, admixture and selection can generate heterogenous genomic landscapes of divergence which make inferences about the regions that are responsible for species boundaries problematic. However, as the signal of admixture and selection on genomic loci varies with recombination rate, their relationship can be used to infer their relative importance during speciation. Here, we explore patterns of genomic divergence, admixture and recombination rate among hybridizing lineages across the Heliconius erato radiation. We focus on the incipient species, H. erato and H. himera, and distinguish the processes that drive genomic divergence across three contact zones where they frequently hybridize. Using demographic modeling and simulations, we infer that periods of isolation and selection have been major causes of genome-wide correlation patterns between recombination rate and divergence between these incipient species. Upon secondary contact, we found surprisingly highly asymmetrical introgression between the species pair, with a paucity of H. erato alleles introgressing into the H. himera genomes. We suggest that this signal may result from a current polygenic species boundary between the hybridizing lineages. These results contribute to a growing appreciation for the importance of polygenic architectures of species boundaries and pervasive genome-wide selection during the early stages of speciation with gene flow.