Rapid diversification underlying the global dominance of a cosmopolitan phytoplankton

Genetic and morphological delineation between G. huxleyi strains

We first assessed genetic variability through analysis of genomic polymorphism to determine whether distinct genetic lineages exist in G. huxleyi and to test whether these relate to morphotypes. We used 2,086,643 high-quality biallelic single nucleotide polymorphisms (SNPs) retrieved from the 47 clonal culture strains with the best genome sequence coverage (>20×). A principal component analysis (PCA) and a discriminant analysis in principal component (DAPC) both delineate three well-defined genetic groups, with the distribution of strains being unequal and with no overlap on the principal components (Fig. 1a; Supplementary Fig. S3a,b). With regards to population structure, the DAPC analysis suggested that 3 clusters (K = 3) can be used to depict a genotype membership matrix for each strain (Fig. 1b; Supplementary Fig. S4). As such, it confirmed the three-lineage delineation proposed by the PCA, while illustrating no admixture between lineages.

Phylogenetic inference based on alignments with higher mapping coverage only (47 strains) or including sequences with lower mapping coverage (59 strains) all supported segregation of strains into three main lineages, which we term clades A1, A2 and B, with A1 and A2 being more closely related to each other than to B (Fig. 1b; Supplementary Fig. S5a, b). This delineation is congruent with previous studies on the phylogeny of the Gephyrocapsa genus [17, 46, 65]. These clades also correspond to differences in morphotypes (Fig. 1b, c). All strains in clade A1 produce unambiguous A-group coccolith morphotypes (type A and type R). Similarly, all strains in clade B produce unambiguous B-group coccolith morphotypes (type B and type O). Clade A2 is less distinctive, with strains producing lightly calcified type A coccoliths. Some of these strains could be classified as type B/C [66] or C (both regarded as B-group morphotypes), but distinctive by the lower elevation of distal shield elements and by greater degree of calcification of the central area grid (which is reduced and sometimes absent in morphotypes B/C and C). At a finer level, clade A1 is composed of four sub-clades, which we term A1a, A1b, A1c, and A1d. Strains in sub-clades A1a, A1c and A1d all produce coccoliths with type A morphologies and distinctive degrees of calcification: strains in the sub-clade A1a form relatively lightly calcified coccoliths with regular elements, while strains in sub-clades A1c and A1d produce similar moderately calcified coccoliths, sometimes with conspicuous irregularities (inner tube elements overlapping into the central area). Strains in clade A1b produce distinct coccoliths exhibiting A-group morphology but with heavy calcification, including forms with heavily calcified shields which have been termed type R and also forms with heavily calcified central areas which have been referred to as “type A overcalcified”. Some clade A2 strains produce coccoliths with a similar morphology to strains in A1a, indicative of partially cryptic lineages (Supplementary Fig. S2; Supplementary Table S4).

The congruence between morphotypes and clades is also supported by significant differences in the length of coccoliths measured between some of the clades (Fig. 1d, e). The morphogroups A and B differ significantly, and insignificant comparison relates to the comparison of sub-clades against the clade A2, which reinforces the closest relationship between A1 and A2. We denote also that the case of A1a and A2 demonstrating no significant difference in coccolith length concurs with the cryptic delineation mentioned above.

Based on the clustering analyses and the phylogenetic reconstructions, we tested whether different groupings are distinct species with regards to the null hypothesis “G. huxleyi is a single species”, which correspond to the current state of taxonomy. Species delimitation based on comparison of Marginal Likelihood Estimators (MLE) with Bayes Factors (BF) supported the hypothesis that the three lineages depicted by ordination and phylogenetic reconstructions are distinct species as the best model (Table 1).

D-statistics calculated to estimate gene flow reveal a non-significant excess of alleles shared between the three lineages (Fig. 2a; Supplementary Table S5). Fbranch statistics, (f_b) revealed significant signatures of gene-flow between sub-lineages within A1 associated with correlated estimates in relation to A1a, A2 and B (Fig. 2a) [60]. Signatures on the basal branch of diversification in A1 may correspond to genetic exchanges between A1 and B, with gene-flow signatures attributed to A2 corresponding to correlated estimates due to common ancestry. Recent signatures of gene-flow throughout the evolution of A1 are thus likely associated to the common ancestry between A1a, A2 and B during gene-flow events between the sub-lineages, as supported by the non-significant D statistics between the three lineages. Moreover, the phylogenetic network revealed similar convolutions between A1 sub-lineages but clear separation of the main lineages and longer branches in the A2 lineage (Fig. 2b).

