Nutrient-driven genome evolution revealed by comparative genomics of chrysomonad flagellates

Phototrophic eukaryotes have evolved mainly by the primary or secondary uptake of photosynthetic organisms. A return to heterotrophy occurred multiple times in various protistan groups such as Chrysophyceae, despite the expected advantage of autotrophy. It is assumed that the evolutionary shift to mixotrophy and further to heterotrophy is triggered by a differential importance of nutrient and carbon limitation. We sequenced the genomes of 16 chrysophyte strains and compared them in terms of size, function, and sequence characteristics in relation to photo-, mixo- and heterotrophic nutrition. All strains were sequenced with Illumina and partly with PacBio. Heterotrophic taxa have reduced genomes and a higher GC content of up to 59% as compared to phototrophic taxa. Heterotrophs have a large pan genome, but a small core genome, indicating a differential specialization of the distinct lineages. The pan genome of mixotrophs and heterotrophs taken together but not the pan genome of the mixotrophs alone covers the complete functionality of the phototrophic strains indicating a random reduction of genes. The observed ploidy ranges from di- to tetraploidy and was found to be independent of taxonomy or trophic mode. Our results substantiate an evolution driven by nutrient and carbon limitation.

C ompetition for nutrients and organic carbon are among the major ecological selective forces in the evolution of eukaryotes. The presence of phototrophic, heterotrophic, and mixotrophic taxa and evidence for multiple gains and losses of photosynthesis across nearly all major eukaryotic supergroups reflects the significance of nutritional constraints in the evolution of life [1][2][3] . Taxa that have experienced such an evolutionary modification of their basic nutritional strategy must bear hallmarks of this fundamental switch in their genomes. Here we track down the genomic fingerprints of the eco-evolutionary constraints linked to the varying significance of nutrient and carbon shortage.
In this paper, we refer to obligate heterotrophy and phototrophy if not stated otherwise. Recent hypotheses suggest that the shift from photo-to mixotrophy was one way to overcome nutrient limitations while the shift towards heterotrophy was caused by carbon limitations 4,5 . In contrast to former individual case studies, here we use the parallel evolution of heterotrophic Chrysophyceae from phototrophic and mixotrophic ancestors in order to separate general directions and constraints in genome evolution from those related to the nutritional mode.
For phototrophic organisms, nutrients usually are a limiting factor. Phagotrophic uptake of bacteria and with that of nutrients and organic compounds could overcome limitations of essential nutrients such as nitrogen and phosphorus 6 . Alleviating the nutrient limitation may drive taxa into another challenge as nutrient shortage may simply be replaced by prey/carbon shortage. Phototrophs benefit from large cell sizes increasing the surface to maximize sun radiation, whereas heterotrophs profit from small cells, which enhance prey effectivity and enable feeding on ultramicrobacteria 4,5 . Genome sizes decrease from phototrophic to heterotrophic chrysophytes 5 . This reduction minimizes the cost of nitrogen and phosphorus by saving nucleotides. Furthermore, the nucleotides thymine and cytosine differ by one nitrogen, making AT pairs cheaper than GC. Consequently, a low GC content could be caused by nutrient limitation. The changing relevance of either nutrient or carbon limitation alters constraints in genome evolution which should be reflected by the incorporation of nucleotides with different costs, the loss of obsolete genes, and the evolution of new gene functions for a predatory lifestyle.
The Chrysophyceae within the Ochrophyta (Stramenopiles) are especially suited to investigate the evolutionary significance of this shift of the nutritional mode since the loss of photosynthesis occurred several times independently within this group [7][8][9][10] . Using transcriptome sequencing, Graupner et al. 11 analyzed plastid-targeting and -encoded genes of heterotrophic chrysophytes compared to photo-and mixotrophs in which they identified different stages of plastid pathway reduction and degradation of accompanying structures. Dorrell et al. 10 extended the analysis by including plastid genomes and identified shared losses of function across chrysophytes. To our knowledge, all existing studies so far focused on the plastid and plastid-related functions but none yet on changes in the nuclear genome related to changes in trophic strategies. Findings from transcriptome sequencing of 18 chrysophyte strains provide first insights into molecular changes between trophic modes showing that heterotrophs possess a reduced repertoire of genes related to photosynthesis but an increased or upregulated repertoire of pathways associated with food uptake and motility 12 . Apart from this, changes in the nuclear gene content, genome size, GC content, ploidy, and further genomic features connected to trophy have not been analyzed. Therefore, in this study, we examine genomes of 16 chrysophytes including phototrophic, mixotrophic, and heterotrophic lineages to investigate the impact of the nutritional shift and its drivers.
We hypothesize that: (1) The nuclear genome is reduced in size from photo-over mixo-to heterotrophic species, either as a result of nutrient limitations or as a size adaptation to feed on small bacteria. Likewise, the ploidy in heterotrophs is lower. (2) Accompanying the genome reduction and specialization of heterotrophs, we will find a loss of functional genes, especially no longer required genes for photosynthesis, biosynthesis of certain amino acids, and plastid-targeting genes, as well as a decrease in intergenic regions, reflected by an increased gene density. (3) The GC content of photo-and mixotrophic species is lower than in heterotrophic species if nutrients were the limiting factor during evolution.

