Spatiotemporal Dynamics of Non-Ecological Speciation in Rubyspot Damselflies (Hetaerina spp.)

Abstract

Non-ecological speciation is a common mode of speciation, which occurs when allopatric lineages diverge in the absence of pronounced ecological differences. Yet, relative to other speciation mechanisms, non-ecological speciation remains understudied. Numerous damselfly clades are characterised as non-adaptive radiations (the result of several rounds of non-ecological speciation without subsequent divergence), but there are few damselfly lineages for which we have a detailed understanding of the spatiotemporal dynamics of divergence. Recent phylogeographic analyses demonstrate that American rubyspot damselflies (Hetaerina americana sensu lato) actually comprise at least two cryptic lineages that coexist sympatrically across most of Mexico. To broaden our understanding of the dynamics of diversification to other rubyspot lineages, we investigated the phylogeographic history of smoky rubyspot damselflies (Hetaerina titia) using genomic data collected across Central and North America. Unexpectedly, we found evidence of reproductive isolation between the highly genetically differentiated Pacific and Atlantic lineages of H. titia in a narrow secondary contact zone on the Isthmus of Tehuantepec, Mexico. We then fit models of historical demography to both H. americana sensu lato and H. titia to place these comparisons in a temporal context. Our findings indicate that Pacific and Atlantic lineages of H. titia split more recently than the broadly sympatric lineages within H. americana sensu lato, supporting key assumptions of the non-ecological speciation model and demonstrating that these two pairs of sister lineages are at different stages of the speciation cycle.

Keywords: Odonata; Zygoptera; ddRAD; population genomics; speciation.

FIGURE 1

Ancestry estimates for 205 Hetaerina titia with a dataset of 3819 unlinked biallelic autosomal SNPs. SNPs were generated by mapping ddRAD reads to the draft genome of H. titia . LEA was run for 20 repetitions and an alpha value of 100. (a) The mean estimate of ancestry proportion for all samples within each sample site of H. titia across Central and North America, (b) Isthmus of Tehuantepec and Belize and (c) Costa Rica. Within panels a, b and c, the area of each pie chart is proportional to the number of samples from each site and then coloured by the mean proportion of estimated ancestry (either South Atlantic, North Atlantic or Pacific) across all samples from each site. (d) Estimate of ancestry for each individual. Samples are ordered by drainage, then country and then latitude. Rivers and drainage basins from Hydrosheds. Topography data from the R package elevatr. The black boxes shown in panel (a) are the bounding areas for panels (b) and (c). The black box in panel (b) is the bounding box for Figure 6b and the five samples from the site with an identified hybrid individual are underlined in panel (d).

FIGURE 2

Principal component analysis of 205 Hetaerina titia with a dataset of 3819 unlinked biallelic autosomal SNPs that were generated by mapping ddRAD reads to the draft genome of H. titia . Percentages indicate how much variation is explained by each component and colour indicates the highest assigned ancestry population from sNMF for each individual. The single point directly between the main Pacific and Atlantic cluster is the putative F₁ hybrid. A PCA plot for H. americana sensu lato showing broadly the same level of differentiation between samples, without any individuals showing introgression, is included in Figure S10.

FIGURE 3

The maximum likelihood tree for (a) Hetaerina titia and (b) H. americana sensu lato calculated using RAxML with 3020 SNPs for H. titia and 3949 SNPs for H. americana and H. calverti and mapped onto the genome of H. titia. Scale bar indicates the mean number of substitutions per SNP site. Due to exclusion of invariant sites and differences in the total number of SNPs used in each analysis, scale bars should not be used to compare phylogenetic distances within H. titia to distance within H. americana sensu lato. The nodes marked with a blue star ‘*’ indicate a bootstrap support value (out of 100) of < 95%. The tree tips are coloured according to the species and the max sNMF ancestry assignment (K = 3). The geographical location of each sample is shown in the bottom two maps. Each pie chart shows the number of samples assigned to each ancestry cluster from each sample site, split between H. titia and H. americana sensu lato ( H. americana /calverti).