Comparison of pairwise differentiation estimates per sites (F_ST), for synonymous sites, supported A1 and A2 to be more closely related (F_ST = 0.153) than A1 and B (F_ST = 0.191) as depicted by the phylogenies (trees and network). Although, the differentiation between A2 and B was more marked than any other comparison (F_ST = 0.276; Fig. 2c; Supplementary Tables S6-7) which could relate to the long branches observed in the phylogenetic network (Fig. 2b). These results rather suggest a rapid divergence within the clade A2, likely since its emergence.

These analyses of genetic divergence and phylogenetic relationship therefore indicate that G. huxleyi, usually considered a single species, has differentiated into at least three reproductively isolated species, A1, A2 and B. These results expand on recent phylogenomic and population genetic results [17, 25], and support the binary morphogroup classification [28]. Within these lineages, morphotypes are structured into distinct clades, providing clear evidence that G. huxleyi morphotypes correspond to distinct genotypes.

Environmental drivers of diversification in the three species of Gephyrocapsa

We established a timeline of diversification by comparing the fossil record attributed to G. huxleyi with genomic divergence time based on a molecular clock reconstruction and joint-site frequency spectrum (JSFS) modelling (Fig. 3; Supplementary Table S7, Supplementary Fig. S5). According to sedimentary records, first occurrences (FO) of G. huxleyi-like ancestors occurred synchronously across low latitude sites during glacial stage MIS8 around 290 ka [12]. We used this date to calibrate our chronogram (Fig. 3a; Supplementary Table S7).

Our divergence time reconstruction indicates that modern G. huxleyi populations originated from the divergence of A and B clades in the MIS6 glacial period preceding the end of the Pleistocene (152 (124–185) ka in our study and ~170 (88–223) ka in [17]), which corroborates the time range deduced from sedimentary observations for the appearance of B-like morphotypes [67, 68]. Our JSFS analysis supported this divergence, but as an episode of geographic separation (vicariance) and niche partitioning between A1 and B (138 (126–152) ka; Fig. 3b). Sharper isolation of fronts and stronger physical structuring during this glacial period could have accounted for this divergence in isolating distinct populations. In contemporary oceans, the partially cryptic A1a and A2 clades have sympatric distributions associated with low latitude and warm water masses. Sub-clades A1c and A1d together with clade B form another sympatric assemblage, but associated with temperate and subpolar waters. There is limited overlap between the distributions of these two groupings, which therefore form allopatric assemblages defined by their ecologies. Although, the distribution of strains belonging to sub-clade A1b does overlap with those of A1a, A2 and B (Fig. 4a). Accordingly, and given the basal divergence of A1a in A1, it is likely that A1, A2 and B lineages have been diverging with distinct preferences to sea surface temperature, nutrients concentrations and carbonate chemistry, as suggested by the redundancy analysis (RDA; Fig. 4b, c).

Further diversifications followed during the MIS5 interglacial period, with a significant discrepancy between the chronogram and JSFS-modelled divergence. While the chronogram suggests diversification occurred within B (115 (78–170) ka and 92 (68–125) ka) and between A1 and A2 (103 (88–119) ka), the genetic model indicates that A1 and A2 diverged around the same time as A2 and B (respectively 102 (93–112) and 98 (89–107) ka; Fig. 3b; Supplementary Table S7). This result reinforces the view that the earliest divergence was between A1 and B and indicates a convoluted divergence for A2 (cf. fbranch and phylogenetic network; Fig. 2a, b), which could be the result of hybridisation during early stage of speciation. It confirms gene flow signatures associated with potential interaction between A1 and B prior to diversification within a branch of A1 (Fig. 2a). The three divergences within G. huxleyi tested with JSFS retrieved a similar mode of speciation as for the more general Gephyrocapsa group, that has undergone speciation events followed by occasional gene flow during secondary contacts [25]. In G. huxleyi, putative secondary contacts occurred recently during MIS1 with an extremely reduced measure of gene-flow (M < 1; 10⁻⁶ < M < 10⁻⁵; M: Migration, number of individuals from a population that exchanged gene with another), confirming non-significant D-statistics found between the lineages. Overall, these major divergence events were not restricted to a particular environmental scenario, but rather to the fluctuating state of MIS5 in terms of niche expansion and compression.

