Association mapping of colour variation in a butterfly provides evidence that a supergene locks together a cluster of adaptive loci

Abstract

Supergenes are genetic architectures associated with discrete and concerted variation in multiple traits. It has long been suggested that supergenes control these complex polymorphisms by suppressing recombination between sets of coadapted genes. However, because recombination suppression hinders the dissociation of the individual effects of genes within supergenes, there is still little evidence that supergenes evolve by tightening linkage between coadapted genes. Here, combining a landmark-free phenotyping algorithm with multivariate genome-wide association studies, we dissected the genetic basis of wing pattern variation in the butterfly Heliconius numata. We show that the supergene controlling the striking wing pattern polymorphism displayed by this species contains several independent loci associated with different features of wing patterns. The three chromosomal inversions of this supergene suppress recombination between these loci, supporting the hypothesis that they may have evolved because they captured beneficial combinations of alleles. Some of these loci are, however, associated with colour variations only in a subset of morphs where the phenotype is controlled by derived inversion forms, indicating that they were recruited after the formation of the inversions. Our study shows that supergenes and clusters of adaptive loci in general may form via the evolution of chromosomal rearrangements suppressing recombination between co-adapted loci but also via the subsequent recruitment of linked adaptive mutations. This article is part of the theme issue 'Genomic architecture of supergenes: causes and evolutionary consequences'.

Keywords: association study; cluster of adaptive loci; divergence hitchhiking; inversion; multivariate association; wing colour pattern.

Figure 1.

Genetic architecture and wing pattern diversity in H. numata. (a) Genetic architecture of the H. numata mimicry supergene P characterized by three polymorphic inversions of respective size 400 kb, 200 kb and 1150 kb. (b) Schematic diversity of wing patterns of H. numata in our dataset. (c) Two-dimensional approximation of the morphological space representing the phenotype diversity observed in H. numata. The dotplot displays results from a principal component analysis (PCA) (the first two components are displayed here) computed on wing pattern variations as obtained using colour pattern modelling (CPM) [32]. For display purposes, butterflies were manually classified into mimetic forms based on the literature [33]; different forms are depicted by different colours. The butterflies sampled for this study represent the commonest forms observed in H. numata. Different supergene genotypes are depicted by different symbol shapes. Results for PC 3 and PC 4 are presented in the electronic supplementary material, figure S8. PCAs computed on specimens with the same supergene arrangement and on specific parts of wing pattern are presented in the electronic supplementary material, figure S9. (Online version in colour.)

Figure 2.

Genome-wide association of genetic and wing pattern variation. (a) Multivariate association study using as phenotype the first six principal components describing wing pattern variations (presented in figure 1c and electronic supplementary material, figure S8) and all specimens regardless of their genotype at the supergene. One major peak of association in noticeable on chromosome 15 corresponding to the supergene. One minor peak can be seen on chromosome 7. This is due to an assembly error in the reference genome Hmel2, and in reality, this region lies within the supergene region (see §4). (b) Focus on the peak of association on chromosome 15, corresponding to the position of the three polymorphic chromosomal inversions P₁, P₂ and P₃. (Online version in colour.)

Figure 3.

Distinct regions within the supergene are associated with variation in wing pattern features. (a,d) Multivariate association studies computed on Hn123 and Hn0 specimens, respectively, on the entire wing pattern variation (hindwing and forewing together, here using four principal component as multivariate phenotypes). The plotted p-value is the statistical p-value from the multivariate test of association. 1 × 10⁻⁶ permutations were performed for each variant. Variants highlighted in orange are variants with an empirical p-value < 1 × 10⁻⁶ (i.e. for which no permutation resulted in a lower statistical p-value). The positions of inversion breakpoints are represented by the dotted vertical lines. Associations computed with different numbers of phenotypic principal components are presented in electronic supplementary material, figures S14 and S15. (b,e) Density of significantly associated variants in multivariate analyses (with empirical p-value < 1 × 10⁻⁶) along the chromosome 15. Analyses were computed in 10 000 bp overlapping sliding window (with 100 bp overlap). All significantly associated variants in one or more of the multivariate association analyses (using 2, 3, 4, … , 6 phenotypic principal components; electronic supplementary material, figures S14 and S15) were used. (c,f) Density of significantly associated variants in univariate analyses (with empirical p-value < 1 × 10⁻⁶) along chromosome 15. Analyses computed in 10 000 bp overlapping sliding windows (with 100 bp overlap). All significantly associated variants in one or more of the univariate association analyses (focusing on different part of the wing and using the first, second or third using phenotypic principal component as phenotype; electronic supplementary material, figures S16–S19) were used. (g,h) Phenotypic effects of the top variant from each of the 15 regions that displayed a clear enrichment in significantly associated variants ((b–f) coloured arrows) to the wing pattern in Hn123 or Hn0 specimens, respectively. Heatmaps from blue to red represent, for every pixel, the strength and direction of association of the derived allele, that is ho allelic change at a given genetic position affects this wing area. Overall effects are shown as well as colour-specific effects, the latter representing the extent to which allelic change is associated with the presence or absence of each colour at this wing area. Because blue and red represent the direction of the association, opposite directions (i.e. red and blue values) in the same wing area in two colour-specific heatmaps indicate that the focal locus is associated with a change from one colour to the other in this area. For instance, if the effect of a genetic variant on a given wing area is highlighted in blue when looking at the orange pattern but in red when looking at black pattern, that means that change at this variant is associated with a switch from orange to black at this wing area.

Figure 4.

The effect of four selected variants on Hn123 wing pattern variation. Representation of the association of some genetic variants with specific wing pattern variation. See figure 3g,h and electronic supplementary material, figures S14–19 for additional representations of the association with specific aspect of wing patterns. The first principal components of analyses computed on different parts of the wing are used as proxy for the phenotype (y-axis): PCA computed on hindwings only (a), on the middle part of the forewings only (b) and on the tips of the forewings only (c). See electronic supplementary material, figure S9 for another representation of the principal components. Each dot is an individual specimen. Instead of annotating y-axis with eigenvalues (values of individual on principal components), schematic butterflies with average phenotypes along the principal components are displayed. Boxplot elements: central line, median; box limits, 25th and 75th percentiles; whiskers, 1.5 × interquartile range. (a) Effect of the most strongly associated variants in regions 6 (cortex, intron 2) and 2 (wash, intron 2) on the amount of black on hindwings. (b) Effect of the most strongly associated variants in region 1 (intergenic HMEL032679-cortex) on the forewing middle part. (c) Effect of the most strongly associated variants in region 9 (intergenic HMEL022251-HMEL032698) on the tip of the forewings. (Online version in colour.)