Population genomics of parallel hybrid zones in the mimetic butterflies, H. melpomene and H. erato

Abstract

Hybrid zones can be valuable tools for studying evolution and identifying genomic regions responsible for adaptive divergence and underlying phenotypic variation. Hybrid zones between subspecies of Heliconius butterflies can be very narrow and are maintained by strong selection acting on color pattern. The comimetic species, H. erato and H. melpomene, have parallel hybrid zones in which both species undergo a change from one color pattern form to another. We use restriction-associated DNA sequencing to obtain several thousand genome-wide sequence markers and use these to analyze patterns of population divergence across two pairs of parallel hybrid zones in Peru and Ecuador. We compare two approaches for analysis of this type of data-alignment to a reference genome and de novo assembly-and find that alignment gives the best results for species both closely (H. melpomene) and distantly (H. erato, ∼15% divergent) related to the reference sequence. Our results confirm that the color pattern controlling loci account for the majority of divergent regions across the genome, but we also detect other divergent regions apparently unlinked to color pattern differences. We also use association mapping to identify previously unmapped color pattern loci, in particular the Ro locus. Finally, we identify a new cryptic population of H. timareta in Ecuador, which occurs at relatively low altitude and is mimetic with H. melpomene malleti.

Figure 1.

(A) Distribution in South America of the subspecies included in this study. (B) Maximum likelihood phylogenies with approximate likelihood branch supports. Co-mimics from outside the focal hybrid zones are connected with dotted lines. Focal hybrid zone individuals are shown in color. (Blue) H. m. plesseni and H. e. notabilis; (purple) Ecuador hybrids; (dark red) H. m. malleti and H. e. lativitta; (red) H. m. aglaope and H. e. emma; (orange) Peru hybrids; (yellow) H. m. amaryllis and H. e. favorinus. Additional populations are in black. Country abbreviations: (Ec) Ecuador; (FG) French Guiana; (Co) Colombia; (Pa) Panama.

Figure 2.

Population structure at each of the hybrid zones using the reference aligned data. (A) Sampling locations with altitude in meters, sample size in parentheses, and pie charts of the proportion of individuals of each type sampled from each site. Colors are the same as in Figure 1, except black indicates H. timareta in Ecuador. (B) Structure analysis with k = 2 (H. timareta individuals excluded). Each individual is shown as a horizontal bar with the allelic contribution from population 1 (gray) and population 2 (black). (C) Principal components analysis. (D) Distribution of F_ST values from BayeScan.

Figure 3.

Association mapping (A,D) and outlier analysis (B,E) for H. melpomene (A–C) and H. erato (D–F) in Peru. Each phenotype used for the association mapping is shown in a different color as illustrated in C and F. For clarity, only the top 20 associated SNPs are shown for each phenotype. The outlier analysis results show F_ST values for all SNPs, with significant outliers shown in red. Results from the de novo assembled data are shown as crosses (and in orange for the outlier analysis) and positioned based on the top BLAST hit to the H. melpomene genome; those that were not confidently or uniquely assigned to these positions are shown as stars (e.g., those at the end of chromosome 10 in D). (Unmapped) Scaffolds of the H. melpomene reference genome that were not assigned to chromosomes in v1.1 of the genome assembly.

Figure 4.

Association mapping (A,D) and outlier analysis (B,E) for H. melpomene (A–C) and H. erato (D–F) in Ecuador. See Figure 3 legend for further information.

Figure 5.

Venn diagrams of SNPs detected in the de novo assembled (blue and green) and reference aligned (yellow and red) data by BayeScan outlier detection (red and blue) and association mapping (yellow and green), for each of the four populations.