Regulation of sedimentation rate shapes the evolution of multicellularity in a close unicellular relative of animals

Significant increases in sedimentation rate accompany the evolution of multicellularity. These increases should lead to rapid changes in ecological distribution, thereby affecting the costs and benefits of multicellularity and its likelihood to evolve. However, how genetic and cellular traits control this process, their likelihood of emergence over evolutionary timescales, and the variation in these traits as multicellularity evolves are still poorly understood. Here, using isolates of the ichthyosporean genus Sphaeroforma-close unicellular relatives of animals with brief transient multicellular life stages-we demonstrate that sedimentation rate is a highly variable and evolvable trait affected by at least 2 distinct physical mechanisms. First, we find extensive (>300×) variation in sedimentation rates for different Sphaeroforma species, mainly driven by size and density during the unicellular-to-multicellular life cycle transition. Second, using experimental evolution with sedimentation rate as a focal trait, we readily obtained, for the first time, fast settling and multicellular Sphaeroforma arctica isolates. Quantitative microscopy showed that increased sedimentation rates most often arose by incomplete cellular separation after cell division, leading to clonal “clumping” multicellular variants with increased size and density. Strikingly, density increases also arose by an acceleration of the nuclear doubling time relative to cell size. Similar size- and density-affecting phenotypes were observed in 4 additional species from the Sphaeroforma genus, suggesting that variation in these traits might be widespread in the marine habitat. By resequencing evolved isolates to high genomic coverage, we identified mutations in regulators of cytokinesis, plasma membrane remodeling, and chromatin condensation that may contribute to both clump formation and the increase in the nuclear number-to-volume ratio. Taken together, this study illustrates how extensive cellular control of density and size drive sedimentation rate variation, likely shaping the onset and further evolution of multicellularity.


