The genomics and evolution of inter-sexual mimicry and female-limited polymorphisms in damselflies

Abstract

Sex-limited morphs can provide profound insights into the evolution and genomic architecture of complex phenotypes. Inter-sexual mimicry is one particular type of sex-limited polymorphism in which a novel morph resembles the opposite sex. While inter-sexual mimics are known in both sexes and a diverse range of animals, their evolutionary origin is poorly understood. Here, we investigated the genomic basis of female-limited morphs and male mimicry in the common bluetail damselfly. Differential gene expression between morphs has been documented in damselflies, but no causal locus has been previously identified. We found that male mimicry originated in an ancestrally sexually dimorphic lineage in association with multiple structural changes, probably driven by transposable element activity. These changes resulted in ~900 kb of novel genomic content that is partly shared by male mimics in a close relative, indicating that male mimicry is a trans-species polymorphism. More recently, a third morph originated following the translocation of part of the male-mimicry sequence into a genomic position ~3.5 mb apart. We provide evidence of balancing selection maintaining male mimicry, in line with previous field population studies. Our results underscore how structural variants affecting a handful of potentially regulatory genes and morph-specific genes can give rise to novel and complex phenotypic polymorphisms.

Fig. 1. The evolution of female-limited colour polymorphisms in Ischnura damselflies.

a, A previous phylogenetic study and ancestral state reconstruction proposed that the genus Ischnura had a sexually dimorphic ancestor, with O-like females (red circle). The O morph is markedly different from males, having a bronze-brown thorax and faint stripes, instead of the black thoracic stripes on a bright blue background of males. b, Male mimicry (A females, blue circle) has evolved more than once, for instance, in an ancestor of the (expanded) clade that includes the common bluetail (I. elegans, outlined with solid line) and the tropical bluetail (I. senegalensis, outlined with dashed line). c, I. elegans is a trimorphic species, due to the recent evolution of a third female morph, I (yellow circle). In I. elegans, morph inheritance follows a dominance hierarchy, where the most dominant allele produces the A morph and two copies of the most recessive allele are required for the development of O females. In contrast, the O allele is dominant in I. senegalensis. Terminal nodes in the phylogeny represent different species. Grey triangles represent other clades of Ischnura, which are collapsed for clarity. Damselfly images adapted from ref. under a Creative Commons licence CC BY 4.0. Source data

Fig. 2. Morph determination in I. elegans is controlled in a ~1.5 mb region of chromosome 13.

a, SNP-based GWAS in all pairwise analyses between morphs. Genomic DNA from 19 wild-caught females of each colour morph and of unknown genotype was extracted and sequenced for these analyses. Illumina short reads were aligned against an A morph genome assembly, generated from nanopore long-read data (Extended Data Fig. 1). b, A closer look at the SNP associations on the unlocalized scaffold 2 of chromosome 13, which contained all highly significant SNPs. Transcripts expressed in at least one adult of both I. elegans and I. senegalensis are shown at the bottom (see also Fig. 6). Grey transcripts are shared by all morphs, whereas blue transcripts are uniquely present in A or A and I samples (see ‘Shared and morph-specific genes reside in the morph locus’). The y axis in a and b indicates unadjusted −log₁₀ P values calculated from chi-squared tests. c, F_ST values averaged across 30 kb windows for the same pairwise comparisons as in the SNP-based GWAS. The dashed line marks the 95th percentile of all non-zero F_ST values across the entire genome. The red double arrow shows the region of elevated divergence between O and both A and I samples (∼50 kb–0.2 mb). The blue double arrow shows the region of elevated divergence between A and both O and I samples (∼0.2 mb–1.5 mb). d,e, Population-level estimates of Tajima’s D (d) and π (e) averaged across 30 kb windows. The shaded area contains the 5th–95th percentile of all genome-wide estimates. Source data

Fig. 3. Female morphs of I. elegans differ in genomic content.

a,b, Number of significant k-mers (below the 5% false-positive threshold; Methods) associated with pairwise genome-wide analyses and mapped to the unlocalized scaffold 2 of chromosome 13, in the A-morph assembly (a) and the I-morph assembly (b). c,d, Standardized read depths along the unlocalized scaffold 2 of chromosome 13, relative to background coverage of the A-morph assembly (c) and the I-morph assembly (d). Solid lines (yellow, blue and red) show short-read data (19 samples per morph) and black dashed lines represent long-read data (1 sample per morph).Grey areas show regions of genomic content present in A and I individuals, but absent in all but one O sample. Note that different regions of the scaffold are plotted for the two assemblies (see main text). Source data

