A widely diverged locus involved in locomotor adaptation in Heliconius butterflies

Abstract

Heliconius butterflies have undergone adaptive radiation and therefore serve as an excellent system for exploring the continuum of speciation and adaptive evolution. However, there is a long-lasting paradox between their convergent mimetic wing patterns and rapid divergence in speciation. Here, we characterize a locus that consistently displays high divergence among Heliconius butterflies and acts as an introgression hotspot. We further show that this locus contains multiple genes related to locomotion and conserved in Lepidoptera. In light of these findings, we consider that locomotion traits may be under selection, and if these are heritable traits that are selected for, then they might act as species barriers.

Fig. 1. The geographical distribution and genome-wide phylogeny of Heliconius butterflies.

(A) The distributions of sampled subspecies/species are shown in different colors, along with wing pattern images. (B) The maximum likelihood phylogenetic tree was constructed on the basis of 102.5 Mb of genome-wide SNP data. The color ranges represent major clades. The scale bar represents the percentage of substitutions per site. The round labels indicate the selected subspecies used in downstream analyses, among which the three H. melpomene subspecies form co-mimicry pairs with selected H. erato subspecies and H. timareta subspecies.

Fig. 2. Genomic divergence during early speciation.

Pairwise F_ST values are calculated across the genome in 50-kb sliding windows with a step size of 20 kb, between subspecies of H. timareta, H. cydno and H. pachinus, and H. melpomene; between the eastern and western populations of H. melpomene; and between H. m. rosina and H. c. galanthus. Twenty-kilobase sliding windows are used for the calculation of F_ST values between subspecies of H. erato. Results of comparisons within each clade are grouped in the same color. The red dashed lines represent the top 1% threshold across the genome. The L locus and other wing patterning loci are labeled with gray bars. In H. erato, the orthologous loci of Yb and B/D are known as Cr and R/D, respectively.

Fig. 3. Local patterns of LD, genetic diversity, and phylogenetic trees along the L locus.

(A) Pairwise LD, measured as r², is estimated among biallelic SNPs along the L locus and its adjacent regions for all the H. timareta samples, with the L locus indicated by bold black bars. (B) Values of nucleotide diversity (π) are calculated along the L locus for the subspecies of H. timareta, H. cydno and H. pachinus, and H. melpomene. (C) Maximum likelihood phylogenetic trees are constructed for every 20-kb window along the L locus and adjacent regions. The H. timareta subspecies are grouped polyphyletically in the window of topology B, with signatures of introgression between H. t. thelxinoe and H. numata and between H. t. timareta and H. m. malleti. The eastern and western melpomene populations join the silvaniform and cydno-timareta clades, respectively, in the window of topology A, with signatures of introgression between eastern melpomene and H. elevatus and between western melpomene and H. pachinus, respectively.

Fig. 4. Correlations between genetic differentiation and geographical distance in H. cydno subspecies.

Point plots show the correlation of pairwise genetic differentiation and geographical distance in H. cydno subspecies. The distribution center of each H. cydno subspecies was measured as the geographic midpoint of sampled individuals, and pairwise geodesic distances between H. cydno subspecies were calculated. F_ST for specific loci, L (A) and two wing patterning loci Yb (C) and Br (D), and genome-wide F_ST (B) are plotted against geographical distances. F_ST of loci L and Yb are significantly and positively correlated with distances, as well as the genome-wide mean F_ST. The differentiation of L shows the strongest evidence for differentiation by distance with a slope eightfold greater than the genome-wide mean F_ST and twofold greater than F_ST of Yb in the linear fitting. **P < 0.01; *P < 0.05.

Fig. 5. Functional characterization of the L locus in Heliconius butterflies.

(A) Boxplots show the WBFs of H. c. galanthus, H. m. malleti, and H. m. rosina. Center line indicates median; box extends to upper and lower quartile, and whisker extends to the minimum and the maximum. ***P < 0.001; *P < 0.05 (Tukey’s test). H. cydno beat wings significantly slower than H. melpomene in free-moving flight [one-way analysis of variance (ANOVA), F_2,26 = 3.37, P < 0.001]. The table summarizes the differences of WBF and the genetic differentiations at the L locus between pairs of races. (B) Expression levels of Col4α1, Vkg, Oseg4, and na in western and eastern melpomene races. Col4α1 and Vkg have significantly lower expression in western melpomene in comparison with that of eastern melpomene. The boxplots summarize the normalized counts (DESeq2’s median of ratios). ***P < 0.001, **P < 0.01, and *P < 0.05. (C) Tissue-specific expression levels of Col4α1, Vkg, Oseg4, and na in H. m. rosina. Boxplots show the ranges of fragments per kilobase of transcript per million mapped reads (FPKM). Different letters (a, b, c, and d) indicate statistically significant differences across groups (P < 0.05; Scheffe’s test). FW, forewing; HW, hindwing.

Fig. 6. Functional implication of the orthologous genes of the L locus in D. melanogaster.

(A to C′) The adult Drosophila scutellum, showing two pairs of mechanosensory bristles, is required for wing-wing coordination during flight (A). The adult wing blade contains five longitudinal veins (L1 to L5) (B) and one dorsal anterior marginal vein [(B) and enlarged in (B′)], embedded with tracheal tubes for air exchange. Notably, wing mechano- and gustatory (chemosensory) sensilla, specified by Senseless at the third-instar imaginal discs (C and C′), are closely associated with the marginal vein (B′). The expression of w RNAi by ap-Gal4 does not produce any obvious defect (A to C′). (D to O′) Knocking down either Col4α1, Vkg, na, or Oseg4 produces similar defects associated with locomotion. The scutellum structure is altered, and the number of mechanosensory bristles is reduced (D, G, J, and M). Ectopic veins, hairs (E, H, K, and N) (blue arrowheads), and sensory bristles [insets in (E) and (H)] are produced on wing blade. The pattern of gustatory bristles is misarranged (E′, H′, K′, and N′) (red arrows), mechanosensory bristles are occasionally lost (K′ and N′) (purple arrows), and Senseless is altered consistently (F′, I′, L′, and O′) (yellow arrows). The effects of knocking down the four candidate genes on flight ability (N = 80) are shown in (P) and (Q). The KD flies fail to fly in response to the dropping stimulus (P). The ascent speed of KD flies is slower than that of non-KD individuals (Q). Different letters (a and b) indicate statistically significant differences across groups (P < 0.05) (P and Q). Scale bars, 100 μm.

Fig. 7. Multiple alignments of the L locus.

The multiple alignment of the L locus between several insect genomes suggests that it is a conserved genomic region only in Lepidoptera. Locally collinear blocks are shown as rectangular blocks, with four orthologous genes labeled by horizontal bars in different colors. The interchromosomal boundaries are indicated by red vertical bars.