Introduction
The emergence of multicellularity from single-celled life represents a major transition, which has occurred many times independently across the tree of life [1][2][3][4][5][6][7][8]. Multicellularity can arise either by aggregation of single cells that come together or from single cells that are maintained together clonally after division [9][10][11]. The unicellular and intermediate multicellular ancestors, which led to present-day multicellular organisms, have long been extinct [3], obscuring direct investigation of how multicellular life has emerged. However, several strategies have been used to study the emergence of multicellularity, including the use of experimental evolution (EE) approaches and the investigation of novel non-model organisms at pivotal positions in the tree of life.
Alternatively, non-model organisms with key evolutionary positions can be used to better understand the emergence of multicellularity. In particular, the study of unicellular holozoans (Fig 1A), the closest unicellular relatives of animals, revealed that these organisms contain a rich repertoire of genes required for cell adhesion, cell signaling, and transcriptional regulation and that each unicellular holozoan lineage uses a distinct developmental mode that includes transient multicellular forms [2,[23][24][25][26][27]. For instance, the choanoflagellate Salpingoeca rosetta can form clonal multicellular colonies through serial cell division in response to a sulphonolipid of bacterial origin [28][29][30], whereas the filasterean Capsaspora owczarzaki can form multicellular structures by aggregation [31]. Ichthyosporeans display a coenocytic life cycle unique among unicellular holozoan lineages and pass through a short and transient clonal multicellular life stage prior to the release of newborn cells [32][33][34][35][36]. Despite their pivotal phylogenetic position, their rich "animal" genetic toolkit and the capacity to undergo transient multicellularity, to date, unicellular holozoans lineages and EE, have never been combined.
Beside the formation of new biological structures and increases in organismal size, the emergence of multicellularity is frequently accompanied with an increase in sedimentation rate. Indeed, S. cerevisiae snowflake yeast, multicellular C. reinhardtii, and S. rosetta colonies sediment faster when compared to their unicellular counterparts [15,16]. Such correlation has been described in several marine phytoplankton species where multicellular life stages show faster sedimentation than unicellular ones [38][39][40][41]. This phenotype has a large impact on where in the water column these microbes proliferate [38,[42][43][44] and thus is presumed to be under strong genetic control and selective pressure. Despite its capacity to affect the depth at which marine species flourish, the role of sedimentation rate, the potential impact of its variation, and its connection to the emergence of multicellularity has not been systematically analyzed across unicellular marine organisms, including in the closest unicellular relatives of animals with transient multicellular life stages. Here, using EE, we characterize how regulation of sedimentation rate can influence the emergence of stable multicellular life forms in the ichthyosporean Sphaeroforma genus, a close unicellular relative of animals.
Due to the nature of the coenocytic life cycle (Fig 1B), which is associated with an increase in the number of nuclei and coenocyte volume, we expected to observe an increase in cellular sedimentation rates over time [44,47,48].To better understand this relationship, throughout this study, we conducted overlapping experiments characterizing cell volume and sedimentation rates of Sphaeroforma species for cultures growing at 17˚C for 72 hours. For certain replicates of this core dataset, we included measurements of various genetic variants, the temperature dependence of phenotypes, fitness, as well as the speed of nuclear duplication. Overall, the measurements we report are highly reproducible with >95% variance explained for replicate measurements across phenotypes (S1A Fig and S1-S4 Data) (see Methods). This reflected a high heritability of the different phenotypes in a given environment.
To begin, we measured the sedimentation rate of 5 different Sphaeroforma species that have been isolated from different habitats either free-living or host-derived, namely S. arctica, S. sirkka, S. napiecek, S. gastrica, and S. nootakensis, using a quantitative sedimentation rate assay based on changes in optical density over time [32,45,[49][50][51][52]. According to these estimates, cell sedimentation rates varied greatly, from between 0.4 to 125 μm per second (i.e., up to 0.45 meters per day) ( Fig 1E). This broad variability over the life cycle and among Sphaeroforma species suggested that appreciable changes in size and/or cellular density should accompany different stages of the life cycle.

The transient multicellular life stage of S. arctica is associated with an increase in sedimentation rate
To better understand the cellular basis of sedimentation rate variation, we focused on S. arctica, the most studied Sphaeroforma species to date [32,34,50]. Using fixed cell imaging, we observed that synchronized cultures of S. arctica undergo a complete life cycle in about 48 hours and can reach up to 256 nuclei per coenocyte before undergoing cellularization and releasing newborn cells (Fig 1F). Prior to this release, all cellular division occurs in highly spherical mother coenocytes. Consistent with previous results [34], nuclear division cycles were periodic during coenocytic growth and occurred on average every 9 to 10 hours, as measured indirectly from changes in DNA content ( Fig 1G). Average cell volume increased throughout the coenocytic cycle to reach its maximum value at 36 hours prior to the release of newborn cells (Fig 1H). Similarly, the sedimentation rate increased up to >300-fold the initial value after 36 hours ( Fig 1I). Altogether, we observe that every cycle of nuclear division is associated with a significant increase in nuclear content, volume, and sedimentation rate with a distinct peak just prior to cell release. As the transient multicellular life stage of S. arctica occurs during the latest stage of cellularization and ahead of cell release [32], our results suggest that it is tightly associated with an increased sedimentation rate.

Experimental evolution of fast-settling mutants
Given how cell size and division across the cell cycle are regulated, we reasoned that these variable traits should also be heritable and, hence, evolvable. To test this, we conducted an evolution experiment to generate mutants with increased sedimentation rates [15,16,18]. Briefly, 10 independent cultures of S. arctica (S1 to S10) originating from the same ancestral clone (AN) (Fig 2A) were diluted in fresh MB medium at 17˚C (Fig 1C). Selection was performed every 24 hours by allowing the cultures to sediment in tubes for 2 minutes before transferring and propagating the fastest-settling cells in fresh medium (Fig 2A). EE was continued for 56 transfers (8 weeks) or about 364 nuclear generations, and a frozen stock was conserved every 7 transfers (1 week) (see Methods) (Fig 2A, S1 Table). Most lineages exhibited a distinct clumping phenotype to varying degrees across all evolved populations (Fig 2B). This phenotype has never been observed in any Sphaeroforma culture, despite variation of culturing conditions across labs, and the hundreds of routine passages without selection, thus suggesting that the appearance of such phenotype is linked to our selection regime. The clumps of cells were not maintained by ionic forces or protein-dependent interactions, however could be separated by mild sonication without leading to cell lysis (S2B Fig). To assess when sedimentation rates increased during the selection process for each evolved population, we synchronized cultures using mild sonication before dilution in fresh medium and allowed them to undergo a complete life cycle before measuring the sedimentation rate ( Fig 2C). We observed that populations S1, S4, and S9 had the highest sedimentation rates at the end of the evolution experiment ( Fig 2C). We observed a dramatic increase in sedimentation rate already after 14 transfers for population S1 (approximately 91 generations) (Fig 2C). To assess and compare the variability in sedimentation rates among evolved lineages, we isolated and characterized a single clone from each evolved culture. Our results show that evolved clones all settled significantly faster than the common ancestor but that there was stark variation in sedimentation rates at the end of the EE (Fig 2D). Isolates from lineages S1, S4, and S9 settled the fastest upon sedimentation (S1 Video), while clones from lineages S2, S3, S5, and S7 were intermediate, and lineages S6, S8, and S10 settled the slowest (Fig 2D). Taken together, our results show that we can rapidly and routinely evolve fast-settling mutants in S. arctica using EE but that the outcome is not uniform across lineages.
To assess whether fast sedimentation is adaptive in our EE setup, we independently competed 2 fast-sedimenting isolates, S01 and S03, against their common ancestor (AN). For this relative fitness assay, we mixed monocultures of either S01 or S03 with AN at an approximately 1:1 ratio and subjected them to the same selective regime as during EE. Parallel cultures of each 1:1 mix grown and maintained without sedimentation selection served as experimental control (see Methods). Our results show that the proportion of clumpy S01 and S03 cells both significantly increased over time (a log2 selection coefficient of 0.23 and 0.13 per day for S01 and S03, respectively, p < 1e-7), suggesting that this trait was adaptive in our experimental setup (Figs 2E, S1A, and S2D; S4 Data). In contrast, S01 and S03 cells decreased in frequency relative to the AN strain in the control growth experiment without sedimentation (a log2 selection coefficient of −0.13 and −0.15 per day for S01 and S03, respectively, p < 1e-7), suggesting that the heritable changes that these strains bore were deleterious in normal laboratory propagation conditions (Figs 2E and S2D).
We next sought to characterize the cellular mechanisms giving rise to variation in cellular sedimentation rates in evolved clones. Using quantitative microscopy, we show that clump perimeter and average number of cells observed per clump correlated highly with sedimentation rates (Figs 2E, S2C, and S2E). Indeed, across all experiments reported in this study, we found that using Stokes' law, a single globally fitted density parameter and our size measurements, could explain the majority of variance in sedimentation rate (R 2 = 0.69, RMSE (as a proportion of the range in observations) = 18%).

Fast-sedimenting mutants form clonal clumps
Earlier, we defined 3 distinct developmental stages of S. arctica life cycle: (i) cell growth with an increase in coenocyte volume; (ii) cellularization, which coincides with actomyosin network formation and plasma membrane invaginations; and (iii) release of newborn cells (Fig 1B) [32]. A key developmental movement named "flip," defined by an abrupt internal morphological change in the coenocyte, can be used as a reference point to characterize life cycle stages. Prior to this event (pre-flip), actomyosin-dependent plasma membrane invaginations occur, while afterwards (post-flip), the cell wall is formed around individual cells prior to their release from the coenocyte [32]. Using time-lapse microscopy at 12˚C, we observed that clump formation in fast-settling mutants coincides with cell release (Fig 3A, S2 Video).
Importantly, no new cell aggregation processes were detected after cell release, suggesting that all clumps formed either prior to or concomitantly with the release of newborn cells (S2 Video). To examine whether clumps formed due to defects at the level of either plasma membrane or cell wall, we stained plasma membranes using FM4-64 and cell walls using calcofluorwhite. We found that the process of plasma membrane invaginations during cellularization appears to be unchanged prior to flip (Fig 3B, S3 Video) and that all newborn cells, even in the clumps, have a distinctive cell wall surrounding them (Figs 3C and S3A). These results show that clumps are formed post-flip in fast-settling mutants.
As no cells appear to aggregate after release of newborn cells, our results suggest that clumps are maintained together in a clonal form. To further confirm this, we sonicated both the ancestor (AN) and S1 clumps and stained them separately with 2 distinct fluorescent dyes prior to mixing them in a 1:1 ratio (S3B  and each only contained a single AN cell trapped inside a smaller S1 clump. In a control experiment with two differentially stained S1 cultures, we observed an almost identically small number (i.e., approximately 1.7%) of mixed clumps (S3B Fig); hence, it appears that mixes of cells happen only sporadically at low frequencies by random association, irrespective of genotype. Altogether, our results suggest that the evolved clump phenotype is not a result of spontaneous cell aggregation and instead arises from incomplete detachment between cells.
We next examined the kinetics of cellularization and the process by which cells propagate as part of clumps. We counted the number of cells detaching from clumps at cell release and observed variable detachment across all 3 fast-settling mutants. The clumpiest isolate, population S1, was the least prone to cell detachment ( Fig 3A and 3D, S2 and S3 Videos). Intriguingly, despite exhibiting similar sedimentation rates, S4 and S9 clumps showed a higher detachment frequency than S1 (Fig 3A and 3D). Image analysis of time-lapse movies also showed that life stage durations varied among fast-settling mutants at 12˚C, with S1, S4, and S9 initiating cellularization and undergoing flip, respectively, 5.5, 7.5, and 10.5 hours earlier than the ancestor (S3C and S3D Fig). Post-flip duration also varied significantly among mutants, with S1 and S9 requiring more time to release cells compared to S4 (S3C and S3D Fig). While many aspects of the replication cycle dynamics were variable for these mutants, the duration of cellularization was fairly invariant (approximately 9 hours) (S3C and S3D Fig). Finally, measurements of coenocyte volume show that S4 and S9 coenocytes undergo flip at substantially smaller volumes (approximately 1AU : PerPLOSstyle; alwaysuseperiods; notcommas; toindicatedecimalpoints:Hence; }1 .8 and approximately 3.3 times smaller, respectively) compared to AN or S1 ( Fig 3E). These results show that, despite their shared capacity to clump and similar sedimentation rates, fast-settling mutants exhibit significant variability in their life cycle dynamics, with S4 and S9 mutants initiating cellularization earlier, dispersing from clumps with a higher frequency, and undergoing flip and cell release at smaller coenocyte volumes compared to the S1 mutant.

Increased nuclear number-to-volume ratio leads to faster sedimentation
Above, we observed that S4 and S9 mutants can sediment as fast as S1 despite their smaller coenocyte volumes, suggesting an alternative regulation mechanism of sedimentation rate. Across all experiments reported in this study, we found that approximately 31% of the variance in observations of sedimentation rate could not be explained by cell size only, suggesting that cellular density might also contribute to this variation [37,40,41]. From Stoke's law, we calculated that excess cellular density, i.e., cellular density minus that of distilled water (1,000 kg/ m 3 ), might vary between 40 and 300 kg/m 3 for S. arctica wild-type and evolved clones across their life cycle-the upper limits between the densities of pure protein and pure cellulose. Values reached approximately 650 kg/m 3 for wild S. nootakensis (Snoo) soon after cellularization, approaching the excess density of pure nucleic acid ( Fig 1D). During cell cycle stages prior to cellularization, when cells were most spherical (<36 hours), excess density varied from 40 to 200 kg/m 3 across Sphaeroforma isolates.
To better characterize the relationship between sedimentation rate, cell cycle, and size, we performed higher resolution measurements of sedimentation rates over the complete life cycle of the ancestor and all 3 fast-settling mutants at 12˚C and 17˚C. Consistent with their capacity cells are still separated by a cell wall inside fast-settling clumps, suggesting that clumps are formed post-flip. Bar, 50 μm. (D) Number of cells detaching per clumps at cell release and for the following 3 hours, measured from timelapse movies (n = 50 coenocytes at cell release per strain). (E) Coenocyte volume at flip, measured from time-lapse movies, shows that S4 and S9 are significantly smaller when compared to AN and S1 (n = 50 coenocytes at flip per strain). The data underlying this figure may be found at https://doi.org/ 10 to form clumps, we observed that the sedimentation rate of all fast-settling mutants increases during growth but, unlike the ancestor, does not recover after cell release to their original levels (Figs 4A and S4A). Interestingly, we noticed that individual S4 and S9 coenocytes sediment faster (approximately 2AU : PerPLOSstyle; alwaysuseperiods; notcommas; toindicatedecimalpoints:Hen .5× and approximately 1.6×, respectively) than S1 or AN even before clump formation (24 to 36 hours time points) (Figs 4A and S4A). Such increase in sedimentation rate was not due to a rise in cell size or change in cell shape as both S4 and S9 exhibit smaller cell perimeters throughout the cell cycle (Figs 4B, 4C, S4B, and S4C). Rather, excess cellular density estimations show that both S4 and S9, even prior to cell release and clump formation, tend to be on average 3× denser when compared to the ancestor (Figs 4D and S4D). Altogether, these results show that both cell size and cell density contribute to sedimentation rate variation in S. arctica.
As cell size and nuclear division cycles are decoupled in S. arctica [34], we reasoned that increased cell density in S4 and S9 could be caused by an acceleration of nuclear divisions leading to a rise in the number of nuclei per volume. Using DAPI staining to label nuclear DNA, we observed that S4 and S9 undergo nuclear duplication faster (approximately 2 hours) than both AN and S1 (Figs 4E, 4F, and S4E-S4H). By carefully examining the volumes of coenocytes containing the same number of nuclei at the single-cell level, we show that for the same nuclear content, S4 and S9 tend to be 30% to 45% smaller in volume when compared to the ancestor (Figs 4G and S4I). Consequently, both S4 and S9 exhibited the highest number of nuclei per volume (nuclear number-to-volume ratio) (Figs 4H and S4J). Taken together, these results argue that cell density can contribute an appreciable amount to cellular sedimentation rates (up to approximately 50 μm/s) and that mechanistically this could arise by faster nuclear doubling times relative to cell size.

Evolved genetic variation correlating with fast sedimentation
Up to now, our results suggest that S. arctica mutants evolved faster sedimentation using 2 strategies: (i) clump formation; and (ii) increased nuclear number-to-volume ratio. We found that sedimentation rate variation was highly heritable, persisting for 780 generations of passaging for all 10 isolates without selection for sedimentation phenotype, suggesting that the phenotypes have a genetic basis. To test this, we resequenced the whole genomes of both the ancestral clone (AN) and 1 evolved clone per lineage (S1 to S10) obtained at the conclusion of the evolution experiment (Week 8) with high (>30-fold) coverage of the very large genome size of S. arctica, at 143 Mbp. Following very careful variant filtering (S2 and S3 Tables), we identified a total of only 26 independently evolved variants with an average of 2.6 mutations per clone (range 1 to 5 per clone) (S4 Table). Of the 26 variants, 24 (approximately 92.3%) were SNPs (11 coding, and 13 intergenic or intronic), and 2 were insertions (1 coding, 1 intergenic) (S4 Table). Beside the fact that 2 of these variants are identical SNPs at the same position (Sarc4_g3900) and have independently evolved in 2 different backgrounds (S4 Table), no other form of genetic parallelism was detected. While it is very likely that some of these mutations hitchhiked in the background of beneficial driver mutations, we found that the coding SNP variants were clearly skewed toward nonsynonymous changes (9:2), with a cumulative dN:dS ratio of 1.32. This indicates the general presence of positive selection and, hence, adaptative evolution driving this molecular pattern. Based on COG assignments (S6 Table), 4 of the changes are orthologous to genes implicated in signal transduction, and 2 are related to genes with functions in DNA binding or chromosome condensation.
In the absence of molecular genetic tools and functional knowledge of many of the mutated gene targets, we set out to better understand how the distinct genetic variants could have influenced sedimentation rates and clump formation by examining the predicted expression shows that fast-settling mutants increase their size upon cell release. Every trace represents an independent experiment (n > 180 measurements per time point for each independent experiment). (C) Average perimeter of fast-settling cells and clumps at 24 and 60 hours, respectively, show that S4 and S9 cells and clumps have a smaller size when compared to S1. Every square represents an independent experiment, and the white circle represents the median (n > 180 coenocytes per time point for each independent experiment). (D) Excess cellular density of fast-settling individual coenocytes (before cellularization) and clumps (after cell release) at 24 hours and 60 hours, respectively, show that S4 and S9 single cells are denser when compared to S1 and AN. Every square represents an independent experiment, and the white circle represents the median. (E) Quantification of mean DNA content per time point for fast-settling mutants grown in MB at 17˚C. Every trace represents an independent experiment (n > 400 coenocytes per time point for each independent experiment). (F) Nuclear doubling time, calculated by linear regression of mean nuclear content at time points from 0 hour to 24 hours. Every square represents an independent experiment, and the white circle represents the median (n > 400 coenocytes per time point for each independent experiment). (G) Boxplots of cell volume measurements of DAPI-stained fixed cells. For 1-, 4-, 16-, and 64-nuclei cells. Cells with 1 nucleus represent newborn cells at the end of the experiment (n > 80 coenocytes per DNA content). (H) Boxplots of nuclear number-to-volume ratio of DAPI-stained cells show significant increase for S4 and S9 fast-settling mutants. Every square represents an independent experiment, and the white circle represents the median (n > 600 coenocytes per strain). (I-K) Temporal transcript abundance of genes mutated in fast-dynamics of mutated genes across the cell cycle. The data were derived from a recently published time-resolved transcriptomics dataset of the S. arctica life cycle [32]. Of the mutated genes, 9 showed no expression during the native life cycle, whereas 12 displayed dynamical expression during cellularization, and the remaining 5 genes were more or less stably expressed (Figs 4I-4K and S4K, S5 Table). We also annotated all mutation-associated genes based on a recent comprehensive orthology search [53] (S6 Table).
The fastest-settling and clumpiest mutant isolated from the population S1 bore a synonymous mutation of a homolog of increased sodium tolerance 1 (Ist1) superfamily and as such likely impacts gene expression rather than gene function. This gene shows a dynamic expression during cellularization and codes for a conserved protein involved in multivesicular body (MVB) protein sorting ( Fig 4I, S6 Table) [54,55]. In humans, hIST1, also known as KIAA0174, is a regulator of the endosomal sorting complex required for transport (ESCRT) pathway and has been shown to be essential for cytokinesis in mammalian cells [56]. Similarly, Ist1 orthologs in both budding and fission yeasts play a role in MVB sorting pathway and, when deleted, exhibit a multiseptated phenotype consistent with a role in cytokinesis and cell separation [57,58].
In the clone derived from population S4, we observed 5 distinct mutations, 2 nonsynonymous, 1 intergenic, and 2 intronic SNPs. Among the 2 nonsynonymous SNPs, one causes a E90G change in a homolog of human Kanadaptin (SLC4A1AP), which may play a role in signal transduction [59,60]. Among the noncoding SNPs, 1 mutation is found in an intron of the 7-dehydrocholesterol reductase (DHCR7), expressed during cellularization and known to be key in the cholesterol biosynthesis pathway [61,62], and the second intronic mutation codes for a STE20-like kinase (SLK), which plays numerous roles in cell cycle signaling and actin cytoskeleton regulation (Fig 4J) [63][64][65][66].
Finally, among the 5 mutations discovered in the clone derived from population S9, 2 mutations are in transcription factors that are continually expressed during the cell cycle: an intronic SNP in a basic helix-loop-helix (bHLH) transcription factor, and the sole nonsynonymous SNP leading to A923V change in a gene predicted to encode a nucleotide binding C2H2 Zn finger domain [67]. A third mutation was found in an intron of the highly and dynamically expressed homolog of the regulator of chromosome condensation 1 (RCC1; Fig 4K) [68][69][70]. RCC1 is a chromatin-associated protein implicated in several processes including nuclear formation, mRNA splicing, and DNA replication [69,[71][72][73][74]. Thus, it may contribute to the accelerated nuclear duplication cycle observed in S9 by impacting cell cycle progression. Altogether, the mutations identified in both S4 and S9 may affect both cellularization and cell separation.
Among the variants detected in evolved clones with intermediate-settling phenotype, we highlight an intergenic mutation 118 bp downstream of Dynamin-1 known to be essential for cytokinesis across different taxa [75][76][77], and a nonsynonymous mutation in a protein similar to Fibrillin-2 (Sarc4_g7365T), which is an extracellular matrix (ECM) glycoprotein essential for the formation of elastic fibers in animals (S4K Fig) [78][79][80]. Altogether, our results across all isolates suggest that a large mutational target affects cellular sedimentation and multicellularity.

Sedimentation rate variation across Sphaeroforma species is driven by cell size and density
Lastly, we examined whether the variation in cell sedimentation observed across different Sphaeroforma species (Fig 1E) could also be explained by clumping or increased nuclear settling phenotype across the native life cycle of S. arctica. The data underlying these graphs may be found in S5 number-to-volume ratio. To do so, we investigated the life cycle dynamics, coenocyte volume, and nuclear duplication time among distinct Sphaeroforma sister species. To date, 6 different Sphaeroforma species have been isolated either in a free-living form or derived from different marine hosts: S. arctica, S. sirkka, S. napiecek, S. tapetis, S. gastrica, and S. nootakensis (Fig 5A) [32,46,[49][50][51][52]. Using previously established growth methods for S. arctica combined with live and fixed imaging, we first observed that all sister species but S. tapetis show a synchronized coenocytic life cycle (Figs 5B and S5A, S4 Video). Similar to abovementioned results with fastsettling mutants, sedimentation rate variations (Fig 1E) could be explained by variations in both cell size and cellular density. Indeed, we first observed that both S. gastrica and S. nootakensis occasionally form clumps, exhibit a lower frequency of cell detachment, and thus have an increased cellular density after cell release compared to the other Sphaeroforma species (Figs 1C, 5B, and 5C, S4 Video). Additionally, despite their increased sedimentation rate ( Fig  1C), we found that all sister species exhibited approximately 20% to 45% smaller coenocyte size prior to cell release when compared to S. arctica, which reveals an increased cellular density (Figs 5D, S5B, and S5D, S4 Video). Similar to the fast-settling mutants S4 and S9 above, the increase in cellular density was associated with an acceleration of the nuclear division cycles and the subsequent rise in the nuclear number-to-volume ratio (Figs 5D-5F and S5D-S5F). Notably, S. sirkka and S. napiecek, both previously isolated as free-living [51,52], exhibit an increase in nuclear number-to-volume ratio but no ability to form clumps. Altogether, our results show that, similarly to experimentally evolved strains, fast sedimentation variation could occur by both clump formation and/or increase in the nuclear number-to-volume ratio for Sphaeroforma species (Fig 5H) and thus might represent in itself a widespread and highly variable phenotype in the marine habitat.

Discussion
Here, we performed, for the first time, EE on one of the closest unicellular relatives of animals, and we demonstrate that, under suitable selection pressure, the ichthyosporean S. arctica can evolve stable multicellularity. In particular, we observed the independent rise of clump formation, and faster settling phenotypes across populations within less than 400 generations. The precise detectable onset of the phenotypes varied across lineages and occurred as early as approximately 91 generations in lineage S1. We further show that faster sedimentation phenotypes are highly heritable and, according to direct competition assays with the common ancestor, are expected to be adaptive under the environmental conditions of the experiment. Our results add to previous observations of the rapid emergence of multicellularity in yeast and green algae, which all can evolve multicellular clump-forming structures within short evolutionary timescales [15,16,18]. As ichthyosporeans proliferate through an uncommon coenocytic life cycle [32,33,35,45,46,51], our results show that this evolutionary process of selection for faster sedimentation is accessible at microevolutionary timescales across taxa and organisms with highly diverged modes of proliferation.
In this study, we show that all fast-settling S. arctica cells increased their cell size by increasing cell adhesion post-cellularization, leading to the formation of clumps (Fig 5H). Such results are analogous to cell cluster formation in snowflake yeast and C. reinhardtii, which arises through incomplete separation of mother and daughter cells [15][16][17]. Altogether, these results suggest that regulation of sedimentation rate can constrain unicellular species to generate multicellular cell phenotypes by increasing their cell adhesion efficiency. However, we found that approximately 31% of the variance could not be explained by cell (clump) size only. Indeed, 2 fast-settling mutants (S4 and S9) exhibited an increase in sedimentation rate prior to clump formation, which was associated with an accelerated nuclear division cell cycle leading to an increase in the number of nuclei per unit volume (Figs 4G and 5H). Previous results have shown that, in S. arctica, both nuclear duplication cycles and cell size are uncoupled [34]. Our results support these findings and indicate that nuclear division cycles and cell size could be regulated separately, allowing adaptive change in either, and independently of one another.
By analyzing the genomes of evolved isolates, we identified an overabundance of nonsynonymous mutations, indicating positive selection, with at least 1 and up to 5 independent mutations in each lineage. Many of the better-characterized genes that carry mutations in either coding or intergenic regions are dynamically expressed during cellularization. Several mutations were found in genes coding for cytoskeletal regulators and cytokinesis proteins, which is consistent with previous studies in which cytokinesis-deficient mutants were often associated with an incomplete cell separation across taxa [81][82][83][84][85]. Other mutations were found in genes involved in cell signaling, plasma membrane remodeling, and chromatin condensation regulators reflecting a large and accessible mutational target affecting sedimentation rate phenotypes. Cytokinesis defects may lead to incomplete cell separation, which may explain part of the cell clumping phenotype observed in the fast sedimentors. Given that similar phenotypes emerged independently multiple times with no genetic overlap, we conclude that the mutational target for these traits could be quite large for S. arctica, opening the possibility for variation and evolution in multicellularity-related phenotypes.
We found that closely related species exhibit widespread variation in both clump formation and nuclei-to-volume ratio. Given the evolvability of the trait and the few natural samples, we examined it is difficult to argue which cellular states affecting sedimentation rates are ancestral. For example, considering the species tree and the trait of clumpiness (Fig 5A and 5B, S4 Video), only 3 Sphaerforma species exhibited this behavior under our experimental conditions: S. gastrica, S. tapetis, and S. nootkatensis. Despite characterizing for the first time all the available isolates of Sphaeroforma sp., the number of studied isolates remains limited (6 species), and, thus, we cannot distinguish whether this trait is "ancestral" or "novel." Rather, our experimental data suggest that clumpiness and density could evolve over short microevolutionary timescales such as those we have measured in the lab and this in both directions. Evolved lines with highly clumping phenotypes had the competitive edge over their nonclumping common ancestor in relative fitness assays. In contrast, during exponential growth and in the absence of selection, evolved clumpy line were selectively disfavoured relative to their common ancestor, which suggests that clumpy phenotypes are prone to rapid replacement by less clumpy or dense morphs when environments change.
Indeed, marine organisms have evolved various passive or active means of maintaining their position in the water column, for example, using motility and/or ingenious approaches to regulate buoyancy [44,[86][87][88][89][90]. Sphaeroforma species are remarkably spherical, immobile, lack flagella, and yet exhibit a substantial increase in cell size and density over the life cycle, thus representing a challenge to maintaining buoyancy in marine habitat. This work establishes that Sphaeroforma's cell size and density are subject to tight cellular control and are highly evolvable traits. Taken together, these observations suggest that sedimentation rate is a highly evolvable trait, which itself likely shapes the gain and loss of multicellularity.

Culture conditions
All Sphaeroforma sp. cultures were grown and synchronized as described previously for S. arctica [32,34]. Briefly, saturated cultures in MB (Difco BD, NJ, USA; 37.4 g/L) were diluted into fresh medium at low density (1:250 dilution of the saturated culture) and grown in rectangular canted neck cell culture flask with vented cap (Falcon; ref: 353108) at either 17˚C or 12˚C, resulting in a synchronously growing culture. Saturated culture of Sphaeroforma sp. are obtained after 3 weeks of growth in MB.

Experimental evolution
Ten replicate population (S1 to S10) of genetically identical S. arctica (AN) were first diluted 250-fold in 5 ml of MB and grown in rectangular canted neck cell culture flask with vented cap (Falcon; ref: 353108) at 17˚C. Cells were grown at 17˚C rather than the previously used 12˚C in order to increase growth rates and accelerate evolutionary outcome. Every 24 hours, the entire population was transferred into a 15-ml falcon tube and allowed to sediment for 2 minutes on the bench. A 5-ml pipette was then positioned vertically and used to collect 500 μl of cell culture from the bottom of the falcon tube. The cells, still vertically positioned in the pipette, were then allowed to sediment once more for 15 seconds before the transfer of a single drop, equivalent to 20 μl, into 5 ml of fresh MB (approximately 250× dilution). Every 7 transfers, a frozen fossil was conserved by adding 10% of DMSO to 1 ml of culture and preserved at −80˚C. Single clones of each replicate population (S1 to S10) were obtained at the end of week 8 by serial dilutions. Since S. arctica grows as a coenocyte, a temporal generation is not defined by a complete coenocytic cycle but is equivalent to the doubling of the number of nuclei. We estimated the nuclear doubling time by measuring the number of cells and the number of nuclei per cell for each transfer separately (S1 Table). Briefly, the entire EE experiment comprised approximately 364 generations, in which all populations underwent a total of 28 complete coenocytic cycles. Each coenocytic cycle included two 24 hours subpassages and comprised a total of approximately 13 doublings (S2A Fig, S1  Table). Importantly, the effective population size was kept in a very narrow range across subpassages (at approximately 10 5 ) and thus over the entire experiment (S2A Fig, S1 Table). Therefore, in evolutionary terms, the population size was consistently high enough to favor of natural selection over random evolution throughout the course of the experiment. This assumption was reconfirmed genetically by deriving a dN:dS ratio >1 from sequencing data (see main text).

Sedimentation rate measurements
Sedimentation rate was measured for Sphaeroforma sp. every 12 hours for a total of 72 hours unless indicated otherwise. To ensure reproducibility and homogeneous results, saturated Sphaeroforma sp. cultures were sonicated prior to the dilution in fresh MB media (250-fold dilution) using a Branson 450 Digital Sonifier (3 pulses of 15 seconds, 10% amplitude). For each measurement, obtained either from different stages of the cell cycle or from different Sphaeroforma species, 1 ml of cell culture was added into a disposable plastic spectrophotometer cuvette (semi-micro, 1.5 ml) and homogenized by vortex. Optical density (OD 600 ) was measured using an Eppendorf Biophotometer (Model #1631) at T = 0, corresponding to the first time point after placing the cuvette in the spectrophotometer. The OD 600 was then continuously measured every 30 seconds for 3 minutes while cells were slowly sedimenting in the cuvette. To ensure that OD 600 measurements stayed within the detection limits of the spectrophotometer, early life stages (T0-T48) were not diluted in the cuvette, whereas later life stages (T60-T72) were diluted 1/100 in fresh MB media.
For assessing clump dissociation in S2B Fig, AN or S1 cultures were incubated for 2 hours at 37˚C in MB, artificial sea water (ASW) (Instant Ocean, 36 g/L) with different salt concentrations (18 g/L for 0.5X and 72 g/L for 2X) to assess any effect of electrostatic forces, phosphatebuffered saline (PBS) 1X (Sigma-Aldrich) with either Proteinase K at 200 μg/mL final concentration (New England Biolabs, Ipswich, MA, USA) to assess for any protein-dependent effect or sonication using a Branson 450 Digital Sonifier (3 pulses of 15 seconds, 10% amplitude). Only sonication resulted in dissociation of the clumps and a reduction in sedimentation rates.

Maximum likelihood estimation of sedimentation velocity and cellular density based on OD600 sedimentation rate assay
To briefly summarize, we related our OD600 A.U. sedimentation rate measurements and radius measurements (S1 Data) to previously published datasets [37] and [91] (S2 and S3 Data) to gain an estimate of our sedimentation rate measurements in metric units, which allowed maximum likelihood estimation of cellular densities based on these estimates (S1 Data).
We started by estimating the average radius of cells from perimeter measurements as We estimated sedimentation velocity in our measurements by assuming OD600 (OD) changes at a constant proportional rate of change with respect to time t = 0 and t = 1 by This yielded a rate of change in arbitrary distance units per second, which we next sought to relate to metric distance units. For this, we turned to the datasets from [40], compiled along with data from other studies by [37] in his Appendix For each of these observations, we calculated seawater density in experimental assays as SC + salinity where salinity constant SC = 35.16504/35 and the salinity is the salt content in g�l −1 [91]. We estimated the specific gravity or density of media used in each sedimentation rate measurement based upon the dataset published by [91] (S3 Data), using a second-order polynomial function of their density measurements as a function of salinity and temperature using R's lm() function. The error of this estimate is extremely low (<4e-3 kg�m −3 ) and was not propagated downstream.
We next calculated the density p p of cells across observations based on sedimentation velocity V in m�s, cell radius R in m, and media density p f by rearrangement of the terminal velocity equation where dynamic viscosity μ = 0.00109 Pa�s and gravitational acceleration g = 9.780 m�s 2 . This yielded a median phytoplankton excess density (p p -1,000) of 139 kg�m −3 with a range of 30 to 1,300 kg�m −3 . Some excess density estimates exceeding protein (220 kg�m -3 ) and cellulose (500 kg�m −3 ) have previously been suggested to arise by calcide in diatoms [92]. We concluded that this method of estimating cell density yielded similar values to those published by, for example, [92] and [40]. For downstream analyses, we excluded obvious outliers, including measurements of samples with densities exceeding that of cellulose, and measurements of D. rex, a large diatom, which, in the Smayda dataset, exhibited extraordinarily low sedimentation rates given their size.
To estimate the sedimentation velocities of our dataset, we assumed that our data (S1 Data) would fall within the typical measurement in the Smayda dataset (S2 Data). For the outlierexcluded subset of the Smayda dataset, we calculated an expected sedimentation velocity for what we would measure in our experimental setup based on the specific gravity of the seawater formulation we used in our measurements (37.4 g/l or a solvent density of 1,028.9 kg�m −3 ). We then used least-squares minimization to estimate 2 parameters: a scalar S of the AU velocity measurements, dODdt, which could best match our dataset with these expected velocities, and an average density parameter p p_hat , which could predict these values and the Smayda dataset's values based on Stokes' law. For this, we minimized the loss function where V 0 is either the seawater density-corrected velocity from Smayda or the dODdt parameter we calculated above, and H is the one-hot binary scalar in [0,1] corresponding to whether the i'th data point in our N observations was from Smayda's or our dataset, respectively. We used R's optim() function with the "L-BFGS-B" method with initial parameters values of S = 0.03 and p p_hat = 100 and repeated this fit for 500 10-fold (10% out-of-bag) bootstrap samples of individual observations across our full dataset and Smayda's to gain an estimate of the error on the parameters. S and p phat fits were positively correlated across testing/training folds (Pearson's ρ = 0.97); however, mean values were fairly limited, with S = 0.028 ± 0.002 and p p_hat = 78 ± 6.4 kg�m −3 (with error equal to the standard deviation of estimates across 500 folds). These narrow estimates indicated that the fits were reasonably well defined by the underlying dataset. The means and standard deviations across predictions of out-of-bag samples are what we have reported as means and standard error in S1 Data and propagated along with interreplicate and batch error reported in the figures and text.

Percentage of clumps measurements
To measure percentage of clumps across all mutants and transfers in S2C Fig, we first sonicated saturated S. arctica cultures using a Branson 450 Digital Sonifier (3 pulses of 15 seconds, 10% amplitude). Cells were then diluted in fresh MB (1:250 dilution), and cell concentration was then measured using a hemocytometer. Approximately 20 cells were then transferred per well in a 96-well plate. Cells were then monitored every 24 hours using an inverted optical microscope, and the percentage of clumps observed after cellularization was measured manually using a tally counter. This experiment was performed 3 independent times, and error bars are standard deviations.

Head-to-head competition
To perform the fitness experiment (head-to-head competition), we used saturated and sonicated cultures of the ancestor (AN), S01 and S03 evolved isolates. Following the measurement of the cell concentration using a neubauer chamber, we independently mixed both S01 and S03 with the AN at an approximately 1:1 ratio and subjected them to either the same selective regime as the evolution experiment (see above) or a control culture that was simply maintained in exponential growth without sedimentation selection. The experiment was maintained for 72 hours with passages every 24 hours. To quantify the outcome of the experiment, every 24 hours for 72 hours (7.5 to 12.5 generations), after dilution of the competition culture, we counted the number of clumpy and nonclumpy cells at each passage. For this, a fraction of the cultures was sonicated and highly diluted to obtain fewer than 10 cells per well of a 96-well plate. These cells were allowed to undergo a full life cycle in order to count by microscopy the proportion of cells, which bore the clumpy phenotype. The high penetrance of the clumpy phenotype in these clones therefore allowed not only a quantification of the clumpy phenotype itself in competition, but also the measurement of overall clonal fitness of these derived isolates.
Each data point for the fitness assay was a count of cells with an associated binomial error. In Figs 2E and S2D, this error is represented as 95% confidence intervals calculated using the binom.logit() function from the R package "binom" (Sundar Dorai-Raj (2014). binom: Binomial Confidence Intervals for Several Parameterizations. R package version 1.1-1. https:// CRAN.R-project.org/package=binom).
The fitness measurements began from highly synchronized cultures, and, therefore, reporting a per-generation fitness value was clouded by a sudden jump in cellular counts that occurred at 48 hours. Therefore, for simplicity and given the clear time-dependent changes in clumpy phenotype frequencies, we report fitness in terms of log2 fractional change per unit time ("log2 selection coefficient" in days). To obtain these numbers, we first used R's rbinom() function to resample the dataset 500 times, with each resampling using each sample's original measurement of clumpy cell count, total cell count, and proportion clumpy. For each resampling, after calculating the new value for proportion clumpy, we logit2 transformed the data and rescaled to the initial measurement with the following equation: where p is the proportion of the population that was clumpy, the tx underscore is each time point, and the t0 underscore is the time point of the initial proportion measurement in units days. The slope of the least-squares fit for this quantity with respect to time was obtained by R's lm() function, yielding the per day change for each resampled biological replicate (the "log2 selection coefficient"). The mean and standard deviation of this coefficient across all resampled replicates were taken as the fitness point estimate and standard error, respectively, and the 95% confidence intervals for Fig 2D is that standard error times 1.96.

Heritability measurements
To obtain a coarse understanding of how heritable our phenotypes were, we used a statistical genetic approach to quantify the heritability of traits, or what proportion of the variance in phenotype is due to an individual's inherited state. Briefly, the total Variance of a measured phenotype V_P(Total) can be variance partitioned between Environment, Genetics (or other inherited factors), and technical error or stochastic events (epsilon): The different partitioning terms in V_P(T) can be expressed as a proportion of V_P(T): And, therefore, we can refer to the proportion of variance due to genetics as pV_P(G).
In the lab, we eliminated environmental variance with tightly controlled experimental conditions such as temperature and the number of hours a measurement was taken after cell cycle synchronization. This allowed us to define heritability therefore as the environment-controlled genotypic variance: For this, we determined H for DNA content, perimeter length, sedimentation rate, and fitness phenotypes by calculating the variance explained by mean phenotype values within distinct genotype/environment combinations (S1A Fig). The results show that H exceeds 95% across phenotypes, and across the entire dataset, H exceeded 99% of the total phenotypic variance (ANOVA F = 1,118 on 252 and 735 DF, p = 0). This means that for a typical individual genotype in a given environment, we could predict its average phenotypic measurement with >97% accuracy.

Microscopy
Microscopy of live and fixed cells was performed using a Zeiss Axio Observer Z.1 Epifluorescence inverted microscope equipped with Colibri LED illumination system and an Axiocam 503 mono camera. An EC Plan-Neofluar 40x/0.75 air objective was used for images of fixed cells, and an N-ACHROPLAN 20x/0.45na Ph2 air objective was used for all live imaging, unless indicated otherwise.

Cell fixation and staining
Throughout this study, saturated Sphaeroforma cultures were mildly sonicated prior to diluting them 250X in fresh marine broth to initiate a synchronized culture. To assess for any temperature dependency, cultures were grown at both 17˚C and 12˚C, and measurements were conducted every 12 hours for a duration of 72 hours. For every time point, cells were fixed using 4% formaldehyde and 250 mM sorbitol for 30 minutes before being washed twice with PBS. For nuclei staining, cells were centrifuged at 1,000 rpm for 3 minutes after fixation and washed again 3 times with PBS before adding DAPI at a final concentration of 5 μg/mL to 5 μl of concentrated sample. DAPI-stained samples were imaged to measure DNA content and coenocyte size. It is important to note that results obtained from fast-settling mutants prior to cell release correspond to measurements of unicellular coenocytes (24 hours for 17˚C and 36 hours for 12˚C), whereas results collected after cell release correspond to measurements of multicelled clumps (48 hours for 17˚C and 72 hours for 12˚C). For cell wall staining, cells were incubated with Calcofluor-white (Sigma-Aldrich) at a final concentration of 5 μg/ml from a 200X stock solution prior to fixation. Cells were then fixed as previously mentioned and concentrated before being disposed between slide and coverslip.

Live cell imaging
For live cell imaging, saturated cultures were diluted 250× in fresh MB medium inside a μ-Slide 4-or 8-well slide (Ibidi) at time 0. To ensure oxygenation during the whole period of the experiment, the cover was removed. To maintain the temperature at 17 or 12˚C, we used a P-Lab Tek (Pecon GmbH) Heating/Cooling system connected to a Lauda Ecoline E100 circulating water bath. To reduce light toxicity, we used a 495-nm Long Pass Filter (FGL495M-ThorLabs). For plasma membrane live staining (Fig 3B, S3 Video)

Image analysis
Image analysis was done using ImageJ software (version 1.52) [93]. For nuclear content distribution across Sphaeroforma sp.'s life cycle, fixed and DAPI-stained coenocytes were imaged, and the number of nuclei per coenocyte was counted using the ObjectJ plugin in imageJ. To compute nuclear duplication times, log2 of geometric mean of DNA content was calculated as: Note that for S. tapetis, nuclear doubling times could not be computed due to the asynchrony in growth. For measurements of cell volume in live and fixed cells, we used the oval selection tool to draw the contour of each cell and measured cell perimeter. As cells are spherical, we computed cell volume as V = 4/3πr 3 where r is the cell radius. For measurements of clumps perimeter, we transformed the images into binaries to ensure later segmentation. We then used the particle analysis function in ImageJ with a circularity parameter set to 0.15 to 1 to measure cell perimeter. For nuclear number-to-volume ratios, the number of nuclei was divided by the coenocyte volume measured as previously described for fixed cells. All figures were assembled with Illustrator CC 2020 (Adobe). Several figures were generated using ggplot2 in R version 4.0.5 [94].

Sphaeroforma arctica genome sequencing and assembly
Genomic DNA was extracted for the ancestral strain (AN) and a single clonal isolate from each evolved population (S1-S10) using QIAamp DNA Blood Midi Kit (Qiagen) following the manufacturer's recommendations from 50 mL culture incubated at 17˚C for 5 days in 75 cm 2 flasks. The Qubit (Invitrogen) quantification ranged between 3 and 13 μg of genomic DNA in total. All of the subsequent steps were performed by the CRG Genomics Unit (Barcelona): Sequencing libraries were prepared from the pure high molecular weight DNA using TruSeq DNA HT Library Preparation kit (Illumina HiSeq Sequencing v4 Chemistry). A paired-end library with a target insert size of approximately 500 bp was sequenced on an Illumina HiSeq2500 platform in paired-end mode, with read lengths of 125 bp. The resulting paired raw read files were demultiplexed by the sequencing facility and data stored in 2 separate, gzipcompressed FASTQ files of equal sizes. Genome sequencing data have been deposited in NCBI SRA under the BioProject accession PRJNA693121.

Bioinformatic analyses of the genomes
Data processing. On average, each paired-end sequencing library contained approximately 58.8 million reads of 125 bp sequence lengths (S2 Table), equaling approximately 7.35 billion base pairs (Gbp). From these data, we carefully removed adapter sequences and reads shorter than 50 bp from the raw read data using trimmomatic v0.36 [95], yielding an average of approximately 4.36 Gbp of filtered sequence data per genome (S2 Table). The quality of both raw and trimmed sequencing data was assessed in FastQC v0.11.7 [96]. In more detail, based on the FastQC output for raw reads, we initiated trimming by calling the following parameters: ILLUMINACLIP This translates into the following trimming steps: • cut adapters and other Illumina-specific sequences using a combined file of default adapters ("Nextera.fa" AND "TruSeq3-PE2.fa") to catch as many spurious contaminations during library prep as possible, with seed mismatches = 3; palindrome clip threshold = 25; simple clip threshold = 10; • end-clipping of the final 15 bp of all reads, due to evidence for elevated adapter content; • quality-clip all bases on leading ends as long as bases were of lower quality than Q < 30; • removing all bases on trailing ends as long as bases were of lower quality than Q < 25; and, • finally, conducting a sliding window approach, where the reads were trimmed once the average quality within a window of 4 consecutive bases falls below a threshold of Q <28.
The trimming output consisted of 4 FASTQ files for each genome, of which 2 files contained intact paired-end reads, and another 2 files containing all unpaired reads for each end separately.
Repeat-masking of reference genome. We relied on the latest (i.e., fourth) assembly version of the S. arctica reference genome (Sarc4; [32]) for variant detection and annotation. Initial runs of the read-mapping steps revealed a high proportion of tightly clustered variants (approximately 20.4%) concentrated in certain regions of the genome (S3 Table). We investigated this phenomenon and determined that the issue was caused by repetitive stretches of sequence and thus decided to mask all of the potentially problematic repeat regions. For this, repetitive regions were screened for and properly annotated in RepeatMasker v4.1.1 [97] relying on the 20181026 release of GIRI RepBase database for annotations [98] and applying the slow high-sensitivity search mode with the following parameters: -pa 10 -s -gff -excln -species Opisthokonta.
• Variant calling. We called variants using CLC's "Fixed Ploidy Variant Detection" assuming a haploid genome (Ploidy = 1) for S. arctica. We further required a variant probability of at least 90%, ignoring positions in excess of 90× coverage, broken read pairs and nonspecific matches. All variants needed to be covered by at least 10 variant-bearing reads and a minimum consensus of 80% (i.e., at least 10/12 variant supporting reads). We also applied the following read quality filters: neighborhood radius = 10; minimum central quality = 20; minimum neighborhood quality = 20), and read direction filters (direction frequency = 0.05; relative read direction filter = yes (significance = 0.05); read position filter = yes (significance = 0.05).
• Variant filtering. (1) To consider only variants that have emerged through the course of evolution, we automatically removed mutations present in the ancestor (AN). (2) We then went on to manually curate all mutation predictions. We did this by aligning the mapping tracks (profiles) of all resequenced genomes (AN, S1-S10) and screening all initial candidate mutations. The overwhelming majority of variant calls were visually shared across all evolved clones but not called universally due to extremely low statistical support, or low local coverage in some of the genomes. The final dataset hence only contained 25 predicted variants (one shared among 2 lineages, hence, representing 26 independent mutations) and exported these for each of the 10 evolved single clones as Variant Call File (VCF) format.
• Variant annotation. We annotated filtered variants from converted VCF files with breseq v0.33.2 [99] for convenience. More specifically, we converted VCF files into breseq's Genome Diff (GD) file format using the command "gdtools VCF2GD". We then annotated all genomes in GD format jointly by running the command "gdtools ANNOTATE" and specifying original S. arctica genome assembly (Sarc4) in GenBankFormat (GBK) as reference. Mutations were tabulated and sorted. Finally, we counted and categorized both observed mutations and nonsynonymous and synonymous sites at risk across the reference genome using the command "gdtools COUNT -b" for statistics and the calculation of the compound dN/dS ratio.
Supporting information S1 Data. Our data used for estimating sedimentation rate in metric units and density in mass per volume. "V_au" is sedimentation rate in AU OD600 units per second. "R_meters" is the mean radius of cells for the genotype and environment. "Species," genotype, temp, and hours_growth are information for aggregating the data based on genotype and environment. "V_mu" is the maximum likelihood estimate of cellular sedimentation rate-the mean of all out-of-bag samples used in fitting-in meters per second. "V_se" is the standard deviation of cellular sedimentation rate across all out-of-bag samples in the fitting. "pp_mu" the maximum likelihood estimate of density-the mean of all out-of-bag samples used in fitting-in kg/m 3 . "pp_se" the maximum likelihood estimate of density-the standard deviation of all out-of-bag samples used in fitting-in kg/m 3 .
(XLSX) S2 Data. Smayda dataset used for calibrating our dataset. Adapted from [37] Appendix Table 1. "classification" species are from his table annotation. "salt_percent" is the percent salt reported in the table. "V_meters_per_second" is the velocity reported, converted from meters per day to meters per second. "Pf" is the density of the media based on the reported temperature and salt concentration, estimated by the data from Millero and Huang. "R_meters" is the mean radius in meters. When a range was given, this is the average of that range. "Pp_" is the density of the sample in kg/m 3 . "Pp" is the excess density kg/m 3 . "V_exp" is the expected velocity in the seawater used in our experiments (37.4 g/l). (XLSX) S3 Data. The Millero and Huang data [92] used for estimating seawater density based on salinity and temperature. "temp_C" temperature in degrees Celsius. "salinity" salt in g/l. "density_kgm3" excess density in kg/m 3 . "density" density in kg/m 3 . (XLSX) S4 Data. Fitness experiment data used for heritability measurements. "batch" is an arbitrary alphabetical letter referring to the date the replicate was performed (either D or E). "rep" is the replicate within batch/day (rep1 or rep2). "geno" is the genotype information (which genotype was being competed against the "AN" ancestor, or "alone" if alone in competition). "time" is the time the sample was collected in units' hours. "sel_cond" is the condition of the selectioneither the initial population ("init"), no selection for sedimentation rate ("NO_SEL"), selection for sedimentation rate ("YES_SEL"). "NOT_CLUMPY" is the count of microcolonies with not clumpy phenotype. "YES_CLUMPY" is the count of microcolonies with clumpy phenotype. "cell_density" is the measured cell density at time of dilution after period of growth. "N" is the total population growth calculated from cell_density per fold dilutions (1/250 per day). "gen" is the number of generations of growth based on cellular counts. Note for both "N" and "gen" that due to the synchronized coenocytic growth the cultures, they underwent a large increase in cell count after 48 hours. "date" is the date the replicate was performed (April or May 2019). "temp" is the temperature of the competition. The data underlying this figure may be found in  (TIF) S1 Video. Video of Sphaeroforma arctica AN and fast-settling (S1, S4, and S9) cultures sedimenting. Time interval between frames is 0.5 second. The movie is played at 7 frames per second (fps). We can observe rapid cell sedimentation in S1, S4, and S9 when compared to AN. Note that S1 clumps are bigger and sediment faster than the 2 other mutants. The movie was acquired for cultures pregrown for 72 hours at 12˚C and obtained with a mobile phone (Samsung A20).