The first divergences within B and A2 are associated with sub-clades composed of strains originating from distinct oceans, suggesting patterns of geographic isolation, as a result of migrating fronts in relation with gradual ice growth that followed the interglacial maximum. Distinct morphotypes are associated with geographic separations within B, morphotype B being found in the Atlantic Ocean and morphotype O in the Pacific Ocean, while morphological variations are less consistent in A2 (Supplementary Figs. S1–2 and Supplementary Tables S1-2 and S4). Within A1, intra-specific divergences then occurred through different events of vicariance during the MIS4-2 glacial period, establishing genetically distinct populations along a latitudinal gradient. A significant eccentricity minimum associated with these events could account for interactions between newly diversified populations within A1 due to stronger compression of ecological niches [69]. In this scheme, gene-flow between emerging populations could have played a role in the adaptive process associated with expanded habitat to higher latitudes. Gene flow events within A1 redistributed allelic composition associated with common ancestors of the three lineages into newly formed populations, contributing to potential adaptation to environmental variability along latitudinal gradients. The 400 ky cycle of the absolute eccentricity minimum, which had not occurred since the origin of the G. huxleyi lineage, may be related to major events of diversification within Gephyrocapsa [69].

Early diversifications within A1 fit with the timing of the diachronic acme (≥50% dominance in the total fossil coccolithophore flora) along a latitudinal gradient (85 ka in low latitudes, 73 ka in transitional latitudes, 61 ka in high latitudes of the North Atlantic Ocean) [70, 71], and were also associated with a notable increase in coccolith size (to >4 µm) until the last glacial maximum (LGM) [72,73,74]. Physiological experiments have demonstrated that G. huxleyi is extremely competitive in certain nutrient limitation scenarios [75, 76]. Therefore, G. huxleyi might have benefitted from further fertilisation of the ocean linked to reduced sea level during this glacial period [77], which also contributed to increase seawater alkalinity levels, favouring calcification [78] (as may be now the case in the Black Sea [79], for example). Such bloom-forming conditions would account for modulation of the life cycle toward clonal reproduction [28], leading to increased genetic diversity through fixation of heterozygous substitutions [80]. In this context, distinct populations, perhaps even strains/individuals, may reach reproductive isolation faster, in association with increased probabilities for incompatibilities in offspring during secondary contact [47, 81]. This micro-evolutionary pattern integrates well with previously described patterns of speciation [17, 25], in which macro-evolutionary size variations observed in the fossil record were caused by repeated species radiations rather than fluctuations in the relative abundance of large- and small-celled species. By contrast, interglacials, and especially the current MIS1, may witness the selective impacts of an increase in atmospheric CO₂ [74] (i.e., global warming), as attested by the gradual reduction of coccolith size and abundance in the fossil record [72,73,74].

Taxonomic implications

Based on the results presented herein, we propose a formal taxonomic reassessment that integrates morphological, phylogenetic, admixture, and ecological information relative to the genus as diagnostic features. The species G. huxleyi (=clade A1) originated relatively recently along with two other species, leading us to split the entity G. huxleyi into three species by emending the former G. huxleyi, erecting G. pseudohuxleyi sp. nov. (=clade A2), and reinstating G. pujosae comb & stat. nov. (=clade B). We believe that this new nomenclature will be useful for future studies of assemblages using (meta)genomic comparison, coupled or not with electron microscopy, taking into account that some populations are (pseudo)cryptic and others may not calcify. This proposal reflects current knowledge of the process of speciation and may change with future evolution of concepts, in the same manner as taxonomic considerations for this complex have evolved over the 20^th century [82]. For practical reasons, particularly for studies that do not employ genomic or electron microscopy analyses (as is currently the case for example for most micropalaeontological investigations), these three species can be accommodated into a “superspecies” concept under the name G. huxleyi, as has implicitly been the case for years under the name Emiliania huxleyi.

