Extraordinarily rapid speciation in a marine fish

Abstract

Divergent selection may initiate ecological speciation extremely rapidly. How often and at what pace ecological speciation proceeds to yield strong reproductive isolation is more uncertain. Here, we document a case of extraordinarily rapid speciation associated with ecological selection in the postglacial Baltic Sea. European flounders (Platichthys flesus) in the Baltic exhibit two contrasting reproductive behaviors: pelagic and demersal spawning. Demersal spawning enables flounders to thrive in the low salinity of the Northern Baltic, where eggs cannot achieve neutral buoyancy. We show that demersal and pelagic flounders are a species pair arising from a recent event of speciation. Despite having a parapatric distribution with extensive overlap, the two species are reciprocally monophyletic and show strongly bimodal genotypic clustering and no evidence of contemporary migration, suggesting strong reproductive isolation. Divergence across the genome is weak but shows strong signatures of selection, a pattern suggestive of a recent ecological speciation event. We propose that spawning behavior in Baltic flounders is the trait under ecologically based selection causing reproductive isolation, directly implicating a process of ecological speciation. We evaluated different possible evolutionary scenarios under the approximate Bayesian computation framework and estimate that the speciation process started in allopatry ∼2,400 generations ago, following the colonization of the Baltic by the demersal lineage. This is faster than most known cases of ecological speciation and represents the most rapid event of speciation ever reported for any marine vertebrate.

Keywords: Baltic Sea; ecological speciation; evolution; genomics; rapid speciation.

Fig. 1.

Map of the sampling locations showing modeled mean bottom salinity in the Baltic Sea from Bendtsen et al. (62). Areas theoretically suitable for pelagic spawning are colored in green, yellow, and red. Locations are color-coded based on expected spawning behaviors. ÅLA, Åland Archipelago; BAR, Barsebäck; BOR, Bornholm Basin; DAB, Dabki; GDY, Gdynia; GUL, Gullmaren; HAN, Hanko (Western Gulf of Finland); HEL, Helsinki (Gulf of Finland); IRB, Irbe Strait; KAT, Kattegat; KVÄ, Kvädöfjärden; LAT, Latvian Sea; and ÖLA, Öland Island. More than 90% (Table S1) of samples from BOR (pelagic habitat) and ÖLA (demersal habitat) were ready for spawning at the time of capture, hence these locations were used as reference pelagic and demersal populations, respectively.

Fig. 2.

(A) PCA of allele frequency data from all loci. Codes are as follows: DEM, putative demersal locations; NS, North Sea samples; PEL, putative pelagic locations within the Baltic; TZ, transition zone between the Baltic Sea and the North Sea (Barsebäck samples). (B) Relationship between geographic and genetic distance among sampling locations. “Within” refers to pairwise comparisons among locations belonging to the three major genetic clusters identified by PCA, fastSTRUCTURE, and K-means clustering. (C) Individual ancestries from fastSTRUCTURE analysis (K = 3).

Fig. S1.

PCA and distribution of Baltic Sea individuals along the first PC axis based on allele frequency data from loci from the core F_ST distribution estimated by OutFLANK (A and B) and from a dataset excluding all loci identified as outliers by any of the outlier tests performed (C and D). E and F represent the relationship between geographic and genetic distances among sampling locations from the same two “neutral” datasets. “Within” refers to pairwise comparisons among locations belonging to the three major genetic clusters identified by PCA, fastSTRUCTURE, and K-means clustering.

Fig. S2.

Results from genotypic clustering analyses. A and B represent the Akaike information criterion (AIC) and Bayesian information criterion (BIC), respectively, for runs of K-means clustering using increasing numbers of K, and suggest K = 3 and K = 2, respectively, reflecting higher penalties of BIC for increasing model complexity. C and D represent prediction error from fivefold cross-validation for the fastSTRUCTURE analyses and plots of the log-marginal likelihood lower bound (LLBO) at increasing number of K, using both the simple and logistic models. E shows the assignment probability and estimated ancestry from DAPC and fastSTRUCTURE analyses when K = 2 and K = 3, respectively.

Fig. 3.

Frequencies of outlier loci in each sampling location. SNPs are identified by the locus number (above) and the SNP position within the locus in base pairs (below). The first eight loci are outliers between pelagic and demersal flounders, and the last three are outliers between North Sea and Baltic pelagic locations. Individuals in Irbe were grouped separately in pelagic (above) and demersal (below) groups based on PCA, fastSTRUCTURE, and K-mean clustering results. Blue alleles are characteristic of pelagic Baltic populations; red alleles are characteristic of demersal Baltic populations. Black alleles are characteristic of North Sea populations.

Fig. S3.

BAYESCAN analyses for outliers between pelagic and demersal flounders (A) and between North Sea and Baltic Sea pelagic locations (B). Outliers that were also identified by all other outlier tests (Fdist2, OutFLANK, and FLK) are represented by filled red circles. C–H show observed frequencies (gray histograms) of F_ST and expected distribution of F_ST under neutral expectations (black lines) as estimated by OutFLANK. (C) F_ST distribution between pelagic and demersal Baltic flounders; (D) right tail of the distribution, showing the eight outliers that were jointly identified by all tests. (E) F_ST distribution when comparing pelagic locations in the Baltic with North Sea; (F) right tail of the distribution, showing the three outliers that were jointly identified by all tests. (G) Distribution of global F_ST between pelagic flounders locations within the Baltic Sea and (H) comparisons among locations of demersal flounders within the Baltic Sea grouped according to sampling location. No statistically significant deviations from expected distributions were found among different sampling locations within the two Baltic populations. BAYESCAN analyses also failed to reveal any outliers in the same comparisons.

Fig. 4.

Unrooted ML tree based on the concatenated sequences from 2,051 loci for a subset of 57 individuals. Black filled circles represent putative pelagic individuals. Red filled circles represent putative demersal individuals. Yellow circles represent individuals from a putative demersal location (Irbe), which clustered in all previous analyses with pelagic individuals. Branch support values represent ML bootstrap support and Bayesian clade credibility values, respectively. Only the single node that shows high support (>80) in both ML and Bayesian analyses is labeled.

Fig. 5.

Scenarios tested via ABC-RF. Scenario 1: double-colonization scenario with early invasion by demersal flounders. Scenario 2: single colonization scenario with early invasion of pelagic flounders. Scenario 3: single-colonization scenario with early invasion of demersal flounders. Scenario 4: double-colonization scenario with early invasion by pelagic flounders. Each lineage is color coded: black lines represent the North Sea lineage (NS), blue lines represent the Baltic pelagic lineage (PEL), and red lines represent the Baltic demersal lineage (DEM). Broken lines represent possible demographic changes in the ancestral population (black broken line) as well as potential bottlenecks at the time of invasion (blue and red broken lines). The blue dotted line with arrowhead represents secondary introgression among the two Baltic species. All divergence scenarios were tested against each other under an orthogonal combination of all other demographic parameters (changes in N_e and introgression). Full details of all models tested via ABC-RF are given in Supporting Information. The most likely model was scenario 1, including a demographic expansion in the ancestral population and a mild bottleneck in the demersal lineage.

Fig. S4.

PCA in the space of summary statistics, showing datasets simulated from the prior distribution of the parameters (green open circles), from the posterior predictive distribution (green filled circles), as well as the observed dataset (yellow circle). A–D are plots of the PCA using different combinations of the first five principal components, which cumulatively represent 89% of the total variation. The observed dataset is always located within a small cluster from the posterior predictive distribution.