Results
Genome assembly and gene content. Sequencing generated at least 60 million 150 bp-long paired-end Illumina reads for each strain (see Table 1) and a total of 700,000 PacBio reads (details see Table S1). The comparison of binning strategies has shown that MetaBAT removed much more reads than MaxBin2 resulting in one-third smaller assembly sizes (details see Table S2).
Since the steps for filtering prokaryotic reads were identical resulting in similar qualities, we chose MaxBin2 for binning. On average 61.2% (median) of the sequencing reads were used for the assembly (see Table 1). Despite similar assembly sizes and N50 values between axenically and non-axenically cultured strains (see Table 1, 2), the BUSCO analysis revealed partly incomplete genomes. The completeness of recovered genes in the draft genomes was on average 59% (see Table 3). Assembly size and with that detection of genes and completeness of the genomes seems to be affected by the presence of bacteria; i.e., the quality of the assemblies was higher for axenic strains. However, the following analyses of ploidy, GC content, and gene density are independent of genome completeness. We predicted 23,000 to 120,000 genes for the analyzed chrysophyte species (see Table 4). In assembled genomes, parts of long genes can possibly lie on several contigs and be counted multiple times thus overestimating the number of genes. We, therefore, clustered similar or duplicated genes to orthologous groups (OG) obtaining numbers in the range from 17,000 to 60,000, comparable to gene numbers obtained from transcriptome sequences of Chrysophyceae (8275 to 72,269 12 , Ochromonas sp.: 19,692 genes 13 ). Between 14 and 25% of these orthologous groups could be annotated with a KEGG Orthology (KO) ID for the assignment to the metabolic pathway. It is not surprising that additional use of PacBio sequencing and axenic strains lead to better assembly qualities (see supp. Table S3, Fig. S1). However, the omission of axenic and PacBio sequenced species leads to similar results concerning gene and pathway composition for the different nutritional modes. The pathway compositions regarding the KEGG functional groups differed mostly between single strains, instead of nutritional modes (see Fig. 1). Investigated pathways were almost complete with more genes present in the genome compared to the transcriptome (see Fig. 2, also Figs. S7 to S12 and Supplementary Data 1).
Pan genome and core genome. To expand the analysis of genes and pathways in single species, we compared the sets of orthologous groups (OGs) between the trophic modes. The pan genome, the entire gene set of a group of species, as contrasted to the core genome, the intersection of genes within the group of species. The pan genome of all investigated phototrophic species contained 47,750 orthologous groups, while the pan genome of     Genome size. Considering only genomes with high completeness, hetero-and mixotrophic species did not differ largely in assembly size, ranging at around 67 Mbp. We observed significant differences in assembly size between all three groups (ANOVA p-value < 0.05), but these differences were mainly due to the substantially larger genomes of the phototrophic species (see Table 2). Due to the few species in this group, we refrain from drawing a definite conclusion here.
Gene density. In general, there is a correlation between genome size and gene density in eukaryotes 14 . We expected to find a correlation between nutritional mode and gene density. Especially, heterotrophs should have a high gene density, because of the strong selection towards smaller cells. The average gene density ranged from 1031 (phototroph), over 1049 (mixotroph) to 1520 (heterotroph) genes/Mb. A slight trend towards higher gene density with advanced heterotrophy was visible, but these differences were not significant (p-value = 0.32, see Fig. 4).
GC content. The mean GC content of the phototrophic (43.5%) and mixotrophic group (43.4%) was similar, whereas heterotrophs (51.6%) had a strongly increased GC content (ANOVA: p-value = 0.0041, posthoc test mixo-versus heterotrophs p-value = 0.0053, see Fig. 4). We separated the total GC content based on coding sequences, introns, and non-coding sequences (see Table S4). The mean GC content of noncoding regions was smallest (35.5%), intron GC content was in between (41.4%) and coding sequences had the highest GC content (54.8%) (pairwise t-tests, each p-value ≤ 0.001). Using only the third position of the coding sequences (GC3) resulted in a GC value of 52.6%. The coding, GC3, intron, and non-coding GC content were independent of the nutritional mode (ANOVA: p-value > 0.05). Significant dependencies were also not apparent between assembly size, or total genome size, and the GC content (p-values > 0.5).