Gephyrocapsa huxleyi (Lohmann) Reinhardt emend. Bendif, Probert, Beaufort, Rickaby & Archontikis

Description

Coccoliths with moderately elevated distal shield (2–4 μm length) and elements of variable width (0.05–0.25 μm); inner tube with variable width, sometimes irregular, sometimes irregularly extended on the central area; central area sometimes with a grill of curved rods, sometimes thick lath-lick element forming a solid plate with irregular holes, sometimes strainer-like grill with regular holes, sometimes closed. Comprise previously described morphotypes A, over-calcified and R.

Genetic diagnosis

genetically distant from other species of Gephyrocapsa by genome sequences. Admixture pattern distinct from G. pujosae comb. nov. and G. pseudohuxleyi sp. nov. and forms the phylogenetic clade A1 (Fig. 1).

Basionym

Pontosphaera huxleyi Lohmann 1902 p. 130, pl. 4 Figs. 1–6, pl. 6 Fig. 69 [19].

Synonyms

Hymenomonas huxleyi (Lohmann) Kamptner [83]; Coccolithus huxleyi (Lohmann) Kamptner [84]; Emiliania huxleyi (Lohmann) Hay & Mohler [26]; E. huxleyi var huxleyi Medlin & Green [30].

Lectotype

Lohmann 1902 p. 130, pl. 4 Figs. 1–6, pl. 6 Fig. 69 [19].

Holotype

Type specimen represented by metabolically inactivated strain RCC1853 cryopreserved at the Roscoff Culture Collection (RCC; roscoff-culture-collection.org). Note: corresponds to isolate collected in the Ionian Sea, which is the type habitat of the original observation.

Habitat

Present in all oceans, in water with monthly sea surface temperature ranging from 0 to 25 °C.

Gephyrocapsa pseudohuxleyi sp. nov Bendif, Probert, Beaufort, Rickaby & Archontikis

Description

Coccoliths with moderately elevated distal shield (2–4.5 μm length) and elements of variable width (0.05–0.15 μm); inner tube with variable width, sometimes irregular, central area sometimes with a grill of curved rods, sometimes thick lath-like solid plates (sometimes with irregular holes). Regarded as a variant of morphotype A.

Genetic diagnosis

genetically distant from other species of Gephyrocapsa by genome sequences. Admixture pattern distinct from G. pujosae comb. nov. and G. huxleyi and forms the phylogenetic clade A2 (Fig. 1).

Holotype

Metabolically inactivated strain PLYM217 cryopreserved at the RCC as RCC1731.

Etymology

Based on the contraction of “pseudo-“ and huxleyi to highlight the partially cryptic relationship with the former species at the coccolith morphology level.

Habitat

Present in all oceans at low latitudes in water with monthly temperature ranging from 15 to 30 °C; can be found in temperate water during summer.

Gephyrocapsa pujosae Verbeek comb. & stat. nov. emend. Bendif, Probert, Young, Beaufort, Rickaby & Archontikis

Description

Coccoliths with distal shield elevated (2.5–5 μm length), composed of wide or narrow shield elements (0.05–0.12 μm width). Central area with thin lath like elements forming sometimes a thin solid plate, sometimes absent. Comprise previously described morphotypes B, B/C, C and O.

Genetic diagnosis

genetically distant from other Gephyrocapsa species by genome sequences. Admixture pattern distinct from G. pujosae comb. nov. and G. huxleyi and forms the phylogenetic clade B (Fig. 1).

Basionym

Emiliania pujosae Verbeek 1990 p. 23-24 pl. 1 Figs. 4–9.

Synonym

Emiliania huxleyi var. pujosae Medlin & Green 1996 [29, 30].

Holotype

Metabolically inactivated strain PLY92D cryopreserved at the RCC as RCC174.

Habitat

Present in all oceans, temperate and high latitudes, found in water with monthly sea surface temperature ranging from 0 to 17 °C.

Source: Ecology - nature.com