Genomic hotspots for adaptation: the population genetics of Müllerian mimicry in the Heliconius melpomene clade

Abstract

Wing patterning in Heliconius butterflies is a longstanding example of both Müllerian mimicry and phenotypic radiation under strong natural selection. The loci controlling such patterns are "hotspots" for adaptive evolution with great allelic diversity across different species in the genus. We characterise nucleotide variation, genotype-by-phenotype associations, linkage disequilibrium, and candidate gene expression at two loci and across multiple hybrid zones in Heliconius melpomene and relatives. Alleles at HmB control the presence or absence of the red forewing band, while alleles at HmYb control the yellow hindwing bar. Across HmYb two regions, separated by approximately 100 kb, show significant genotype-by-phenotype associations that are replicated across independent hybrid zones. In contrast, at HmB a single peak of association indicates the likely position of functional sites at three genes, encoding a kinesin, a G-protein coupled receptor, and an mRNA splicing factor. At both HmYb and HmB there is evidence for enhanced linkage disequilibrium (LD) between associated sites separated by up to 14 kb, suggesting that multiple sites are under selection. However, there was no evidence for reduced variation or deviations from neutrality that might indicate a recent selective sweep, consistent with these alleles being relatively old. Of the three genes showing an association with the HmB locus, the kinesin shows differences in wing disc expression between races that are replicated in the co-mimic, Heliconius erato, providing striking evidence for parallel changes in gene expression between Müllerian co-mimics. Wing patterning loci in Heliconius melpomene therefore show a haplotype structure maintained by selection, but no evidence for a recent selective sweep. The complex genetic pattern contrasts with the simple genetic basis of many adaptive traits studied previously, but may provide a better model for most adaptation in natural populations that has arisen over millions rather than tens of years.

Figure 1. Annotation of the HmB and HmYb gene regions.

Sequenced H. melpomene BAC clones across the 721 kb HmB locus and the 330 kb HmYb locus as defined by recombination mapping. In total 24 predicted genes were identified in the HmB region and 20 in the HmYb region. AEHM-22C5 and AEHM-7G5 are overlapping clones from different fingerprinted contigs, and both contain part of the Mad gene. Also see Table S1 for HmB annotation and ref. 65 for HmYb annotation.

Figure 2. Locations of field collected samples used for population genetic analysis.

H.m. rosina, H.m. melpomene, and H.m. amaryllis display the medial red forewing band controlled by a dominant allele at HmB. In addition, H. pachinus, H.m. rosina, and H.m. amaryllis display the yellow hindwing bar controlled by a recessive allele at HmYb.

Figure 3. Genotype-by-phenotype associations across three phenotypic comparisons.

(A) The HmYb region between H. m. melpomene and H. m. rosina (triangles) and between H. m. aglaope and H. m. amaryllis (circles), (B) The HmYb region between H. pachinus and H. cydno, and (C) the HmB region between H. m. aglaope and H. m. amaryllis. Significance cut-off was calculated using the Bonferroni correction, applied across all 866 SNPs tested.

Figure 4. Linkage disequilibrium (LD) for combined H. m. aglaope and H. m. amaryllis populations.

Pairwise LD between sites across the genomic region is shown as a heatmap, with greyscale indicating estimated r² values. Darker colours indicate a stronger correlation between SNP variants at the two sites, across the individuals sampled. The highest LD estimates are all within gene regions, or between closely linked loci. The largest block of LD is across the region with strong genotype-by-phenotype association in the HmB region, across the Slu7, Kinesin and GPCR genes (A, lower left). Comparisons are shown for a) all sites, b) sites associated with phenotype in the HmB locus, and c) sites associated with phenotype in the HmYb locus.

Figure 5. Mean and standard error of relative gene expression for Mad, Slu7, Kinesin, and GPCR.

For the ‘Between wing segments’ experiment, H. melpomene cythera forewings were dissected into P (proximal), B (band) and D (distal) segments of three pupal developmental stages; EP (early pupae), PO (pre-omochrome) and OO (only omochrome). For the “Between Races” experiment, forewings (Fw) and hindwings (Hw) of H. m. malleti (left side) and H. m. cythera (right side), were dissected from Larvae (5^th instar) and EP (early pupae).