Stable species boundaries despite ten million years of hybridization in tropical eels

Abstract

Genomic evidence is increasingly underpinning that hybridization between taxa is commonplace, challenging our views on the mechanisms that maintain their boundaries. Here, we focus on seven catadromous eel species (genus Anguilla) and use genome-wide sequence data from more than 450 individuals sampled across the tropical Indo-Pacific, morphological information, and three newly assembled draft genomes to compare contemporary patterns of hybridization with signatures of past introgression across a time-calibrated phylogeny. We show that the seven species have remained distinct for up to 10 million years and find that the current frequencies of hybridization across species pairs contrast with genomic signatures of past introgression. Based on near-complete asymmetry in the directionality of hybridization and decreasing frequencies of later-generation hybrids, we suggest cytonuclear incompatibilities, hybrid breakdown, and purifying selection as mechanisms that can support species cohesion even when hybridization has been pervasive throughout the evolutionary history of clades.

Fig. 1. Genomic and morphological variation in tropical eels.

a Distribution of Anguilla species in the Indo-Pacific. The color and position of dots within hexagons indicate species presence within the region covered by the hexagon, according to the Global Biodiversity Information Facility database and our own collection. Sampling locations are indicated with black dots. Numbers following location names specify the number of samples taken. Stacked bars indicate the species identities of individuals, according to mitochondrial and morphological species assignment. b Morphological variation among the four species A. marmorata (n = 100), A. megastoma (n = 30), A. obscura (n = 30), and A. interioris (n = 1). Dots represent individuals and are colored according to mitochondrial species identity. c Genomic principal component analysis (PCA) based on 155,896 variable sites. Specimen IDs are given for individuals with intermediate genotypes. The cyan circle indicates a cluster of 11 individuals mitochondrially assigned to A. marmorata (SAA16011, SAA16012, SAA16013, SAA16027, SAW17B27, SAW17B49, VAG12012, VAG12018, VAG12019, VAG13071, VAG13078), in addition to the highlighted VAG12044 that is mitochondrially assigned to A. megastoma. d Time-calibrated phylogeny based on 5000 transition sites. Each individual tree shown in gray represents a sample from the posterior tree distribution; a maximum-clade-credibility summary tree is shown in black. Color code in b, c, and d is identical to a. PC, principal component; AD, distance between the dorsal fin and the anus; PDH, predorsal length without head length; TL, total length.

Fig. 2. Contemporary hybridization among tropical eels.

a–d Genomic variation inferred in a PCA of all individuals of A. marmorata, A. megastoma, A. obscura, and A. interioris (Supplementary Fig. 5c), shown separately for four hybridizing species pairs: A. marmorata and A. megastoma (a), A. marmorata and A. obscura (b), A. megastoma and A. obscura (c), and A. marmorata and A. interioris (d). Individuals with intermediate positions are marked in gray with specimen IDs; species color code is as in Fig. 1. e Ancestry painting for 20 hybrids between A. marmorata and A. megastoma. The top and bottom horizontal bars represent 302 sites that are fixed for different alleles between the two species; all other bars indicate the alleles at each of those sites. White color indicates missing data. Heterozygous alleles are shown with the top half in each bar matching the second parental species and vice versa. f Ancestry painting for three contemporary hybrids between A. marmorata and A. obscura, based on 742 sites fixed between these two species. g Ancestry painting for one hybrid between A. megastoma and A. obscura, based on 525 fixed sites. h Ancestry painting for one hybrid between A. marmorata and A. interioris, based on 429 fixed sites. i Histogram of heterozygosity observed in hybrids. j Histogram of the proportions of hybrid genomes derived from the maternal species (according to mitochondrial sequence data). k Histogram of the relative morphological similarities between hybrids and the maternal species, measured as the relative proximity to the mean maternal phenotypes, compared to the proximity to the mean paternal phenotype. l Comparison of the proportions of hybrids’ genomes derived from the maternal species and the similarity to the mean maternal species’ phenotype. The dotted line indicates a significant positive correlation between the two measures (two-tailed t test; t = 3.1, p = 0.008, R² = 0.381). PC, principal component; mito., mitochondrial genome; AD, distance between the dorsal fin and the anus; PDH, predorsal length without head; TL, total length.