FIGURE 4

(a) Estimates of divergence dates (million years ago—mya) between populations of Pacific Hetaerina titia (titia‐Pac), Atlantic H. titia (titia‐NAtl and titia‐SAtl), H. americana (americana‐N and americana‐S) and H. calverti (calverti) calculated using SNAPP analysis in Beast. Node labels indicate the mean estimated divergence date with 95% highest posterior density in blue. All branches had a posterior distribution of 1. Tree plotted in R using the packages treeio and ggtree. Input data was 552 autosomal SNPs called using a de novo method of SNP calling (ipyrad). (b, c, d) The prior and posterior distribution (where applicable) of divergence times between the major lineages. The three different histograms denote the posterior distribution of the divergence times using three different SNP datasets, those mapped the draft genome of H. americana (dark grey), mapped the draft genome of H. titia (grey) and de novo SNP calling. The prior, where applicable, is denoted by the grey density distribution. (e) Comparison between the divergence times of H. americana and H. calverti and Atlantic and Pacific H. titia . Each point is the divergence times from a tree in the posterior distribution, the black line indicates values where the divergence times between the lineages are equal.

FIGURE 5

The estimated divergence times (Mya = million years ago) and effective population size (theta—Ne in millions of individuals) from G‐PhoCS analysis of Hetaerina titia and H. americana. Migration rate is the number of individuals per generation with vertical arrows indicating direction of migration (from and to). All models ran for 1,000,000 iterations with 10% burn in. Blue bars show 95% highest posterior density for each divergence date (a) Model estimates for H. titia with no migration bands. The estimated divergence time for Atlantic and Pacific H. titia was 2.72 mya (2.65–2.80 mya HPD) and divergence time between Northern and Southern Atlantic clusters was estimated as 0.59 mya (0.55–0.62 mya HPD). (b) Model estimates for H. americana and H. calverti with no migration bands. The divergence time for H. calverti and H. americana was estimated to be 4.7 mya (3.47–4.68 mya HPD). The Northern and Southern H. americana clusters diverged 1.60 mya (1.55–1.66 mya). (c) Model estimates for H. titia demography with migration bands between Northern and Southern Atlantic H. titia. The divergence time between Northern and Southern Atlantic H. titia was 0.87 mya (0.76–1.01 mya HPD) and the split between Atlantic and Pacific H. titia was 2.71 mya (2.64–2.79 HPD). (d) Model estimates for H. americana and H. calverti with migration bands between North and Southern H. americana . The divergence time for H. calverti and H. americana was estimated to be 4.68 mya (4.57–4.79 mya HPD). The Northern and Southern H. americana clusters diverged 1.75 mya (1.68–1.83 mya). G‐PhoCS runs presented here are conducted on the RAD loci mapped de novo using ipyrad.

FIGURE 6

Hybrid zone between Pacific and Southern Atlantic Hetaerina titia in the Isthmus of Tehuantepec. (a) Genotypes for 914 autosomal SNPs and 19 sex‐linked SNPs that had a < 0.8 allele frequency difference between Pacific and Atlantic individuals (calculations excluded samples from CUAJ01/02). Each sample is positioned along the y‐axis with each SNP ordered by the position along each chromosome along the x‐axis. The F₁ hybrid is sample CUAJa02. Each SNP is coloured by whether they were homozygous for the Pacific allele (0/0—red), homozygous for the Atlantic allele (1/1—beige) or heterozygous (0/1—orange). (b) Sample locations around the Isthmus of Tehuantepec. The Atlantic and Pacific watershed boundary is shown in black (see Figure 1b for a map of wider region and Figure S23 for a map of terrain height rather than a shaded relief). (c) A triangle plot showing the hybrid index, measuring the percentage of ‘parental’ genotype and the heterozygosity of each sample. A theoretical F₁ hybrid would be placed at the top corner of the triangle. SNPs on the X chromosome were excluded when calculating the hybrid index and heterozygosity.