Mimicry diversification in Papilio dardanus via a genomic inversion in the regulatory region of engrailed- invected

Abstract

Polymorphic Batesian mimics exhibit multiple protective morphs that each mimic a different noxious model. Here, we study the genomic transitions leading to the evolution of different mimetic wing patterns in the polymorphic Mocker Swallowtail Papilio dardanus. We generated a draft genome (231 Mb over 30 chromosomes) and re-sequenced individuals of three morphs. Genome-wide single nucleotide polymorphism (SNP) analysis revealed elevated linkage disequilibrium and divergence between morphs in the regulatory region of engrailed, a developmental gene previously implicated in the mimicry switch. The diverged region exhibits a discrete chromosomal inversion (of 40 kb) relative to the ancestral orientation that is associated with the cenea morph, but not with the bottom-recessive hippocoonides morph or with non-mimetic allopatric populations. The functional role of this inversion in the expression of the novel phenotype is currently unknown, but by preventing recombination, it allows the stable inheritance of divergent alleles enabling geographic spread and local coexistence of multiple adaptive morphs.

Keywords: Batesian mimicry; butterflies; genomic rearrangement; polymorphism; supergene.

Figure 1.

Phenotypic variation in Papilio dardanus and samples used. (a) Seven female forms and a male. (b) Origin of samples for sequencing and population genetic analyses, from four subspecies: P. dardanus polytrophus (Kenya), P. dardanus tibullus (Kenya), P. dardanus cenea (South Africa), and P. dardanus meriones (Madagascar). The specimen of subspecies P. dardanus tibullus was used for the construction of the draft genome sequence. The tree depicts the relationships among these four subspecies and is based on a tree presented in [3]. Three female forms were analysed: hippocoonides, cenea, and ‘male-like’.

Figure 2.

Population genomic analysis of the full engrailed–invected containing scaffold. The exons of orange and engrailed are shown by large arrows and the upstream region marked in grey. Brood 59 and Brood 48: recombination events in pedigree broods. SNPs from RADseq data for two broods are mapped on the engrailed–invected scaffold, shown by black triangles. Red triangles mark the intervals with confirmed recombination events, and the number of recombination events within these intervals are circled. The central band in the figure shows the map of the scaffold with exons (thin vertical red lines) and the upstream region of engrailed (grey). ACT and PD indicate the position of the AFLP markers of [16]. Linkage disequilibrium (LD; Kelly's ZnS statistic), F_st, nucleotide diversity (pi), coverage and p-distance (to the reference genome) for the scaffold, calculated for the f. cenea and f. hippocoonides samples in 5 kb windows. Coverage and p-distance were calculated separately for the four cenea and for four hippocoonides specimens. The p-distance to the reference genome is also given for the P. dardanus meriones sample. Scales are in million base pairs. (Online version in colour.)

Figure 3.

Length and relative position of the inversions in the upstream regulatory region of engrailed. At the top, the map of the engrailed–invected region is shown, with short arrows indicating exons and miR-2768 [39] shown in blue. Below the map is the direction of boundary-defining primers. The grey shading indicates the extent of the 40 kb inversion associated with f. cenea. For each of the forms, dark grey–light grey shading is used to indicate directionality of the 40 kb region. Scale is in base pair. (Online version in colour.)