Fig. 3. Past introgression among tropical eels.

a Likelihood support of individual RAD loci for different relationships of A. interioris: as sister to A. marmorata and A. luzonensis (bottom left), as sister to A. obscura and A. bicolor (bottom right), and as sister to a clade formed by those four species (top). The position of each dot shows the relative likelihood support of one RAD locus for each of the three tested relationships, with a distance corresponding to a log-likelihood difference of 10 indicated by the scale bar. The central triangle connects the mean relative likelihood support for each relationship. A black dot inside that triangle marks the central position corresponding to equal support for all three relationships. Sample sizes (n) report the number of loci that support each of the two competing relationships connected by that edge. b Heatmap indicating maximum pairwise D (above diagonal) and f₄ (below diagonal) statistics (see Table 1). Combinations marked with “x” symbols indicate sister taxa; introgression between these could not be assessed. Asterisks indicate the significance of f₄ values (*p < 0.05; **p < 0.01; ***p < 0.001; not adjusted for multiple comparisons; see Table 1 for precise values), determined through one-sided comparison with coalescent simulations with the F4 software. The cladogram on the left summarizes the species-tree topology according to a and the significant signals of introgression according to b. c, d Comparisons of the maximum D value per species with the species' geographic range or population mutation rate Θ. Geographic range was measured as the number of geographic hexagons (see Fig. 1) in which the species is present, and Watterson’s estimator was used for the population mutation rate Θ. n.s. not significant. e Genomic patterns of phylogenetic relationships among A. marmorata, A. obscura, and A. megastoma, based on WGS reads mapped to the 11 largest scaffolds (those longer than 5 Mbp) of the A. anguilla reference genome. Blocks in light gray show 20,000-bp regions (incremented by 10,000 bp) in which A. marmorata and A. obscura appear as sister species, in agreement with the inferred species tree; in other blocks, A. megastoma appears closer to either A. obscura (gray) or A. marmorata (dark gray).

Fig. 4. Potential causes of cytonuclear incompatibility between A. marmorata and A. megastoma.

a Dot plot visualizing an inversion between A. megastoma scaffold scf7180010919884 (purple) and the homologous A. marmorata scaffold scf7180010922493 (cyan). Each dot represents a tuple of nine nucleotides that are identical or the reverse complement between the two scaffolds. Light gray rectangles mark regions homologous to exons of the zebrafish (Danio rerio) myhc4 gene; order and orientation of these exons are indicated at the top (Supplementary Table 12). The inversion may thus affect the transcript of this gene in A. megastoma. Other detected genomic rearrangements are shown in Supplementary Figs. 21 and 22. b Illustration of fixed sites in coding sequences between A. marmorata and A. megastoma. Light gray rectangles mark regions homologous to exons of the zebrafish ttna gene (Supplementary Table 14). Dark gray lines indicate RAD sequences for 72 and 26 “core” individuals of A. marmorata and A. megastoma, respectively. Nonreference genotypes are colored according to the species in which the genotype is most frequently found. Within exon 191, an adenine to guanosine substitution in A. marmorata at position 32,226 of A. anguilla scaffold scf1929 changes amino acid 5732 of the protein from aspartic acid to glycine. Other protein-coding sites fixed between A. marmorata and A. megastoma are listed in Supplementary Table 13. c Nonsynonymous substitutions between the A. marmorata (inner, cyan circle) and A. megastoma (outer, purple circle) mitochondrial genomes. Gray segments indicate protein-coding genes and nonsynonymous substitutions within these are shown with black lines. All mitochondrial amino-acid changes are listed in Supplementary Table 15.