Genome-resolved biogeography of Phaeocystales, cosmopolitan bloom-forming algae

Abstract

Phaeocystales, comprising the genus Phaeocystis and an uncharacterized sister lineage, are nanoplanktonic haptophytes widespread in the global ocean. Several species form mucilaginous colonies and influence key biogeochemical cycles, yet their underlying diversity and ecological strategies remain underexplored. Here, we present new genomic data from 13 strains, including three high-quality reference genomes (N50 > 30 kbp), and integrate previous metagenome-assembled genomes to resolve a robust phylogeny. Divergence timing of P. antarctica aligns with Miocene cooling and Southern Ocean isolation. Genomic traits reveal metabolic flexibility, including mixotrophic nitrogen acquisition in temperate waters and gene expansions linked to polar nutrient adaptation. Concordantly, transcriptomic comparisons between temperate and polar Phaeocystis suggest Southern Ocean populations experience iron and B₁₂ limitation. We also identify signatures of horizontal gene transfer and endogenous giant virus/virophage insertions. Together, these findings highlight Phaeocystales as an ecologically versatile and geographically widespread lineage shaped by evolutionary innovation and adaptation to contrasting environmental stressors.

Fig. 1. Significance of Phaeocystis spp.

a Global abundance of Phaeocystis is in the range of 1–28 % of total marine eukaryotes, on par with well-known algal groups (e.g., diatoms and coccolithophores) and zooplankton. These figures highlight their importance in marine environments as primary producers. Estimates are based on biomass (MAREDAT), total occurrence in unigene collections, and based on genome-mapping of environmental reads (this study). Box-and-whisker plots within the violin plots here and in b show median, interquartile range, and 1.5*IQR values. b Environmental expression of organosulfur compound (DMSP and DMS) biosynthetic genes shows Phaeocystis-specific expression in comparison to other eukaryotic groups. Data from MATOU. c Generalized life cycle of colony-forming Phaeocystis (P. antarctica, P. globosa, and P. pouchetii) with four main morphotypes; two types of scale-forming haploid flagellates; a diploid cell cluster embedded in extracellular matrix, a large colony, and a naked diploid flagellate. In general, colonies form under nutrient-replete conditions in sufficient light (Brussaard et al.) and enclose photosynthetic non-flagellated cells. Haploid flagellates are associated with colony senescence and decline, are probably involved in sexual reproduction, and represent the life stage that persists through nutrient-deplete conditions. Other Phaeocystis species have been only found as solitary flagellates (P. cordata, P. rex, and P. scrobiculata; reviewed in Andersen et al.). Ploidy and typical morphotype sizes are indicated. d Estimates of lineage divergence times (node bars: 95% HPD) based on the concatenated phylogeny of 17 proteins (10,766 positions). Note the monophyly of polar strains in blue, coinciding with the last Antarctic reglaciation 12 Mya. e Pairwise sequence identity of best blast hits for protein models from Phaeocystis and other algal groups prevalent in the marine environment. Compared to diatoms (Pt × Tp, Phaeodactylum tricornutum vs. Thalassiosira pseudonana), chlorophytes (Cr × Ol, Chlamydonomas reinhardtii vs. Ostreococcus lucimarinus), and pelagophytes (Aa × Ps, Aureococcus anophagefferens vs. pelagophyte CCMP2097), Phaeocystis are a recently divergent group (Pa, Pc, Pg, P. antarctica, P. cordata, P. globosa, respectively). f Protein orthologous group overlap between Phaeocystis reference genomes and other algal groups. Source data are provided as a Source Data file.

Fig. 2. Biogeography of Phaeocystis spp. with respect to size fractions.

a The isolation sites of CCMP accessions and their biogeography, based on genome-wide Tara Oceans metagenomic read mapping, normalized to library size, expressed as reads per million (RPM). Data includes numbers for Phaeocystales MAGs. For each oceanic domain, multi-level pie charts are shown, with inner circles representing the contributions of each Phaeocystis spp. to the total read abundance, and outer circles representing the proportions of reads from different size fractions mapping to these species, clockwise according to the color legend. In most oceanic domains, Phaeocystis spp. occur as flagellates among nano- and picoplankton. b Copies per liter (CPL) abundance of 18S-V4 of Phaeocystis OTUs in the CalCOFI/NCOG data. c Log_e of the ratio of mitochondria- and plastid-mapping reads for Phaeocystis spp. where both organellar genomes are available. Left panel plots globosa genotypes, right all other species. Thicker lines mark the three species with the highest number of mitochondrial-mapping reads, P. antarctica, P. cordata and P. globosa genotype 2. Note a somewhat bimodal distribution for P. cordata, suggesting higher mitochondrial activity is condition specific. Source data are provided as a Source Data file.

Fig. 3. Phaeocystales functional profiles in temperate and polar biotopes.

