Hybridization reveals the evolving genomic architecture of speciation

Abstract

The rate at which genomes diverge during speciation is unknown, as are the physical dynamics of the process. Here, we compare full genome sequences of 32 butterflies, representing five species from a hybridizing Heliconius butterfly community, to examine genome-wide patterns of introgression and infer how divergence evolves during the speciation process. Our analyses reveal that initial divergence is restricted to a small fraction of the genome, largely clustered around known wing-patterning genes. Over time, divergence evolves rapidly, due primarily to the origin of new divergent regions. Furthermore, divergent genomic regions display signatures of both selection and adaptive introgression, demonstrating the link between microevolutionary processes acting within species and the origin of species across macroevolutionary timescales. Our results provide a uniquely comprehensive portrait of the evolving species boundary due to the role that hybridization plays in reducing the background accumulation of divergence at neutral sites.

Figure 1. Five Hybridizing Species of Heliconius in Costa Rica Demonstrate Varying Levels of Genome-wide Differentiation and Gene Flow

(A) Phylogeny of H. cydno, H. pachinus, and H. melpomene, along with their outgroup species, H. hecale and H. ismenius, based on genome sequence data. Their distantly related comimics are shown on the right. (B) Collection sites of individual samples, color-coded according to (A). (C) History of divergence and gene flow among focal taxa based on analysis of genome-wide data using IMa2 (N_e, effective population size; 2Nm, population migration rate). (D) Empirical F_ST distributions among H. cydno, H. pachinus, and H. melpomene, with shading indicating F_ST distributions based on coalescent simulations with and without interspecific gene flow.

Figure 2. Additional Evidence for Gene Flow among Sympatric Species in Costa Rica

(A–D) Sympatric H. melpomene and H. cydno show reduced divergence, measured by both F_ST and d_XY, relative to allopatric comparisons, across two different regions of the genome. Error bars (indicating 95% confidence intervals) and p values are based on bootstrap resampling. (E and F) Furthermore, Patterson’s D statistic is highly elevated in these regions, indicative of biased allele sharing in sympatry due to introgression. Error bars (indicating 95% confidence intervals) and p values are based on bootstrap resampling. *p < 0.05, **p < 0.01, ***p < 0.001.

Figure 3. Signatures of Genomic Differentiation, Focusing on the 12 Regions that Are Divergent between H. cydno and H. pachinus

Known wing color patterning loci (K, Ac, Yb, B/D) are listed, as are genes WntA and Optix. F_ST plots and divergent segment markers are color coded by pairwise comparison.

Figure 4. Z Chromosome and Autosome Divergence in Pairwise Comparisons between Species

Pairwise F_ST represented as boxplots with whiskers between (1) cydno-pachnius (left), (2) cydno-melpomene (middle), and (3) pachinus-melpomene (right) for autosomes versus the Z chromosome, highlighting elevated divergence on the Z chromosome in comparisons with H. melpomene. Similar distributions, separated out by chromosome, are shown in Figure S7.

Figure 5. Dynamics of Genome-wide Divergence during Speciation

(A) Exponential growth in genome-wide divergence compared to linear substitutions as a function of divergence time. Note that d_XY is expressed as the total number of nucleotide substitutions across the genome, rather than a proportion, so the same y axis applies to both the divergence and substitution lines. (B) Venn diagram of the total base-pair overlap between divergent regions in pairwise comparisons.

Figure 6. Divergent Regions of the Genome Exhibit Signatures of Selection and Adaptive Introgression

Each panel shows the mean values of population genetic statistics inside divergent regions (white bars) versus the genomic background (gray bars). Segregating site density (A), π within species (B), derived allele frequency (C), maximum linkage disequilibrium (D), Tajima’s D (E), fraction of fixed differences between species (F), mean pairwise sequence divergence between species (d_XY) (G), and absolute value of Patterson’s D statistic for the four taxon ordering: H. cydno, H. pachinus, H. melpomene, outgroup (H. hecale and H. ismenius) (H). Error bars (indicating 95% confidence intervals) and p values are based on bootstrap resampling. ***p < 0.0001.