Colour polymorphisms are popular study systems among biologists interested in evolutionary dynamics, genomics, sexual selection and sexual conflict. In many damselflies, such as in the globally distributed genus Ischnura (forktails), female colour polymorphisms occur in some species. Female-polymorphic species contain two or three female morphs, one of which is male-coloured (androchrome or male mimic) and co-exists with sexually dimorphic (heterochrome) females. These female colour polymorphisms are considered to be maintained by frequency-dependent sexual conflict, but their macroevolutionary histories are unknown, due to the lack of a robust molecular phylogeny. Here, we present the first time-calibrated phylogeny of Ischnura, using a multispecies coalescent approach (StarBEAST2), incorporating both molecular data and fossil information of 41 extant species (55% of the genus). We estimate the age of Ischnura to be between 13.8 and 23.4 millions of years, i.e. Miocene. We infer the ancestral state of this genus as female monomorphism with heterochrome females, with multiple gains and losses of polymorphisms, evidence of trans-species polymorphisms and a significant positive relationship between polymorphism incidence and current geographic range size. Our study provides a robust phylogenetic framework for future research on the dynamic macroevolutionary history of this clade with its extraordinary diversity of sex-limited female polymorphisms.

General data collection procedures

We have gathered DNA-sequence data, phenotypic data and geographic range size data for 41 species of Ischnura. We use the generic names Ischnura and Pacificagrion, according to the “World Odonata List” (https://www.pugetsound.edu/academics/academic-resources/slater-museum/biodiversity-resources/dragonflies/world-odonata-list2/), which maintains a taxonomic classification of Odonata with currently valid names.

We obtained sequences for five molecular markers including nuclear (H3: 28.86 kb, PMRT: 51.68 kb, 28S: 67.91 kb) and mitochondrial regions (COI: 57.72 kb, 16S: 47.64 kb) from 147 specimens, of which 64 were sequenced for this study (Table S1). Genomic DNA was extracted from 1-2 legs using a Qiagen DNeasy tissue kit and following the manufacturer’s instructions, except that we increased the amount of proteinase K and incubation time for museum samples collected before 2005, in order to obtain a higher DNA yield.  Amplification products were purified using the Affymetrix ExoSAP-IT reagent (Cat. No. 78200) and sequenced on an ABI 3100 capillary sequencer at Lund University. These novel sequence data were complemented with already-published sequences (Table S1), downloaded from the National Centre for Biotechnology Information (https://www.ncbi.nlm.nih.gov/genbank/). There were no indels in the sequences of two molecular markers (COI and H3), which were aligned using MUSCLE. To reduce gap overmatching in ribosomal sequences (16S and 28S) and intronic regions (present in PMRT), we used the phylogenetically informed algorithm PRANK to align these data.

Data on female colour morph status (i.e. female-monomorphic, or female-polymorphic) has been recently compiled for Ischnura damselflies, based on published literature, field observations, examination of museum specimens and online resources. We complemented these data with information on whether female-polymorphic species had two or three morphs (i.e. female-dimorphic or female-trimorphic) (Table S2). Here, species with a single female morph were classified as female-monomorphic with androchrome females,  if they had only one female morph with colouration and patterning similar to that of conspecific males, or as female-monomorphic with heterochrome females, if they had only one female morph with colouration and patterning differing from the conspecific males. If a species contained two distinct female morphs (always one androchrome and one heterochrome), we classified it as female-dimorphic. Finally, species that had three morphs (which were always one androchrome and two heterochrome female morphs), were classified as female-trimorphic.

We complemented these colour-state data with coarse-grained estimates of geographic range size (Table S2). Because accurate geographic distribution data is lacking for many species, we estimated geographic range using occurrence data from administrative regions, mainly countries, and where available, states or provinces (Table S2). Presence data was translated to geographic range size (km²) by adding the areas of all the countries or other administrative regions where a species is known to occur (Table S2), as has been done previously in comparative studies where geographic range is taken as a proxy for population size. The areas of administrative regions were extracted from the GADM database of global administrative areas.

Time-calibrated species tree

We analysed sequence data under a multispecies coalescent model using StarBEAST2. As the model assumes free recombination amongst unlinked loci, we linked the tree models between mitochondrial sequences. We placed a birth-death prior, conditional on the root, on the species tree to allow for extinction, set the gene ploidy to 0.5 for mitochondrial sequences and 2.0 for nuclear sequences and allowed population-size to be analytically integrated to avoid over-parameterization. We used a relaxed uncorrelated log normal (UCLN) clock per locus to allow for variation in evolutionary rates between lineages. Substitution models for each locus were inferred during the analysis using the package “bModelTest” to allow for site model uncertainty. We placed a lognormal prior with a mean of 0.0115 and standard deviation of 1.0 on each relaxed mitochondrial clock and a broader lognormal prior with a mean of 0.0115 and standard deviation of 2.0 on each relaxed nuclear clock to centre clock rates on the suggested insect mitochondrial divergence rate of 2.3% per million years whilst allowing for clock rate uncertainty, particularly in the nuclear loci.

Data for Ischnura fossils were downloaded from the Paleobiology database (www.pbdb.org) on 14^th May 2018. There were three fossils classified as either in or closely related to the genus Ischnura. Two of these fossils were relatively young (7.2 - 5.3 mya), each one consisting of a single wing with missing fragments. The third fossil was preserved in Dominican amber (20.44 – 13.82 mya) and was identified and described as Ischnura velteni based on a combination of wing venation and morphological traits in other body parts. The combination of these characters suggests that I. velteni belongs to the extant genus Ischnura (Bechly, 2000), and we therefore considered it the most appropriate crown fossil with which to calibrate the root age of the Ischnura tree. We placed an exponentially distributed prior on the root age of the species tree with an offset of 13.8 to ensure the most recent common ancestor (MRCA) of Ischnura could not have appeared after the minimum age of this known member of the genus. The mean of the exponential was drawn from a gamma hyperprior with a shape parameter α = 5.0 and an inverse scale parameter β = 2.0. This hyperprior specification centred the root age prior distribution around 10 mya prior to the offset, allowing a moderately high level of uncertainty in the age difference between the fossil and the MRCA of Ischnura.

We ran two separate Markov Chain Monte Carlo (MCMC) algorithms for 200 million iterations with a sampling rate of 20,000. We confirmed that independent chains were stationary and had converged for these and all subsequent analyses in R v.3.6.1 (R Core Team 2019) using the “coda” package. We combined posterior trees using LogCombiner (available as part of BEAST v2.4.7, simultaneously discarding 20% of the posterior as burnin, and summarized them into a maximum clade credibility (MCC) tree using TreeAnnotator (available as part of BEAST v2.4.7). The MCC tree and all further phylogenetic plots were annotated in R v.3.6.1 (R Core Team 2019) using the packages “ape", “ggtree” and “ggplot2”.

Reconstruction of ancestral character states of female colour polymorphism

We used the Multistate-method in BayesTraits version 3.0 to infer female-colour character states at ancestral nodes, given this new Ischnura phylogeny and the female-colour character states of extant species. To account for phylogenetic uncertainty, we ran the model using a sample of 1000 trees from the combined species-tree posterior of the two independent runs. Additionally, we ran the model using reversible jump MCMC to account for uncertainty in the number of model parameters. In this model, the number of different transition rate parameters and their values are sampled in proportion to their posterior probability. Consequently, in each posterior sample a number of transition rates are set to zero, and the non-zero rates belong to one or more categories. We placed an exponential prior on all rate parameters, with the mean seeded from a uniform hyperprior bounded between 1 and 100. We ran two separate chains for 10 million iterations, sampling every 1000 with 20% posterior burnin. Trees were scaled to an average branch length of 0.01 to prevent various rates from becoming too small. Although our main focus is on ancestral state reconstructions, we present the posterior distributions of transition rate parameters in the Supporting Material (Fig. S1).

Linking female colour polymorphism to geographic range size

We investigated if the occurrence of female colour polymorphisms could be influenced by demographic changes by using geographic range size of extant species as a proxy for population size. We used a Bayesian phylogenetic mixed-effect model (BPMM) implemented in the package MCMCglmm v. 2.29 to test for such phylogenetic correlation. Female-morph status was treated as a binary response variable, with species being classified as either female-monomorphic or female-polymorphic. Since the occurrences of female-trimorphic species and female-monomorphic species with androchrome females were low, we classified these species-states as being female-polymorphic and female-monomorphic, respectively. We thus investigated if the probability of a species being female-polymorphic were significantly and positively related to geographic range size (a continuous predictor variable).

We accounted for phylogenetic uncertainty by using a random sample of 2100 species-trees generated in StarBEAST2. To estimate the phylogenetic random effect, each tree was sampled for 10000 iterations, of which only the last iteration was saved to the posterior, and the first 100 were also discarded as burnin. We used a Knonecker prior (mu=0, V= σ2units + π2/3) for the fixed effects. We used a parameter-expanded χ² distribution (V=1 nu=1000, alpha.mu=0, alpha.V = 1) for the prior on the phylogenetic variance components.. The residual variance was fixed to 1, as it cannot be identified for binary response models. As a statistical test of whether the probability of female colour polymorphism increases with range size, we report the fraction of the posterior distribution in which increasing range size does not result in an increasing probability of female colour polymorphism. The results of this analysis were visualized by plotting the expected probability of female colour polymorphism against estimated geographic range size.

Data consists of raw data (phenotypic information, DNA-sequences) collected in the field by us, through museum collections and the litterature. These data should allow for the replication of ourselves and should serve as a baseline for phylogenetic comparative studies of this genus.

For additional details, contact the corresponding author (erik.svensson@biol.lu.se).