Discussion
The comparative genome analysis of chrysophytes revealed genome reduction accompanied by an increase in the genomic GC content and an increase in gene density as a concomitant phenomenon of the evolutionary switch from phototrophy to heterotrophy. Gene reduction occurred in all heterotrophic lineages but different genes were lost in the various phylogenetic lineages presumably reflecting differential evolutionary selective pressures and adaptations to different environments in the distinct lineages. Even though genome reductions evolved independently, gene losses were predominately observed for non-annotated genes; i.e., genes of unknown functions, while genes and pathways of the primary metabolism were kept in all lineages.
We found considerable differences in genome size and ploidy between strains. While genome size reduction seemed to be related to the change of trophy but independent of phylogeny, a variation of ploidy was found in different lineages and seemingly independent of the nutritional mode. From the assembly sizes, we could not definitively assess differences attributable to the different trophic modes, but the two larger phototrophic genomes suggest it (see Table 2). Nucleic acids are amongst the most phosphorus and nitrogen-demanding biomolecules and large genomes are costly to build and maintain, thus under nutrient limitation a reduction in genome size is expected. Further, the genome size has been shown to correlate with cell volume (e.g., ref. 15 ), which is reduced in phagotrophic species preying on ultramicrobacteria. This is supported by a study from Olefeld et al. 5 analyzing 46 chrysophyte strains, in which a reduction in cell volume and genome size from photo-over mixo-to heterotrophs could be clearly shown with flow cytometry. In this study 5 , the size of the phototrophic strains was considerably larger, although it should be noted that this applies to algae with secondary plastids, while species with primary plastid may be much smaller e.g., Ostreococcus tauri 16 .
The study from Olefeld et al. 5 calculated haploid genome size based on the assumption of diploid genomes in chrysophytes. As we show here, ploidy varies between strains and thus variation in ploidy overlays the general trend of genome size reduction. Taking the different ploidies into account the genome size estimates based on flow cytometry 5 and based on assembled genomes (this study) largely correspond except for few strains for which the ploidy level could not be unambiguously determined ( Table 2). In particular, the genomes of phototrophic strains were considerably larger as compared to those of heterotrophic and mixotrophic strains. Since the two methods do not match consistently, we rather rely on the genome size results from flow cytometry. We compared gene density and show statistical results for both methods in the supplement (see supp. Table S5, Fig. S13). In no case, we could find a significant difference between the three trophic groups concerning gene density.
Ploidy varied between diploidy and tetraploidy and was independent of phylogeny and nutritional mode. The independence from phylogeny was not surprising as different ploidy levels have already been demonstrated within one species 17,18 . Apparently, a change in ploidy occurs frequently and is often associated with broader ecological niches and/or invasiveness due to increased heterozygosity and flexibility 19 . In particular, the extra gene copies in polyploid genomes may be beneficial for evolving new functions, novel gene combinations, and modified gene expression. For instance, flowering plant lineages polyploidize at a rate that is about 2-10% of the speciation rate 20,21 . In protists, which replicate predominantly asexual, the rate could even be higher, since the disadvantages of vulnerability to infertility or pairing difficulties in meiosis 22 are of little relevance. Furthermore, polyploidy is advantageous with respect to gene redundancy 22 , especially in predominantly asexual species in order to prevent the accumulation of mutations (Muller's ratchet) 23 . However, not only increases in ploidy but also a decrease of the ploidy level may frequently occur as it has been demonstrated for tetraploid or triploid yeast within a short time period of only 186 generations 24 . Irrespective of the potential advantages of polyploid genomes, an increase in ploidy level counteracts the attempt of reducing genome size since cell size and DNA content usually correlate 25,26 . Polyploid organisms form larger cells 27,28 , but the reason for this correlation is not yet clearly understood. Both, the higher amount of DNA and the higher expression of proteins have been suggested to be causative 29,30 , however, deviations thereof are known even within the same genus 31 . These deviations imply that more factors are decisive for the cell size. Likewise, even though plants often increase their cell size by increasing their ploidy level 28 , preventing the ploidy enhancement did not affect cell size 32 . Summarizing, our data point to a trade-off between the advantages of genome size reduction and polyploidization: Cell size reduction-and with that genome size reduction-is discussed as an adaptation to increase the efficiency in preying on small bacteria, in particular on ultramicrobacteria. At the same time, genome size reduction lowers the costs of DNA buildup and reproduction. In contrast, polyploidization may allow for higher flexibility enabling the taxa to populate broader niches and may prevent or at least slow down the accumulation of mutations. When nutrients are not limiting this trade-off should level off as it has also been shown for cultures with repeated transfer 17 , since selection pressure decrease under laboratory conditions. However, due to the selective advantages of both, genome size reduction and polyploidy, it seems reasonable that polyploidy may also evolve under nutrient limitations despite the extra cost for genome maintenance.
To analyze the gene composition and gene loss we compiled for each nutritional group all combinations of the pan genome and the core genome. The pan genome comprises genes present in at least one representative, while the core genome comprises genes present in all representatives. One major finding of our study is the deviation between the core and pan genome of the heterotrophic taxa. The small core genome in contrast to the large pan genome of the heterotrophs implies that the strong genome reduction in the evolution of the heterotrophs was accompanied by a strong niche specialization. Even though numerous genes have been lost in each heterotrophic lineage the distinct lineages lost different genes. Phylogenomic patterns were consequently only found in closely related lineages, i.e., loss of similar genes and pathways (in particular in the clade comprising Poteriospumella, Poterioochromonas, and Chlorochromonas). Further, our data demonstrate that in all trophic groups some of the genes are not found in at least one of the other trophic groups. However, combining the group of mixotrophic and heterotrophic species, their genes cover almost all genes and pathways found in the phototrophic species. We expected to see a gradual reduction of pathways related to photoautotrophy depending on the stage of nutritional mode. This gene loss pattern has been shown for the plastome, even if exceptions such as genes related to Rubisco could still persist 33 . The study of Graupner et al. 11 suggests that for pathways directly linked to photosynthesis and genes encoded in the plastid a certain sequence of gene losses occurs, accompanying the shift from mixotrophy to heterotrophy. Here we show that this is not the case for the majority of nucleus-encoded genes but gene loss seems to be different for distinct lineages indicating that different strains presumably were exposed to different evolutionary forces. This is supported by an increased fraction of genes related to environmental information processing and secondary metabolism in the photo-and mixotrophic group (see Fig. 1) while the heterotrophic species have a more specialized and unique gene inventory in these functional groups.
We further demonstrate that primary pathways and annotated pathways of the secondary metabolism are almost complete in the genome. Gene reductions concern mostly genes that could not be annotated; i.e., genes of unknown function which presumably are important for niche specialization but not for the basic requirements of the cells. For instance, the pathways for the biosynthesis of amino acids were complete in all investigated strains even though amino acids can, in principle, also be taken up with food in the phagotrophic taxa (see Fig. 2, also Figs. S7 to S12). This completeness of pathways seems to be different for plastidencoded genes as suggested by transcriptome data and supported by the structural reduction as observed in microscopical analyses 9 . The maintenance of the pathways of the primary (and to a large part also of the secondary) metabolism could be advantageous for backing-up nutritional requirements in the case of food shortage. As we analyzed the occurrence of genes, our results do not necessarily imply that all genes are functional. But the fact that pathways were mostly found to be complete supports this assumption. Still, undetected mutations could have rendered genes non-functional or modified their function, for certainty, additional verifications would be needed.
There is a small trend towards higher gene density with advanced heterotrophy, but possibly due to the small sample size, the difference between groups was not significant. The gene densities were generally high compared to gene densities of other organisms. Stramenopiles have on average a gene density of 200-400 genes/ Mb 18 , whereas prokaryotes or archaea typically have around 1000 genes/Mb. As eukaryotes with small genomes are known to have gene densities similar to bacteria 14 and constant gene densities independent of genome size 34 our estimates still seem reliable. However, this high gene density could partly be due to an overestimation in the gene prediction, in particular, due to the merging of repeat regions during assembly and the potential presence of genes on both strands. Since the gene density between strains sequenced only with Illumina is comparable to strains additionally sequenced with Pacbio, a bias due to collapsed repeat regions should be very small.
The GC content is known to be correlated with genome size in bacteria but the relationship for other kingdoms is less clear. Studies analyzing genome size and base composition in available sequenced genomes in various kingdoms found a correlation between average GC content and genome size which was positive in bacteria (specifically Proteobacteria and Actinobacteria), weakly positive in Ascomycota fungi and some plants, but negative in animals and indifferent in two analyzed protistan phyla [35][36][37] . We show that the GC content of chrysophytes did not correlate with genome size. In contrast, the total GC content was significantly different between nutritional modes. The total GC content differs due to several factors, which were not significant individually (e.g., gene density, the difference between nutritional modes in GC content in introns, coding sequences, and non-coding sequences), but contributed jointly to the measured difference. While the mean GC content of the phototrophic and mixotrophic groups was similar, the heterotrophs showed a strongly increased GC content. In general, an increased GC content is associated with higher needs regarding nutrients and energy 38,39 . Adenine and guanine require 5, cytosine 3, and thymine/uracil only 2 nitrogen atoms, thus the difference in costs is observable in the GC content of a genome. Species living in or adapted to nitrogen-limited environments therefore often use nucleotides that require fewer nitrogen such as A and T 40 . In addition, at least in bacteria, organisms show distinct GC patterns that are not explained by phylogeny but by similar environments 41 . It has been shown that plants that require more nitrogen to conduct photosynthesis experience stronger selection to minimize nitrogen biosynthesis costs 42 . As the evolution towards heterotrophy; i.e., the establishment of an alternative nutrient source, is typically related to adaptations to low nutrient availability, a lower GC content should be expected in phototrophic chrysophytes. The evolutionary constraint of nutrient limitation should be already relaxed in mixotrophic species, but as the photosynthetic apparatus is costly and accounts for about half of the cells proteins 43 , mixotrophs may be subject to similar constraints as the phototrophs. Heterotrophic species should predominantly be limited by carbon, not by nutrients, and thus can afford a higher GC content. The increased energy demand of high GC content led to a low GC content in non-coding regions and introns, whereas coding sequences remained GC-rich encoding cheaper amino acids 38 .

Conclusions
We demonstrate that nutritional constraints drive genome evolution and that limitations in nutrient and carbon acquisition leave footprints in genome evolution. Based on comparative genome analysis of lineages that independently and in parallel evolved heterotrophy from mixotrophic ancestors we separated genomic footprints linked to a nutritional switch from random shifts. In particular, genome size reduction and a shift in GC content reflect changes in nutrient limitation during the evolution of obligate heterotrophy from phototrophic and mixotrophic ancestors. Gene losses accompany these changes but vary between lineages indicating a presumably increasing niche separation and differentiation between the distinct heterotrophic lineages. In the broader context of eukaryotic mega evolution, the discovery of disparate gene losses in different lineages and its implication of an accelerated differentiation is intriguing as it adds a new facette on the evolution of eukaryotic diversity. Nutrient and carbon limitation may not only be crucial for the transformation of feeding strategies but furthermore speed-up the evolutionary diversification and thus contribute to rapid radiations on short evolutionary time scales.

Methods
Cultivation and sequencing. We cultivated and sequenced 16 strains (see Table 5) according to Hahn et al. 44 and Majda et al. 18 . Non-axenic heterotrophic and mixotrophic cultures were grown with bacterial food supply (Limnohabitans planktonicus; strain IID5T). Two days before DNA harvesting in these cultures the feeding with bacteria was omitted.
Genome assembly and binning. Unless otherwise stated, the default settings were used for the following programs. An automated workflow with Snakemake 45 processed the sequencing data. Supplementary Fig. S2 gives an overview of the assembly and binning procedure. Reads were quality checked, filtered and adapters were removed by the sequencing provider, further, the read quality was checked by FastQC (v0.11.5; with a Phred quality score criteria of around 10 for PacBio reads and above 20 for Illumina reads 46 ). SPAdes (v3.13.0; with parameters: -meta 47 ) assembled the Illumina reads. If PacBio reads were available they were incorporated in the SPAdes assembly and additionally assembled with CANU (v1.8; with parameters: genomeSize = 100 m, correctedErrorRate = 0.105 48 ), whereby the genome size estimations from Olefeld et al. 5 were used if possible. We renounced from genome size estimation based on k-mers, since it often led to size bias in chrysophytes 18 , especially through read filtering in non-axenic strains. To classify the reads we combined a tool using compositional features and other applying taxonomical methods. The tool MaxBin2 (v2.2.5 49 ) binned the contigs created by SPAdes. Subsequently, Kraken2 (v2.0.7 50 ) classified the bins taxonomically. Therefore, the NCBI database (Release 2018-11-01) and the assemblies of the axenic cultures Chlorochromonas danica, Poterioochromonas malhamensis, and Poteriospumella lacustris were used as a reference. Including the axenic cultures enabled to classify chloroplast bins as eukaryotic instead of prokaryotic and prevented their exclusion. Bins consist out of several contigs. In some cases, the bins contain contigs classified as eukaryotic, as well as prokaryotic or unclassified. Ambiguous cases are very rare (less than 5%), but we have nevertheless defined a criterion for this. Depending on the bin composition, we classified and processed unclassified contigs as follows: 1. If the number of eukaryotes was at least half the number of bacteria, the unclassified were defined as eukaryotic. 2. If the number of unclassified contigs was twice as large as the number of bacteria, the unclassified were defined as eukaryotic. We have defined these points because the number of known bacteria in the NCBI database is much higher than that of related protists.
In addition, the binning tool MetaBAT(v2.12.1 51 ) was tested followed by the same classification and filtering steps. The bins were merged into two files containing either eukaryotic or bacterial sequences. The Illumina reads were aligned with Bowtie2 (v.2.3.0 52 ) against the bacterial contigs. Subsequently, the hits were mapped again with Bowtie2 against the eukaryotic contigs. Unaligned reads were marked as bacterial and excluded from further processing. Normally, reads mapping to both the bacterial and the eukaryotic contigs were randomly assigned. To keep the eukaryotic reads in this case we used incremental mapping.
The contigs of the CANU assembly were only binned by Kraken2, as the long contig length caused sufficient classification accuracy. In a second step, we filtered out all contigs that were not eukaryotic. We repeated the SPAdes assembly (without -meta parameter) with filtered Illumina reads and filtered PacBio reads if available (parameter: -untrusted-contigs). Contigs smaller than 500 bp were discarded.
Gene annotation and clustering. Genes were annotated according to ref. 18  OrthoFinder (v2.2.6; with parameters: -S diamond) clustered the protein sequences of all strains to orthologous groups (OG). The majority of gene annotations within an orthologous group determined the annotation of the whole group. We built the union (all OGs) and intersection (OGs shared among all) between species of one nutritional mode (phototroph, mixotroph, and heterotroph) and between nutritional modes. In addition, for these sets, the composition of KEGG functional groups was determined.
GC content and gene density. The GC content was assessed for the whole genome and for intergenic regions. Significant differences in gene density and GC content between the nutritional modes were calculated using analysis of variance (ANOVA) in R. The gene density (d) was calculated for the number of genes (n) by the formula: Pathway analysis. Present annotated genes were visualized in the KEGG pathway maps by the R package Pathview (v1.24.0 63 ). We focused on the following pathways: biotin metabolism, histidin metabolism, tryptophan metabolism, nitrogen metabolism, lysine biosynthesis, pyhenylalanine, tyrosine, and tryptohan biosynthesis, thiamine metabolism and valine, leucine, and isoleucine biosynthesis. In addition, a comparison was made to determine whether genes from a transcriptome dataset 12 were also present.