Rapid diversification underlying the global dominance of a cosmopolitan phytoplankton

Abstract

Marine phytoplankton play important roles in the global ecosystem, with a limited number of cosmopolitan keystone species driving their biomass. Recent studies have revealed that many of these phytoplankton are complexes composed of sibling species, but little is known about the evolutionary processes underlying their formation. Gephyrocapsa huxleyi, a widely distributed and abundant unicellular marine planktonic algae, produces calcified scales (coccoliths), thereby significantly affects global biogeochemical cycles via sequestration of inorganic carbon. This species is composed of morphotypes defined by differing degrees of coccolith calcification, the evolutionary ecology of which remains unclear. Here, we report an integrated morphological, ecological and genomic survey across globally distributed G. huxleyi strains to reconstruct evolutionary relationships between morphotypes in relation to their habitats. While G. huxleyi has been considered a single cosmopolitan species, our analyses demonstrate that it has evolved to comprise at least three distinct species, which led us to formally revise the taxonomy of the G. huxleyi complex. Moreover, the first speciation event occurred before the onset of the last interglacial period (~140 ka), while the second followed during this interglacial. Then, further rapid diversifications occurred during the most recent ice-sheet expansion of the last glacial period and established morphotypes as dominant populations across environmental clines. These results suggest that glacial-cycle dynamics contributed to the isolation of ocean basins and the segregations of oceans fronts as extrinsic drivers of micro-evolutionary radiations in extant marine phytoplankton.

Fig. 1. Relationship between genetic structure and morphotypes in G. huxleyi.

a Principal component analysis (PCA) based on 2,086,643 SNPs recovered from 47 G. huxleyi genomes; b Relationship between coalescent species phylogeny (ASTRAL tree based on 1000 supergenes) and DAPC clustering; c Correspondence between morphotypes and lineages within G. huxleyi, and sub-lineages within A1 (scale bar = 4 μm). Variable elements in relation to genotypes are highlighted in the schematics under the SEM pictures; d Distribution of coccolith length for 5 randomly chosen strains representing each clade and sub-clade, with a jittered box-plot on the left and a half-violin plot on the right for each group; e Matrix plot of Bonferroni corrected p-value corresponding to the Dunn-test for the comparison of coccolith length measurements between groups.

Fig. 2. Excess of allele sharing and differentiation in G. huxleyi.

a f-branch (f_b) statistics between lineages and sub-lineages. The gradient represents the f_b score, grey blocks represents tests not consistent with the species tree (for each branch on the topology of the y axis, having itself or a sister taxon as donor on the topology of the x axis); asterisks denote block jack-knifing significance at p < 0.05 (after Bonferroni correction); b Phylogenetic network inferred using a subset of 83,563 SNPs across 47 strains; c Combined box and violin plots showing the distribution of genetic differentiation per sites between lineages for synonymous sites.

Fig. 3. Tempo of diversification between and within lineages in G. huxleyi.

a Phylogenetic chronogram of G. huxleyi based on analysis of genome sequence data of 63 strains. The phylogeny was rooted using one strain of G. muellerae and three strains of G. ericsonii/parvula as an outgroup. Every node in the phylogeny has a posterior probability at 1. The dating of ancestral nodes is based on relaxed molecular clock calibrated with the first appearance of G. huxleyi (node 1; 290 ka) in the fossil record. Details for nodes are reported in Supplementary Table S7. 95% Highest posterior density intervals for ages are shown as grey bars. b Visual representation of the parameters inferred for consecutive speciation events between lineages A1, A2 and B. Circles reflects effective population size (Ne) estimated for extant and ancestral species, on nodes and on leaves. Node 2 and 4 are highlighted in bold for correspondence with the chronogram. Divergence time are provided as an interval accounting for μ (=mutation rate) uncertainty [86]. Arrows represent secondary contact events with migration values in italic. All parameter estimates are listed in Supplementary Table 9. g: generation time. c Absolute abundance of G. huxleyi in sites U1475 [87] and SO139-74KL [88]. d Global Δ Sea Surface Temperature [89] (ΔSST; blue line) and LR04 [90] (red line) over 350 ka; (MIS: Marine Isotopic Stage).

Fig. 4. Relationship between genetic lineages and environmental variables.

a Distribution map of lineages and sub-lineages based on strains used in this study. b Redundancy Analysis plot with constrained predicators. c Strain distributions lineages/ and sub-lineages per relevant environmental parameters (significances (p values) of the Dunn test correspond to: *p < 0.05; **p < 0.01, ***p < 0.001).