Fig. 4. SVs differentiate morph haplotypes in the common bluetail damselfly (I. elegans).

a, Alignment between morph-specific genomes assembled from long-read nanopore samples with genotypes Ao, Io and oo. Grey lines represent alignments of at least 5 kb and >70% identity. The black line connects regions of genomic content shared by the three morphs within the morph locus. The red to blue gradient represents a ~20 kb region that carries an inversion signature in A and I females relative to the O haplotype (Extended Data Fig. 2). The blue to yellow gradient represents a ~150 kb alignment between the start of the unlocalized scaffold 2 of chromosome 13 in A and a region ~3.5 mb apart in the I haplotype. Coordinates at the bottom are based on the DToL reference assembly. b, Schematic illustration of the hypothetical sequence of events responsible for the evolution of novel female morphs. First, a sequence originally present in O was duplicated and inverted in tandem, potentially causing the initial divergence of the A allele. Second, part of this inversion was subsequently duplicated in A, in association with a putative TE, leading to multiple inversion signatures in the A haplotype relative to an O reference (Extended Data Fig. 3). Finally, part of the A duplications were translocated into a position ~3.5 mb downstream into an O background, giving rise to the I allele. Currently, A females are also characterized by another region of unique content and unknown origin (question mark). A females show elevated sequence divergence in the internal region of the morph locus that is shared by all haplotypes (dark grey bars; see also black line in a). Coordinates on the O haplotype are based on the (DToL) reference assembly. Grey numbers in IV give the approximate size of genomic sequences in A and I that are absent in O. Damselfly images adapted from ref. under a Creative Commons licence CC BY 4.0. Source data

Fig. 5. A shared genomic basis of A females in I. elegans and I. senegalensis.

a, I. senegalensis is a female-dimorphic species, where one female morph (O-like) is distinctly different from males and resembles O females in I. elegans, and the other female morph (A) is a male mimic. Photo credit: Mike Hooper. b, Standardized read depth of pool-seq samples (n = 30 females of each morph per pool) of I. senegalensis, against the A-morph assembly of I. elegans, calculated in 500 bp windows. The x-axis shows the first 1.5 mb of the unlocalized scaffold 2 of chromosome 13. c, Alignments between morph-specific genomes from a homozygous O-like female of I. senegalensis (top), an Ao female of I. elegans (middle) and a homozygous A female of I. senegalensis (bottom). Lines connecting the assemblies represent alignments of at least 500 bp and >70% identity. The black line connects genomic content in the morph locus, which is shared by the three morphs of I. elegans. In I. elegans, this region is rich in SNPs differentiating A females from the other two morphs (see Fig. 2b). The blue–turquoise gradient connects sequences uniquely present in the A morphs of I. elegans and I. senegalensis. Source data

Fig. 6. The morph locus of I. elegans is situated in the unlocalized scaffold 2 of chromosome 13.

a, Diagram of the ~1.5 mb morph locus on the A-morph assembly, showing from top to bottom: morph-specific read depth coverage; the location of LINE retrotransposons in the the Jockey family; the mapping locations of A-derived reads with a previously detected inversion signature against O females; and transcripts expressed in at least one adult individual of both I. elegans and I. senegalensis. Transcripts plotted in black are present in both the A and O assemblies, while transcripts in blue are located in genomic regions that are unique to the A haplotype or are shared between A and I but not the O allele. b, Functional annotations and sex- and morph-specific expression of transcripts. Square fill indicates whether transcript expression was detected in each group. RNA-seq data for I. elegans comes from whole-thorax samples from sexually immature and sexually mature wild-caught adults (n = 3 females of each morph and 3 males). RNA-seq data for I. senegalensis comes from a recent study in which the abdomen, head, thorax and wings were sampled in two females of each morph and two males (one individual of each group sampled upon emergence and one sampled after two days). Source data

Extended Data Fig. 1. Outline of data and analyses used in this study.