a Venn diagram of Pfams for the top 1000 orthogroups with highest average expression (TPM) in each of the three analyzed biotopes (CCE, California Current Ecosystem; Arctic; SOC, Southern Ocean). b Patterns of significantly different enrichment among the three biotopes/datasets (nodes). Edges represent significant difference; colors correspond to colors in panels c-d. Numbers represent counts of annotated orthogroups and Pfams with respective significance patterns (p-adjusted Mann-Whitney test; Methods). c Relative normalized expression (TPM) of 7,316 orthogroups in three biotopes, colored by significance of expression difference. Circles and triangles mark orthogroups with and without Pfam annotations, respectively. d Relative normalized enrichment of selected Pfams in three biotopes, colored by significance of enrichment difference as in c. e Relative expression of B₁₂-dependent (metH) and B₁₂-independent (metE) methionine synthases (metH/(metH+metE)) globally and in the three biotopes from panels a–d. * – significant difference between biotopes (p-adjusted Mann-Whitney test; box-and-whisker plots within the violin plots show median, interquartile range, and 1.5*IQR values; Supplementary Note 5). f Relative expression of organosulfur biosynthesis enzymes methylthiohydroxybutyrate methyltransferase (dsyB) and DMSP lyase (Alma1) showing a relative enrichment in Alma1 expression near a Svalbard bloom. g Relative expression of iron-indicator markers ferredoxin (Fd, PF00111) and flavodoxin (Fld, PF00258). Flavodoxin, expressed in iron-limiting conditions, is widely utilized by Phaeocystis. h Expression of iron-responsive proteins, colored by oceanic region, and their trend lines in Southern+Atlantic Ocean and other oceanic domains (see legend in panel i). Values normalized to total Phaeocystales expression; error bands represent 95% CI to the corresponding linear regressions. Two-sided Spearman’s rho and p-values are shown in the upper left corner for the Southern+Atlantic Ocean (in blue, n = 49) and the other data (in orange, n = 82). The expression of all genes in panels e–h was normalized to Phaeocystales total. i Gene expansion of selected gene families as a function of latitude (* – genes with adjusted p-value < 0.001; Supplementary Note 5). Gene copies normalized to length and single-copy gene loci. CA, carbonic anhydrase; Nitrite red., nitrite-sulfite reductase; VIT, vacuolar ion transporters. Source data are provided as a Source Data file.

Fig. 4. Functional analysis of fast-evolving Pfam families in Phaeocystis spp. and other algal genomes.

a Global fractions of KOG biological functions in Phaeocystis. Functional groups (letters) correspond to panel b. b For each phylogenetic node, upper and lower semi-circles represent, respectively, the number of gained and lost genes belonging to fast-evolving Pfam families, classified by KOG biological functions. The total number of genes in analyzed algal genomes is shown for comparison, clearly distinguishing haptophytes Emiliania huxleyi and Phaeocystis spp. from other algae in terms of gene richness. Only non-overlapping gene models were used for Phaeocystis (i.e., not isoforms). Based on a multigenic analysis, the ultrametric timetree on the left indicates the approximate divergence times of shown species in millions of years (My). c Average abundance of genes (RPM, reads per million) belonging to fast-evolving Pfams from Tara Oceans metaT data, plotted against the number of gained/lost genes per respective family. Source data are provided as a Source Data file.

Fig. 5. Horizontal gene transfer (HGT) events in Phaeocystis draft genomes.

a The workflow overview (Methods). b–e, The left and right panels summarize the taxonomic composition for all detected HGT events (referred to as total, i.e., including cases where Haptophyta/Phaeocystales possibly donated genes) and those filtered by ingroup monophyly, respectively. Pa - P. antarctica, Pg - P. globosa, Pc - P. cordata. Note that the vast majority (98.2%) of Alveolata sequences detected in the HGT clades belonged to Dinoflagellata. Hundreds of additional candidate HGTs, including viral sequences, were found, though the direction of gene transfer to or from Phaeocystis could not be confidently inferred. b, d, Numbers of genes, where haptophyte sequences were present in the HGT clade, and these events likely pre-date the split of Phaeocystales from other Prymnesiophyceae. c, e Numbers of HGT genes exclusive to Phaeocystales. Stacked bars show the contributions of various lineages to these HGT events. f The number of HGT-originated orthologous groups gained and lost at each lineage of Phaeocystis. g Functional annotations of the total HGT pool. h Phylogeny of ice-binding proteins. Each tip represents a sequence with a taxonomic affiliation according to the color legend, or an unassigned sequence from the Tara MATOU v1 database (dark grey). Sequences with predominantly polar occurrence or found in polar algae are marked in the inner band; clades found in P. antarctica are marked by purple arrowheads. Dataset modified from, redundant sequences removed. Source data are provided as a Source Data file.