For our main study species Ischnura elegans, we obtained short-read genomic data from 19 field-caught females per morph, and long-read genomic data from three females with genotypes Ao, Io, and oo. The long-read samples were used to assemble morph-specific genomes, scaffolded against the Darwin Tree of Life reference assembly. We obtained whole-thorax RNAseq data from females of each morph in both sexually immature and sexually mature colour phases (n = 3 of each morph and colour phase). Immature and mature males (n = 3 of each) were also sampled for whole-thorax RNAseq data. We used short-read pool-seq data (n = 30 individuals of each morph per pool) of the close relative Ischnura senegalensis to investigate whether the female polymorphisms in both species share a genomic basis. We also analysed expression levels of candidate genes in this species, using samples from a previously published study, which produced transcriptomic data from four body parts (head, thorax, wing and abdomen) of each A females, O females and males (n = 2), sampled at adult emergence and two days thereafter. The k-mer based GWAS is reference-free, but significant k-mers were mapped to the morph-specific assemblies to determine their chromosomal context. Damselfly images adapted from ref. under a Creative Commons licence CC BY 4.0.

Extended Data Fig. 2. An inversion signature differentiates A and I individuals from the O morph.

Read mapping and sample coverage at the start of the scaffold 2 of chromosome 13 in a our O assembly and b the DToL reference assembly, showing a signature of a ~ 20 kb inversion in A and I samples. A single O sample also exhibited this signature but was excluded here for clarity (see Supporting Text 3). Note that the first 415 kb of the reference DToL assembly are absent in our scaffolded O assembly, and therefore the x-axis is shifted by 415 kb in b.

Extended Data Fig. 3. The A and I reads mapped to inversion break points on the O assembly (see Extended Data Fig. 2) map to multiple locations on the A assembly.

a Reads from the first inversion breakpoint. b Reads from the second inversion breakpoint. Each row represents a sample and each circle an individual read. The x-axis corresponds to coordinates on the A assembly. Source data

Extended Data Fig. 4. Proportion of TE content in non-overlapping 1.5 mb regions.

The gray dots correspond to genomic windows outside chromosome 13. The main assembly and the unlocalized scaffolds of chromosome 13 are depicted with different colours. The dashed line marks the 95 percentile of TE coverage across all windows. Source data

Extended Data Fig. 5. Linkage disequilibrium (LD) in the genome of Ischnura elegans.

LD estimates are shown for the first 15 mb of each chromosome and all unlocalized scaffolds of chromosome 13. The morph locus is found in the first ~ 1.5 mb of the unlocalized scaffold 2 of chromosome 13, which has a total size of ~ 15 mb. Each dot represent the square correlation coefficient (R²) between two variant sites on the x axis, separated by the number of sites indicated in the y axis.

Extended Data Fig. 6. Evidence of a translocation between A and I haplotypes.

Mapping and coverage of long reads from an Io sample across the first 5.6 mb of the unlocalized scaffold 2 of chromosome 13 in the A assembly, showing a signature consistent with either a 5.54 mb inversion or a translocation of inverted A content. Absence of morph divergence beyond ~1.5 mb on the A assembly supports the translocation scenario.

Extended Data Fig. 7. Structural variants between A and O-like females of I. senegalensis along the morph locus identified in I. elegans.

a Read mapping and sample coverage of I. senengalensis pool-seq data at the start of the unlocalized scaffold 2 of chromosome 13 in the O assembly of I. elegans. The same ~ 20 kb inversion signature is found in A and I samples of I. elegans (see Extended Data Fig. 2). b-c The A-pool reads mapped to break points on the O assembly map to multiple locations on the A assembly. b Reads from the first breakpoint. c Reads from the second breakpoint. Each row represents a pool of I. senegalensis and each circle an individual read. The x-axis corresponds to the A assembly of I. elegans. Source data

Extended Data Fig. 8. Morph divergence using the DToL assembly (O haplotype) as mapping reference.

a SNP-based genome-wide associations in all pairwise analyses between morphs. Genomic DNA from 19 wild-caught females of each colour morph and of unknown genotype was extracted and sequenced for these analyses. Illumina short reads were aligned against the DToL reference assembly. b A closer look of the SNP associations on the unlocalized scaffold 2 of chromosome 13, which contained all highly significant SNPs. The y axis in a and b indicates unadjusted -Log₁₀ P-values calculated from chi-squared tests. c Fst values averaged across 30 kb windows for the same pairwise comparisons as in the SNP based GWAS. The dashed line marks the 95 percentile of all non-zero Fst values across the entire genome. The red double arrow shows the region of elevated divergence between A and both O and I samples. Population-level estimates of d Tajima′s D, and e nucleotide diversity (π) averaged across 30 kb windows. The shaded area contains the 5–95 percentile of all genome-wide estimates. Source data