# A multi-omic resource for exploring microbial eukaryotes in the meromictic freshwater Lake Pavin

Damien Courtine, Cécile Lepère, Ivan Wawrzyniak, Anne Moné, Hermine Billard, Jonathan Colombet, Arthur Monjot, Corinne Cruaud, Corinne da Silva, Jean-Marc Aury, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Damien Courtine, Cécile Lepère, Ivan Wawrzyniak, Anne Moné, Hermine Billard, et al. A multi-omic resource for exploring microbial eukaryotes in the meromictic freshwater Lake Pavin. Scientific Data, inPress, Metagenomics in marine and terrestrial aquatic environments, 10.1038/s41597-026-06573-0. hal-05477144

HAL Id: hal-05477144

https://hal.science/hal-05477144v1

Submitted on 26 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# A multi-Omic resource for exploring microbial eukaryotes in the meromictic freshwater Lake Pavin

Received: 8 October 2025

Accepted: 6 January 2026

Cite this article as: Courtine, D., Lepère, C., Wawrzyniak, I. *et al.* A multi-Omic resource for exploring microbial eukaryotes in the meromictic freshwater Lake Pavin. *Sci Data* (2026). https://doi.org/10.1038/s41597-026-06573-0

Damien Courtine, Cécile Lepère, Ivan Wawrzyniak, Anne Moné, Hermine Billard, Jonathan Colombet, Arthur Monjot, Corinne Cruaud, Corinne Da Silva, Jean-Marc Aury, Didier Debroas & Gisèle Bronner

We are providing an unedited version of this manuscript to give early access to its findings. Before final publication, the manuscript will undergo further editing. Please note there may be errors present which affect the content, and all legal disclaimers apply.

If this paper is publishing under a Transparent Peer Review model then Peer Review reports will publish with the final article.
# A multi-Omic resource for exploring microbial eukaryotes in the meromictic freshwater Lake Pavin.

## Authors list

Damien Courtine<fnref n="1" />*, Cécile Lepère<fnref n="1" />, Ivan Wawrzyniak<fnref n="1" />, Anne Moné<fnref n="1" />, Hermine Billard<fnref n="1" />, Jonathan Colombet<fnref n="1" />, Arthur Monjot<fnref n="1" />,<fnref n="2" />, Corinne Cruaud<fnref n="3" />, Corinne Da Silva<fnref n="4" />, Jean-Marc Aury<fnref n="4" />, Didier Debroas<fnref n="1" /> and Gisèle Bronner<fnref n="1" />.

## Abstract

Although recent advances in high-throughput sequencing have greatly expanded our understanding of microbial diversity and function in aquatic ecosystems, progress in studying freshwater microbial eukaryotes has been more limited, mainly due to their large genomes, immense diversity, and largely uncharacterised physiologies. In this work, we present a comprehensive multi-omic dataset, eukaryote-centred, including targeted-metagenomic (18S rDNA V4 and V9), metagenomic, metatranscriptomic and single amplified genomes (SAGs). Both the oxic and anoxic layers of Lake Pavin (France), a permanently stratified freshwater lake, were sampled at four distinct times throughout 2018, by day and night, targeting microbial eukaryotes of two size classes (0.65-10 µm and 10-50 µm). This dataset comprises 106 eukaryotic metagenome-assembled genomes (MAGs), over 9 million unigenes and 11 SAGs, encompassing several under-represented taxa in public databases (e.g. Perkinsea, Chytridiomycota, Cryptista). Altogether, this dataset represents a resource for exploring the functional diversity and spatio-temporal dynamics of microbial eukaryotes.

## Background and summary

As reservoirs of microbial biodiversity, aquatic ecosystems play a central role in global biogeochemical cycles. However, the specific and functional diversity of microorganisms inhabiting these ecosystems remains largely unknown<fnref n="1" /> because a significant proportion of these microbes are uncultivable using traditional laboratory methods. This is particularly true for microbial

<footnote n="1">¹Laboratoire Microorganismes: Génome et Environnement, Université Clermont Auvergne, CNRS, Clermont-Ferrand, France</footnote>

<footnote n="2">²UMR6538 GEO OCEAN, Université de Bretagne Occidentale, CNRS, IRD, Ifremer, 29280, Plouzané, France</footnote>

<footnote n="3">³Institut de Biologie François Jacob, Genoscope, CEA, Université Paris-Saclay, Evry, France</footnote>

<footnote n="4">⁴Génomique Métabolique, Genoscope, Institut François Jacob, CEA, CNRS, Univ Evry, Université Paris-Saclay, Evry, France</footnote>

<footnote>*Corresponding author, damien.courtine@uca.fr</footnote>
eukaryotes, which have complex metabolisms and whose physiologies may vary upon environmental conditions². Recent advances in high-throughput sequencing technologies and multi-omics approaches have revolutionised the study of these communities. The complementarity of these approaches, whether targeted (amplicon-based) or untargeted (metagenomics, metatranscriptomics), offers a more comprehensive view of microbial communities within an ecosystem. This includes their taxonomic composition, genetic potential, gene expression patterns, and functional activities, thereby revealing the complexity of interactions among these microorganisms.

However, non-targeted approaches are rarely used for the study of microbial eukaryotes, largely due to the complexity involved in analysing their diverse and often large genomes, as well as the challenges in interpreting their functional roles within ecosystems. A few pioneering eukaryote-focused metagenomic studies have been conducted in ocean ecosystems³,⁴, however similar comprehensive investigations in freshwater environments remain scarce⁵,⁶.

Focusing on the meromictic Lake Pavin⁷ and its microbial eukaryotic community, we provide an exhaustive -omic dataset that includes amplicon-based (18S rDNA regions V4 and V9), metagenomic and metatranscriptomic datasets, as well as metagenome-assembled genomes (MAGs) and single amplified genomes (SAGs). The Lake Pavin is located in Auvergne (France), characterised by a mixed oxygenated upper layer (top to *ca.* 60 m) and a permanent anoxic deep layer (*ca.* 60–92 m; Figure 1). Water sampling was carried out in both layers (at 9 m and 80 m) by day and night, in April, June, September, and November, and for two size classes (0.65–10 and 10–50 µm).

Targeted metagenomic (amplicon-based) data are provided as raw reads. Whole genome shotgun metagenomic data, with a total of 8,8 billion read pairs distributed in 57 samples (277 million read pairs on average; Table1), were partitioned and assembled according to size class and oxygen status (Figure 2), leading to 2.8 million contigs of size greater than 2.5 kb, representing 19.9 Gbp. The binning of these contigs produced 5,657 bins, 106 of which were identified as putative eukaryotes, and manually refined using Anvi'o⁸ to obtain MAGs. Gene predictions were used to estimate MAG completeness which is 39.9% on average (Figure 3). Among these MAGs, 52 were classified as heterotrophs, including 25 identified as confirmed or potential parasites (Chytridiomycota, Aphelida, Perkinsea and Microsporidia). Additionally 42 MAGs belonged to pigmented eukaryotes, mostly diatoms (Bacillariophyta) and Chlorophyta (Figure 4). About 10% of the MAGs were assembled from the deep anoxic layer. These MAGs were complemented by 11 SAGs with an average completion of 22.7% (maximum 59.6%, Figure 3) affiliated to Perkinsea, Chytridiomycota, Aphelida and Bacillariophyta. To assess how representative our dataset is, we mapped metagenomic reads from diverse water bodies worldwide onto the Lake Pavin MAGs and SAGs. Clustering based on Aitchison distance calculated from transformed read counts grouped the metagenomes onto three distinct clusters. Lake Pavin MAGs and SAGs are mostly detected in metagenomes associated with lentic ecosystems like lakes or reservoirs (Figure 5a, b) and 25 MAGs, including 22 from pigmented taxa, were found across a variety of ecosystems, such as rivers, ponds or estuaries with some marginal presence in coastal environments (harbour, estuary; Figure 5a, c). These findings underline the specificity of the MAGs and SAGs presented in this
study for freshwater environments, although pigmented taxa appear to inhabit a broader range of
ecosystems compared to non-pigmented taxa.

The metatranscriptomic data, spanning 48 samples with on average 247 million read pairs (Table 1), were processed into unigenes and annotated as previously described in Carradec et al.⁹ and Monjot et al.¹⁰ respectively. Similarity analysis of these data revealed differences in expression profiles between size classes, oxic-anoxic layer and sampling periods (ANOSIM, *p*<0.001), whereas not between day and night samples¹⁰. Coupled with a trait-based approach, amplicon-based and metatranscriptomic data were used to investigate the metabolic potential and genetic signatures of ecologically meaningful functional groups of microbial eukaryotes (e.g. phototrophs, heterotrophs, saprotrophs, parasites, mixoplankton)¹⁰,¹¹. These analyses revealed a remarkable diversity in the anoxic zone of Lake Pavin as well as strong seasonal variations in the oxic layer, particularly among parasitic and mixoplanktonic microorganisms. Parasites benefited from water mixing events in spring and autumn, while mixoplankton, whose metabolic versatility was confirmed by the global analysis of functional genes, took advantage of their phagotrophic abilities during nutrient limitation specific to winter.

Overall, this dataset offers a more comprehensive genetic characterisation of the microbial eukaryotic community within a lake ecosystem featuring contrasting environmental conditions. These data complement marine datasets³,⁴ and offers new resources that are still rare in databases, such as MAGs from freshwater environments and anoxic zones. Moreover, the inclusion of genomes from under-represented taxa in this environment enriches the available resources for studying prokist diversity¹². Additionally, this integrated dataset can serve as a valuable benchmark for the development and validation of future analytical tools and methodologies.

# Methods

## Sample collection

Water samples were collected at the centre of the Lake Pavin (Massif Central, France, 45°29′45″ N, 2°53′17″ E) in April, June, September, and November of the year 2018. For each sampling event, a total of 60 L of water were collected using a Von Dorn bottle from the oxygenated layer (9 m) and the permanently anoxic lower layer (80 m), once during the day and once at night. The samples were pre-filtered sequentially on nylon filter mesh (sifting fabric) with decreasing pore-size (150 µm, 50 µm and 10 µm). Microorganisms retained on the 10 µm filter, representing the large size class (10–50 µm), were resuspended in a final volume of 300 mL of sterile water. Those corresponding to the small size class (<10 µm) were concentrated using tangential flow filtration with a Spectrum KrosFlo PES 0.65 µm column (Spectrum Labs, Rancho Dominguez, CA, USA) and collected in a final volume of 1 L. Given that filtration-based sampling can selectively impact microbial cell recovery due to variations in size, morphology, or adhesion properties, we carefully minimised handling time and optimised resuspension protocols to mitigate such biases. Samples were treated with 0.01% pluronic acid and centrifuged at 10,000 g, 4°C for 30 min. The pellets were then divided for DNA and RNA extraction with the aim to sequence all samples in duplicates.
The pellets were resuspended in 250 µL of RLT buffer for DNA extraction (Qiagen, Hilden, Germany) and with the addition of β-mercaptoethanol for RNA extraction. Samples were immediately frozen in liquid nitrogen.

DNA samples were disrupted by bead beating for 30 s at 30 Hz, three times, with 0.1 g of 0.1 mm glass beads (Sigma-Aldrich, Saint-Louis, MO, USA), followed by 1 min of centrifugation at 800 g. The supernatants were adjusted to 300 µL with extraction buffer (50 mM glucose, 10 mM Ethylenediaminetetraacetic acid (EDTA), 25 mM Tris) and incubated with 150 µL of 10% sodium dodecyl sulfate (SDS) for 5 min. Then, 5 µL of proteinase K (10 mg/mL) and 5 µL of RNAse A (10 U/µL) were added. The mixture was incubated for 2 h at 37°C. Then, 80 µL of cetyltrimethylammonium bromide (CTAB) solution (10%, NaCl 0.7 M) and 100 µL of NaCl 5 M were added before the nucleic acid extraction by phenol-chloroform-isoamyl alcohol (same volume as sample, ratio 25:24:1 v/v/v) and isopropanol precipitation. DNA pellets were washed with 70% ethanol and resuspended in 30 µL of water and stored in LoBinding DNA tubes.

The total RNA extraction were processed with the RNeasy Mini kit (Qiagen) with some modifications. Samples were disrupted by bead beating for 30 s at 30 Hz, three times, and the supernatant was kept. The sample was subjected to a second bead beating round, same conditions, with 400 µl of RLT buffer with 1% 2-mercaptoethanol and 0.1 g of 0.1 mm glass beads (Sigma-Aldrich). Both supernatants were pooled. The end of the protocol was carried out according to the manufacturer's recommendations. Total RNA was eluted in 30 µL of RNAse free water.¹⁰

Several samples, particularly those from the anoxic zone, did not yield sufficient DNA/RNA for sequencing. We therefore had to pool the duplicates, and ultimately obtained 57 DNA and 48 RNA samples (Table 1).

## Libraries preparation and sequencing

The V4 region of the 18S rDNA gene was PCR amplified using the 515F (5'-GTGYCAGCMGCCGCGGTA-3')¹³ and 951R (5'-TTGGYRAATGCTTTCGC-3')¹⁴ primers and subsequently purified as presented in Chauvet *et al.*¹⁵. The V9 region was amplified using the universal eukaryotic forward primer 1389F (5'-TTGTACACACCGCCC-3')¹⁶ and reverse primer 1510R (5'-CCTTCYGCAGGTTCACCTAC-3')¹⁷. Starting from 2.5 to 10 ng of DNA, three PCR reactions were performed for each sample (98 °C for 30 s followed by 25 cycles 10 s-98°C, 30 s-53°C, 30 s-72°C and a final extension 10 min-72°C) with the Phusion High Fidelity PCR Master Mix with GC buffer (Thermo Fisher Scientific, Waltham, MA, USA). Then, both 18S rDNA amplicons (V4 and V9) were processed with the same methods: PCR products were pooled by sample, purified using 1X AMPure XP beads (Beckman Coulter, Brea, CA, USA) and quantified with a Fluoroskan Microplate Fluorometer (Thermo Fisher Scientific). A size profile analysis was performed with a LabChip GX Nucleic Acid Analyzer (Revvity, Waltham, MA, USA). A library was prepared for each amplicon, using the NEBNext DNA Modules Products (New England Biolabs, Ipswich, MA, USA) and NEXTFLEX Unique Dual Index Barcodes (Perkin Elmer, Waltham, MA, USA). Ligated products were amplified using the Kapa Hifi HotStart NGS library Amplification kit (Roche, Basel, Switzerland), followed by 1X AMPure XP (Beckman Coulter) clean up.

Metatranscriptomic libraries were prepared from 100 ng of total RNA with the NEBNext Ultra II
Directional RNA Library Prep Kit (New England Biolabs)¹¹. The metagenomic libraries were prepared using the NEBNext Ultra II kit (New England Biolabs) using 2.4 to 250 ng of DNA. All libraries were sequenced at the Genoscope (Évry, France) with different Illumina systems (Illumina, San Diego, CA, USA): the V4 and V9 amplicons were sequenced on a HiSeq 2500 (2x250 bp) and MiSeq (2x150 bp) respectively, while both metagenomic and metatranscriptomic libraries were sequenced on a NovaSeq 6000 (2x150 bp).

## Meta-omic reads assembly and binning

Metatranscriptomic reads were assembled to form unigenes as described in Carradec et al.⁹, and were further functionally and taxonomically annotated as described in Monjot et al.¹¹. Paired-end reads from 48 samples (Table 1) were assembled on a per sample basis using Velvet¹⁸ v1.2.07 (k-mer size of 89). Contigs smaller than 150 bp were removed, and isoforms were detected using Oases¹⁹ v0.2.08. Contig redundancy was reduced by clustering with CD-Hit-EST v4.6.1 (-id 95 -aS 90), retaining the longest contig as the unigene reference sequence. This unigene catalogue was analysed with TransDecoder²⁰ v5.5.0 to identify putative proteins of a minimum length of 70 amino acids. Predicted proteins were screened against the AntiFam²¹ database v7.0 using hmm-search²² with the --cut_ga parameter and functionally annotated according to KEGG using KoFam-Scan²³ v1.3.0 (HMM profiles release 2022-01-03, e-value < 1e-5) and Pfam-A²⁴ v35.0 using hmm-search and the --cut_ga parameter. The GO terms were extracted from the Pfam annotation using the Gene Ontology cross-reference classification system²⁵. Taxonomic affiliations were inferred by comparing predicted proteins to the database associated to MetaEuk²⁶ using MMseqs²²⁷ (v407b315, mmseqs taxonomy --tax-lineage 1 --lca-mode 2 --max-seqs 100 -e 0.00001 -s 6 --max-accept 100). Proteins and associated unigenes affiliated to Human, Bacteria, Viruses and Metazoans were excluded from the catalogue. Starting from 7.9 billion read-pairs, this pipeline produced in 9,926,643 unigenes and 9,182,313 proteins, of which 30.7 % were assigned a KO annotation, and 30.8% were affiliated to the eukaryotic domain. Overall, 1,280,290 unigenes have both taxonomic and functional annotation. Expression profiles of unigenes were determined by mapping reads from each sample to the unigenes reference sequences using bwa²⁸ v0.7.15 (bwa aln -l 30 -O 11 -R 1; bwa sampe -a 20000 -n1 -N) followed by duplicate removal (samtools rmdup). Read counts were then restricted to alignments with at least 80% of the read length aligned and a minimum of 95% sequence identity.

For metagenome reads, the 57 samples (8.8 billion read pairs; Table 1) were first subjected to reads composition analysis using SimkaMin²⁹ on 50M reads and 10M k-mers per sample. The PCoA on a Bray-Curtis distance of the samples reads composition highlighted sample clustering by size class and oxygen status (Figure 2). Grouping samples according to these parameters, 4 co-assemblies were realised with MEGAHIT³⁰ v1.2.9 and parameters --k-list 3, 47, 59, 79, 99, 119, 141 --min-contig-len 1000. Contigs smaller than 2.5 kb were removed and the assemblies were then processed using the Anvi'o⁸ framework v8 for gene prediction with Prodigal³¹ v2.6.3. The short-reads were then mapped on the contigs with bwa mem³² v0.7.17, results were processed with Samtools³³ to generate BAM files. Alignments were filtered by msamtools³⁴ v1.1.0 to keep reads with more than 95% identity over a minimum of 75 nt. Anvi'o facilities were used to combine all
BAM files within each co-assembly to estimate per sample contigs abundances and run MetaBAT<sup>235</sup> with default parameters for binning. Candidate eukaryotic bins were identified by the presence of 83 single-copy core genes<sup>36</sup> (anvi-estimate-genome-completeness) and were subjected to manual inspection via the Anvi'o interactive interface (command anvi-refine) to remove spurious contigs, based on either “sequence composition” divergence or “differential coverage”. Bins completeness and redundancy were estimated with Anvi'o (anvi-summarize). This process resulted in 106 eukaryotic candidate MAGs of which 103 were retained for advanced gene prediction. As suggested by Schoenle *et al.*<sup>12</sup>, the 5,551 unrefined bins may correspond to organisms that were in association with microbial eukaryotes (e.g. symbiosis, parasitism, predation). It could also be eukaryote genomes with lower completion, that were not detected because of an insufficient coverage of the eukaryote marker gene set during the analysis.

## MAGs gene prediction

As eukaryotic gene prediction remains a challenging task, all current gene prediction tools use intrinsic (e.g. RNA-seq data) or extrinsic (e.g. protein of close organisms) evidences. In order to provide an accurate gene prediction for each MAGs, we searched public databases for proteins from closely related organisms. For this, a MMseqs<sup>27</sup> database was created with proteins from UniRef90<sup>37</sup> and METdb<sup>38,39</sup>. This database was aligned to the MAGs' open-reading frames (min. 50 amino-acids) with mmseqs search (--start-sens 1 --sens-steps 3 -s 5 --translate 1 -e 0.00001 --min-length 50). For each MAG, the list of proteins with a match, hereafter referred as “protein hints”, were extracted with mmseqs convertalis followed by mmseqs createsubdb.
Two *ab initio* gene predictor were used. Augustus<sup>40</sup> v3.5.0 was run with the protein hints aligned on each MAG using miniprot<sup>41</sup> (introns up to 10 kb), while the best value for the --species parameter was estimated from a preliminary taxonomic annotation of the DNA-dependent RNA polymerase A, detected according to Delmont *et al.*<sup>3</sup> methodology. Then, a custom gene model was generated per MAG with the scripts provided by Augustus, which was used to generate a final gene prediction. The genes were also predicted in a single pass by GeneMark-EP+<sup>42</sup> v4.71, with protein hints aligned using ProtHint v2.6.0. Finally, both genes prediction were benchmarked with BUSCO<sup>43</sup> v5.3.2 in “transcriptome” mode against the database eukaryota_odb10. The GeneMark-EP+ gene prediction strategy resulted in MAGs of higher average completeness (Augustus: 34,5 %, GeneMark-EP+: 39,9 %) with marginal increase in duplication rate (Figure 3). Therefore, we retained the GeneMark-EP+ annotation as a reference.

## Single-cell sorting

Focussing on the small size fraction, subsamples of 50 mL were treated with 0.01% pluronic acid and diluted in TE buffer (10 mM Tris-HCl, 1 mM EDTA, pH 8, prefiltered to 0.2 µm), to achieve a flow rate of 1,500 events per second. These subsamples were stained with SYBR Green I (Thermo Fisher Scientific) at 1X final concentration for 10 min at room temperature in the dark. Flow cytometric analysis and cell sorting were conducted using a BD FACSAria Fusion SORP cytometer installed in a Baker BioProtect IV BioSafety cabinet. The instrument was equipped with a 100 µm nozzle and a 1.5 neutral density filter, and operated using BD FACSDiva 9 software
with sterile PBS as sheath fluid. Eukaryotic single cells were identified based on forward scatter (FSC), side scatter (SSC), and SYBR Green I fluorescence intensity, using a 488 nm (50 mW) laser with 502 nm longpass and 530/30 nm bandpass filters. Autotrophic and heterotrophic cells were distinguished by their chlorophyll content, detected using a 405 nm (50 mW) laser with 630 nm longpass and 660/20 nm bandpass filters.

Samples were analysed and sorted under sterile conditions at the lowest flow rate (maximum 1,500 events per second) in three successive stages: first in “continuous yield” mode, then in “continuous purity” mode, and finally in “single cell” mode. To promote cell settling at the bottom of the collection tubes, we sorted one cell together with 1000 Fluoresbrite YG Microspheres (0.50 µm; Polysciences, Warrington, PA, USA) into conical 96-well semi-skirted plates (Bio-Rad, Hercules, CA, USA) containing 5 µL of PBS. This resulted in an average amplification success rate of 30%.

## Single-cell genome amplification, screening, sequencing and assembly

Single-cell whole genome amplification was performed using the Single Cell GenomiPhi DNA Amplification kit (Cytiva, Marlborough, MA, USA) with minor modifications. Specifically, five freeze-thaw cycles (each lasting 4 min) were introduced prior to lysis with the kit's lysis buffer. Additionally, the amplification reaction was extended to 16 h at 30°C, instead of the standard 2 h. The amplified genomic DNA was diluted 100-fold and used as a template for downstream PCR amplification. Nearly full-length 18S rDNA genes were amplified (95°C for 5min followed by 30 cycles 1 min-95°C, 1 min-55°C, 1 min-72°C and a final extension 5 min-72°C) using the eukaryotic primers Euk328f (5'-ACCTGGTTGATCCTGCCAG-3')⁴⁴ and Euk329r (5'-TGATCCTTCYGCAAGGTTCAC-3')⁴⁴. PCR products were then cloned into pCRII-TOPO vectors (Invitrogen, Waltham, MA, USA) and transformed into *Escherichia coli* XL1-Blue competent cells⁴⁵. Clones were randomly picked from different plates and checked for 18S rDNA inserts by PCR amplification using flanking vector primers (M13f and M13r). After purification (PCR purification kit; Qiagen), PCR product were sequenced (Sanger) by Eurofins and were compared to the NT database with the NCBI BLAST web application.

A total of 19 SAGs were selected for further sequencing, with a particular focus on potential parasites such as Perkinsea, Aphelida, and Chytrids. Fragmentation was started with 250 ng of DNA to reach a 100-1000 bp size range using the E220 Covaris (Covaris, Woburn, MA, USA) instrument. Illumina libraries were then prepared using the NEBNEXT End Repair, dA-Tailing and Quick Ligation modules (New England Biolabs) and the NEXTFLEX Unique Dual Index Barcodes (Perkin Elmer). Ligated products were purified by two consecutive 1X AMPure XP clean up (Beckman Coulter), then amplified by 10 PCR cycles using the KAPA HiFi HotStart Library Amplification Kit (Roche), amplification was then purified by 0.8X AMPure XP purification. Libraries were done in duplicates for four SAGs, resulting in 23 items. The libraries were sequenced with an Illumina NovaSeq 6000 instrument (2x150 bp).

The single-cell genomes were assembled with SPAdes⁴⁶ v3.15.0 in single-cell mode (--sc). Contigs > 2.5 kb were kept and compared to the NT database by BLASTN⁴⁷ (-evalue 1.10⁻⁵ - numAlignment 5 -outfmt "6 std staxid") to identify and exclude contaminants (Human DNA, Bacteria, Archaea and Virus). Assembly quality was evaluated using QUAST⁴⁸ and BUSCO⁴³ v5.3.2
(-m genome and the eukaryota_odb10 dataset) and resulted in 14 assemblies corresponding to 11 SAGs that have been included in the dataset. Genes were predicted following the same strategy used for the MAGs with the difference that only GeneMark-EP+ was employed. SAG completion and redundancy were assessed with BUSCO⁴³ v5.3.2 (protein mode) and eukaryota_odb10 dataset.³⁰

## MAGs and SAGs tree based taxonomy and global distribution

The MAGs and SAGs taxonomy was further refined according to a phylogenomic approach. Following the methodology from Duncan et al.⁴, MAGs and SAGs proteins were first searched against the Phylosift⁴⁹ v1.0.1 marker set (phylosift search --isolate --besthit), together with a selection of eukaryote reference genomes downloaded from the NCBI (datasets⁵⁰ command-line tool, list of genomes in Data Record section for details). Fifty Phylosift markers which were present in at least 50% of the genomes were aligned (phylosift align --isolate --besthit). Trimmed alignments (trimAl⁵¹ v1.5 with --automated1) were concatenated for phylogenetic inference with IQ-TREE⁵² v2.2.3 (-m LG+F+I+R10 -B 2000 --alrt 2000 --bnni).

MAGs and SAGs proteomes were also searched for the large and small DNA-dependent RNA polymerases subunits. Considering the DNA-dependent RNA polymerase genes reference dataset from Delmont et al.³, proteins were split into 6 files, one for each subunit⁵³. The alignments of each subunit with MAFFT⁵⁴ v7.525 were cleaned with TrimAL⁵¹ v1.5 (--gt 0.5), and condensed in a HMM profile with hmmbuild²² v3.4. The 6 HMM profiles were searched against the MAGs, SAGs, the reference genomes used for the “Phylosift” analyses, plus the genomic reference database for marine species (METdb)³⁹, and freshwater MAGs from Lake Mendota⁵ (hmmsearch v3.4; e-value threshold of 1e⁻⁵⁰, retaining best-hit). Two sequences of the Malawimonada taxa were found directly in GenBank. Each marker was aligned separately with MAFFT⁵⁴ v7.525, sites with more than 50% gaps were trimmed with goalign⁵⁵ v0.3.7. Alignments were then concatenated and the tree was estimated by IQ-TREE⁵² v2.4.0 with the automatic model selection (-m MPF)⁵⁶. Support was computed with 1,000 replicates for both the ultrafast bootstrap⁵⁷ and the SH-like approximation likelihood ratio⁵⁸. The taxonomy of MAGs and SAGs was inferred from their placement in the trees. Taxonomies of MAGs and SAGs recovered in both approaches were congruent, however the DNA-dependent RNA polymerase-based phylogeny (Figure 3) was overall more resolution.

The presence of our MAGs and SAGs was also investigated in about 3,000 public metagenomes that could contain eukaryotic data. The metagenome reads were filtered with fastp⁵⁹ v0.23.4 authorising no N and an average quality greater than 30, prior their mapping with bwa mem²⁸ v0.7.17, followed by samtools view³³ (v1.16.1, parameters -F4 -q30 -bS) and samtools sort. Aligned reads were then filtered by msamtools filter³⁴ v1.1.0 (-b -l 50 -p 95 -z 80). The coverage was calculated per contig with samtools coverage and merged per assembly afterwards. MAGs or SAGs with fewer than 1,500 reads mapped across all metagenomes were removed from the analysis, as well as metagenomes with less than 100 reads mapped. A MAG or SAG was considered detected in a metagenome when its breadth of coverage was at least 10%. To give some more insights on MAGs and SAGs shared distribution across different ecosystems, we performed metagenome clustering
using read counts. The centered log ratio (CLR) transformation was applied to the mapped read
count per metagenome on MAGs and SAGs (package clr^60 v2.0-8 in R^61 v4.4.1), with a preliminary
replacement of zeros by the count zero multiplicative method (package zCompositions^62 v1.5.0-
4). Aitchison distance (i.e. Euclidian distance of CLR data) was computed and used for the meta-
genomes clustering with Ward.D2 method. The elbow method (fviz_nbclust()^63 v1.0.7) was used
to define optimal number of clusters, which was 3. These clusters have no apparent bias with either
metagenome size or size fraction filtering (Figure S1).

Data Records

All sequencing data were deposited in the ENA/SRA databases under the umbrella BioProject PRJEB61527^65, which gather reads from the 18S rDNA V4 and V9 data (PRJEB61511^66 and PRJEB94150^71 respectively), the metatranscriptomes (PRJEB61515^68), the metagenomes (PRJEB94020^69) and the SAGs reads (PRJEB94149^70).

The other types of data generated for this article are stored in the Zenodo record 17812131^64. It
contains i) the MAGs sequences and genes predictions (both with GeneMark and Augustus) in
GFF3 format, the transcript and protein sequences as fasta files; ii) the SAGs sequences; iii) a table
with MAGs and SAGs metadata (assembly statistics, BUSCO results, affiliation); iv) the DNA-
dependent RNA polymerase proteins HMM profiles from the 6 subunits used to build the dataset
for the phylogeny, as well as their individuals and concatenated alignments; v) the phylogeny built
with Phylosift; vi) a table with the MAGs and SAGs taxonomic annotations, vii) the 4
metagenomic co-assemblies (contigs ≥ 2.5 kb) with the mapping of their reads and the MetaBAT2
bins, available as Anvi'o "contigs" and "profiles databases", viii) the unigenes annotations
(taxonomic and functional) and the table of raw read count per sample, ix) information about the
public metagenomes that mapped on the MAGs and SAGs, in term of read count and genome
detection (breadth of coverage), plus the metadata associated.

Technical Validation

After extraction, the nucleic acids were quantified with the QuBit RNA HS and dsDNA HS assay
kits (Thermo Fisher Scientific) on a Qubit model 3 fluorometer, and the RNA integrity was
controlled with a TapeStation 2100 using the High Sensitivity RNA ScreenTape (Agilent, Santa
Clara, CA, USA) and the absence of DNA was checked by PCR. RNA integrity numbers ranged
from 7.3 to 9.4. For the DNA and RNA duplicates, 7 and 16 samples respectively failed to supply
enough material and thus were pooled for sequencing.

The cell sorting required precautions. All equipment was sterilised by autoclaving or UV decontamination, the cytometer was cleaned and prepared for aseptic sorting according to the manufacturer's instructions. The cytometer performance was verified and adjusted before each use with BD CS&T beads and drop delay was calibrated using BD FACS Accudrops Beads in autodelay mode. All signals were visualised on logarithmic scales, with the threshold set to the minimum value for the SYBR Green I parameter.
The sequencing platform (Genoscope) applied the following filters on the raw data, i) removing adapters and primers on the whole read and low quality nucleotides from both ends (Q-score lower than 20), ii) removing sequences between the second unknown nucleotide (N) and the end of the read, iii) discarding reads shorter than 30 nucleotides after trimming and iv) removing read pairs that come from the low-concentration spike-in library of Illumina PhiX Control. The small contigs (≤ 2.5 kb) were removed from the assemblies as advised by MetaBAT2<sup>35</sup>. Reads were mapped back on the assembly and alignments were filtered to retain those greater than 95% identity on at least 75 nt and the highest scoring hit was kept when multi-mapping occurred. MAGs and SAGs completion and redundancy scores were estimated by BUSCO and the dataset "eukaryota_odb10". The metatranscriptomic libraries were subject to rDNA (5S, 5.8S, 16S, 18S, 23S and 28S) contamination evaluation. The unigenes are the representative sequences from transcripts clustered by CD-Hit-EST >95% identity and >90% of the length of the smallest sequence. For the 18S rDNA V4 metabarcoding, a negative control was sequenced and the raw reads were deposited under the BioProject PRJEB61512.

## Data Availability

All sequencing data were deposited in the ENA/SRA databases under the BioProjects PRJEB61527<sup>65</sup> (Umbrella), PRJEB61511<sup>66</sup>, PRJEB61512<sup>67</sup>, PRJEB61515<sup>68</sup>, PRJEB94020<sup>69</sup>, PRJEB94149<sup>70</sup> and PRJEB94150<sup>71</sup>. The processed data are available in the Zenodo record 17812131<sup>64</sup>.

## Code Availability

Scripts related to the metabarcoding and metatranscriptomics data processing are available in the GitHub repository https://github.com/meb-team/MICROSTORE_MetaT-MetaB_Monjot_2023. Examples of codes and scripts used to process the metagenomic co-assemblies, the gene prediction, and for the mapping of public metagenomes, are available in the repository https://github.com/meb-team/datapaper_eukaryotes_pavin_2025

## References

1. Wu, D., Seshadri, R., Kyrpides, N. C. & Ivanova, N. N. A metagenomic perspective on the microbial prokaryotic genome census. *Sci. Adv.* **11**, https://doi.org/10.1126/sciadv.adq2166 (2025).

2. Wilken, S., Choi, C. J. & Worden, A. Z. Contrasting Mixotrophic Lifestyles Reveal Different Ecological Niches in Two Closely Related Marine Protists. *J. Phycol.* **56**, 52–67, https://doi.org/10.1111/jpy.12920 (2020).

3. Delmont, T. O. *et al*. Functional repertoire convergence of distantly related eukaryotic plankton lineages abundant in the sunlit ocean. *Cell Genomics* **2**, 100123, https://doi.org/10.1016/j.xgen.2022.100123 (2022).
4. Duncan, A. et al. Metagenome-assembled genomes of phytoplankton microbiomes from the Arctic and Atlantic Oceans. *Microbiome* **10**, 67, https://doi.org/10.1186/s40168-022-01254-7 (2022).

5. Krinos, A. I. et al. Time-series metagenomics reveals changing protistan ecology of a temperate dimictic lake. *Microbiome* **12**, 133, https://doi.org/10.1186/s40168-024-01831-y (2024).

6. Boukheloua, R. et al. Global freshwater distribution of Telonemia protists. *ISME J.* **18**, wrae177, https://doi.org/10.1093/ismejo/wrae177 (2024).

7. Lake Pavin. (Springer International Publishing, Cham, 2016). https://doi.org/10.1007/978-3-319-39961-4.

8. Eren, A. M. et al. Anvi’o: an advanced analysis and visualization platform for ‘omics data. *PeerJ* **3**, e1319, https://doi.org/10.7717/peerj.1319 (2015).

9. Carradec, Q. et al. A global ocean atlas of eukaryotic genes. *Nat. Commun.* **9**, 373, https://doi.org/10.1038/s41467-017-02342-1 (2018).

10. Monjot, A. et al. Functional diversity of microbial eukaryotes in a meromictic lake: Coupling between metatranscriptomic and a trait-based approach. *Environ. Microbiol.* **25**, 3406–3422, https://doi.org/10.1111/1462-2920.16531 (2023).

11. Monjot, A., Rousseau, J., Bittner, L. & Lepère, C. Metatranscriptomes-based sequence similarity networks uncover genetic signatures within parasitic freshwater microbial eukaryotes. *Microbiome* **13**, 43, https://doi.org/10.1186/s40168-024-02027-0 (2025).

12. Schoenle, A. et al. Protist genomics: key to understanding eukaryotic evolution. *Trends Genet.* https://doi.org/10.1016/j.tig.2025.05.004 (2025).

13. Parada, A. E., Needham, D. M. & Fuhrman, J. A. Every base matters: assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. *Environ. Microbiol.* **18**, 1403–1414, https://doi.org/10.1111/1462-2920.13023 (2016).

14. Mangot, J.-F. et al. Short-term dynamics of diversity patterns: evidence of continual reassembly within lacustrine small eukaryotes. *Environ. Microbiol.* **15**, 1745–1758, https://doi.org/10.1111/1462-2920.12065 (2013).

15. Chauvet, M., Debroas, D., Moné, A., Dubuffet, A. & Lepère, C. Temporal variations of Microsporidia diversity and discovery of new host–parasite interactions in a lake ecosystem. *Environ. Microbiol.* **24**, 1672–1686, https://doi.org/10.1111/1462-2920.15950 (2022).

16. Amaral-Zettler, L. A., McCliment, E. A., Ducklow, H. W. & Huse, S. M. A Method for Studying Protistan Diversity Using Massively Parallel Sequencing of V9 Hypervariable Regions of Small-Subunit Ribosomal RNA Genes. *PLOS ONE* **4**, e6372, https://doi.org/10.1371/journal.pone.0006372 (2009).

17. López-García, P., Philippe, H., Gail, F. & Moreira, D. Autochthonous eukaryotic diversity in hydrothermal sediment and experimental microcolonizers at the Mid-Atlantic Ridge. *Proc. Natl. Acad. Sci.* **100**, 697–702, https://doi.org/10.1073/pnas.0235779100 (2003).
18. Zerbino, D. R. & Birney, E. Velvet: algorithms for de novo short read assembly using de Bruijn graphs. *Genome Res.* **18**, 821–829, https://doi.org/10.1101/gr.074492.107 (2008).

19. Schulz, M. H., Zerbino, D. R., Vingron, M. & Birney, E. Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. *Bioinformatics* **28**, 1086–1092, https://doi.org/10.1093/bioinformatics/bts094 (2012).

20. Haas, B. J. TransDecoder, https://github.com/TransDecoder/TransDecoder (2018).

21. Eberhardt, R. Y. et al. AntiFam: a tool to help identify spurious ORFs in protein annotation. *Database* **2012**, bas003, https://doi.org/10.1093/database/bas003 (2012).

22. Eddy, S. R. Accelerated Profile HMM Searches. *PLOS Comput. Biol.* **7**, e1002195, https://doi.org/10.1371/journal.pcbi.1002195 (2011).

23. Aramaki, T. et al. KofamKOALA: KEGG Ortholog assignment based on profile HMM and adaptive score threshold. *Bioinformatics* **36**, 2251–2252, https://doi.org/10.1093/bioinformatics/btz859 (2020).

24. Mistry, J. et al. Pfam: The protein families database in 2021. *Nucleic Acids Res.* **49**, D412–D419, https://doi.org/10.1093/nar/gkaa913 (2021).

25. Gene Ontology Consortium. Pfam2GO. https://doi.org/10.5281/zenodo.6399963 (2022).

26. Levy Karin, E., Mirdita, M. & Söding, J. MetaEuk—sensitive, high-throughput gene discovery, and annotation for large-scale eukaryotic metagenomics. *Microbiome* **8**, 48, https://doi.org/10.1186/s40168-020-00808-x (2020).

27. Steinegger, M. & Söding, J. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. *Nat. Biotechnol.* **35**, 1026–1028, https://doi.org/10.1038/nbt.3988 (2017).

28. Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. Preprint at https://doi.org/10.48550/arXiv.1303.3997 (2013).

29. Benoit, G. et al. SimkaMin: fast and resource frugal de novo comparative metagenomics. *Bioinformatics* **36**, 1275–1276, https://doi.org/10.1093/bioinformatics/btz685 (2020).

30. Li, D., Liu, C.-M., Luo, R., Sadakane, K. & Lam, T.-W. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. *Bioinformatics* **31**, 1674–1676, https://doi.org/10.1093/bioinformatics/btv033 (2015).

31. Hyatt, D. et al. Prodigal: prokaryotic gene recognition and translation initiation site identification. *BMC Bioinformatics* **11**, 119, https://doi.org/10.1186/1471-2105-11-119 (2010).

32. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics* **25**, 1754–1760, https://doi.org/10.1093/bioinformatics/btp324 (2009).

33. Danecek, P. et al. Twelve years of SAMtools and BCFtools. *GigaScience* **10**, giab008, https://doi.org/10.1093/gigascience/giab008 (2021).

34. Arumugam, M. msamtools, https://github.com/arumugamlab/msamtools (2025).

35. Kang, D. D. et al. MetaBAT 2: an adaptive binning algorithm for robust and efficient genome
reconstruction from metagenome assemblies. *PeerJ* **7**, e7359, https://doi.org/10.7717/peerj.7359 (2019).

36. Delmont, T. O. Assessing the completion of eukaryotic bins with anvi'o. *Meren Lab* https://merenlab.org/2018/05/05/eukaryotic-single-copy-core-genes/ (2018).

37. Suzek, B. E., Huang, H., McGarvey, P., Mazumder, R. & Wu, C. H. UniRef: comprehensive and non-redundant UniProt reference clusters. *Bioinformatics* **23**, 1282–1288, https://doi.org/10.1093/bioinformatics/btm098 (2007).

38. Keeling, P. J. et al. The Marine Microbial Eukaryote Transcriptome Sequencing Project (MMETSP): Illuminating the Functional Diversity of Eukaryotic Life in the Oceans through Transcriptome Sequencing. *PLOS Biol.* **12**, e1001889, https://doi.org/10.1371/journal.pbio.1001889 (2014).

39. Niang, G. et al. METdb: A Genomic Reference Database for Marine Species. *F1000Research* **9**, https://doi.org/10.7490/f1000research.1118000.1 (2020).

40. Stanke, M. et al. AUGUSTUS: ab initio prediction of alternative transcripts. *Nucleic Acids Res.* **34**, W435–W439, https://doi.org/10.1093/nar/gkl200 (2006).

41. Li, H. Protein-to-genome alignment with miniprot. *Bioinformatics* **39**, btad014, https://doi.org/10.1093/bioinformatics/btad014 (2023).

42. Brůna, T., Lomsadze, A. & Borodovsky, M. GeneMark-EP+: eukaryotic gene prediction with self-training in the space of genes and proteins. *NAR Genomics Bioinforma.* **2**, lqaa026, https://doi.org/10.1093/nargab/lqaa026 (2020).

43. Manni, M., Berkeley, M. R., Seppey, M., Simão, F. A. & Zdobnov, E. M. BUSCO Update: Novel and Streamlined Workflows along with Broader and Deeper Phylogenetic Coverage for Scoring of Eukaryotic, Prokaryotic, and Viral Genomes. *Mol. Biol. Evol.* **38**, 4647–4654, https://doi.org/10.1093/molbev/msab199 (2021).

44. Moon-van der Staay, S. Y. et al. Abundance and Diversity of Prymnesiophytes in the Pico-plankton Community from the Equatorial Pacific Ocean Inferred from 18S rDNA Sequences. *Limnol. Oceanogr.* **45**, 98–109 (2000).

45. Inoue, H., Nojima, H. & Okayama, H. High efficiency transformation of *Escherichia coli* with plasmids. *Gene* **96**, 23–28, https://doi.org/10.1016/0378-1119(90)90336-P (1990).

46. Prjibelski, A., Antipov, D., Meleshko, D., Lapidus, A. & Korobeynikov, A. Using SPAdes De Novo Assembler. *Curr. Protoc. Bioinforma.* **70**, e102, https://doi.org/10.1002/cpbi.102 (2020).

47. Camacho, C. et al. BLAST+: architecture and applications. *BMC Bioinformatics* **10**, 421, https://doi.org/10.1186/1471-2105-10-421 (2009).

49. Darling, A. E. et al. PhyloSift: phylogenetic analysis of genomes and metagenomes. *PeerJ* **2**, e243, https://doi.org/10.7717/peerj.243 (2014).

50. O'Leary, N. A. et al. Exploring and retrieving sequence and metadata for species across the tree of life with NCBI Datasets. *Sci. Data* **11**, 732, https://doi.org/10.1038/s41597-024-03571-0
γ (2024).

51. Capella-Gutiérrez, S., Silla-Martínez, J. M. & Gabaldón, T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. *Bioinformatics* **25**, 1972–1973, https://doi.org/10.1093/bioinformatics/btp348 (2009).

52. Minh, B. Q. et al. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era. *Mol. Biol. Evol.* **37**, 1530–1534, https://doi.org/10.1093/molbev/msaa015 (2020).

53. Werner, F. & Grohmann, D. Evolution of multisubunit RNA polymerases in the three domains of life. *Nat. Rev. Microbiol.* **9**, 85–98, https://doi.org/10.1038/nrmicro2507 (2011).

54. Katoh, K. & Standley, D. M. MAFFT Multiple Sequence Alignment Software Version 7: Improvements in Performance and Usability. *Mol. Biol. Evol.* **30**, 772–780, https://doi.org/10.1093/molbev/mst010 (2013).

55. Lemoine, F. & Gascuel, O. Gotree/Goalign: toolkit and Go API to facilitate the development of phylogenetic workflows. *NAR Genomics Bioinforma.* **3**, lqab075, https://doi.org/10.1093/nargab/lqab075 (2021).

56. Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A. & Jermiin, L. S. ModelFinder: fast model selection for accurate phylogenetic estimates. *Nat. Methods* **14**, 587–589, https://doi.org/10.1038/nmeth.4285 (2017).

57. Minh, B. Q., Nguyen, M. A. T. & von Haeseler, A. Ultrafast Approximation for Phylogenetic Bootstrap. *Mol. Biol. Evol.* **30**, 1188–1195, https://doi.org/10.1093/molbev/mst024 (2013).

58. Guindon, S. et al. New Algorithms and Methods to Estimate Maximum-Likelihood Phylogenies: Assessing the Performance of PhyML 3.0. *Syst. Biol.* **59**, 307–321, https://doi.org/10.1093/sysbio/syq010 (2010).

59. Chen, S. Ultrafast one-pass FASTQ data preprocessing, quality control, and deduplication using fastp. *iMeta* **2**, e107, https://doi.org/10.1002/imt2.107 (2023).

60. Boogaart, K. G. van den, Tolosana-Delgado, R. & Bren, M. Compositions: Compositional Data Analysis, https://doi.org/10.32614/CRAN/package.compositions (2024).

61. R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, Vienna, Austria, 2024), https://www.R-project.org/.

62. Palarea-Albaladejo, J. & Martín-Fernández, J. A. zCompositions – R package for multivariate imputation of left-censored data under a compositional approach. *Chemom. Intell. Lab. Syst.* **143**, 85–96, https://doi.org/10.1016/j.chemolab.2015.02.019 (2015).

63. Kassambara, A. & Mundt, F. Factoextra: Extract and Visualize the Results of Multivariate Data Analyses, https://doi.org/10.32614/CRAN/package.factoextra (2020).

64. Courtine, D. et al. MICROSTORE: 'omic approaches to decipher MICrobial eukaRyOteS funcTiOn in fReshwater lake Ecosystems. Zenodo https://doi.org/10.5281/zenodo.17812131 (2025).

65. ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB61527 (2023).
66. ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB61511 (2023).

67. ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB61512 (2023).

68. ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB61515 (2023).

69. ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB94020 (2025).

70. ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB94149 (2025).

71. ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB94150 (2025).

## Author Contributions

CL, DD, GB – Conceptualisation and funding acquisition; DC, IW, JMA, CDS – Data Curation; CL, AMe, IW, HB, JC, CC – Investigation; CL, AMe, HB, JC, DC, GB, DD – Methodology; DC, JMA, CDS – Software; DC, GB, AMt – Visualisation; DC, GB, CL, DD, CC – Writing original draft; All authors – Writing, review and edition

## Competing Interests

The authors declare no competing interests.

## Acknowledgements

We thank Agnès Vellet, now retired, for her dedicated work on single-cell genomic development. We are also grateful to the AuBi platform and the Mésocentre Clermont-Auvergne of the Université Clermont Auvergne for providing computing and storage resources facilities. This work was supported by the Genoscope, the Commissariat à l'Énergie Atomique et aux Énergies Alternatives (CEA) and France Génomique (ANR-10-INBS-09-08). This project benefited from equipment of CYSTEM platform UCA-PARTNER (Université Clermont Auvergne, Clermont-Ferrand, France).

## Figure legends

Figure 1: Graphical overview of the sampling, data production and analysis. Top, Lake Pavin sampling conditions (2018) and size class filtering strategy. Bottom, bioinformatics analysis for each approach: amplicon based-, metagenomic, metatranscriptomic, and single-cell amplified genome sequencing. Synthetic overview of the data is provided by placing the MAGs and SAGs in an evolutionary framework, and representing their distribution in different aquatic environments via a heatmap. The functional potential of the microbial communities can be derived from the annotation of the unigenes catalogue as found in Monjot et al.¹⁰

Figure 2: Principal Coordinates Analysis of the metagenomes. Relations between the metagenomic samples based on a Bray-Curtis distance matrix of their sequence composition as measured with SimkaMin²⁹ on 50 million reads and 10 million k-mers per
sample.

Figure 3: BUSCO completion and redundancy scores of the MAGs and SAGs based on the eukaryota_odb10 single-copy core genes set.

Completion is the proportion of *unique + duplicated* genes found; redundancy is the proportion of *duplicated* genes.

Figure 4: Phylogenetic tree of the concatenated DNA-dependent RNA polymerase protein sequences from eukaryotes.

The maximum likelihood phylogenetic tree of the concatenated sub-units from the DNA-dependent RNA polymerase, 6 genes in total, including MAGs and SAGs of this work, Tara Ocean MAGs³, METdb transcriptomes³⁹, genomes from the public databases and MAGs from the Lake Mendota⁵, was generated using 7,659 sites in the alignment and the LG + F + I + R10 model. The clade of Opisthokonta was used as the outgroup. Phylogenetic supports were considered high (aLRT ≥ 80 and UFBoot ≥ 95), medium (aLRT ≥ 80 or UFBoot ≥ 95), or low (aLRT < 80 and UFBoot < 95).

Figure 5: Distribution of MAGs and SAGs in water metagenomes

a, Detection of MAGs and SAGs per metagenomic dataset. These datasets were clustered based on the number of read (CLR transformation) mapping on the MAGs and SAGs. Each square summarizes the number of metagenomes that belong to an ecosystem type (bottom) and the number of metagenomes in which at least one MAG or SAG was detected (top). A genome is considered detected in a metagenome when its breadth of coverage is ≥10%. Empty square means no data. b, UpSet plot that describes the clusters composition in term of number of MAG or SAG detected, i.e. at least one genome with a breadth of coverage ≥10% in ≥1 metagenome. c, UpSet plot of the phyla present in the different clusters. for b and c, empty sets are not shown.

## Figure disposition

Figure 1 goes after the sentence "Water sampling was carried out in both layers (at 9 m and 80 m) by day and night, in April, June, September, and November, and for two size classes (0.65-10 and 10-50 µm)"

Figure 2 goes at the end of the section "Meta-omic reads assembly and binning"

Figure 3 goes at the end of the section "MAGs gene prediction"

Figure 4 goes at the end of the section "Single-cell genome amplification, screening, sequencing and assembly"

Figure 5 goes at the end of the section "MAGs and SAGs tree based taxonomy and global distribution"

Table 1 goes at the end of the section "Background and summary"
<caption>Table 1: Summary of data generated per type, in million pairs of reads</caption>

The numbers reflect the deposited data, after quality control from the sequencing platform. The asterisk (*) indicates samples with merged duplicates.

<table><thead><tr><th>Month</th><th>Zon e</th><th>Size fraction</th><th>Period of the day</th><th>Meta-genomic</th><th>Meta-transcriptomic</th><th>Metabarcoding 18S V4</th><th>Metabarcoding 18S V9</th></tr></thead><tbody><tr><td rowspan="6">April</td><td rowspan="2">Oxi c</td><td>Small</td><td>Day</td><td>305</td><td>289</td><td>0.549</td><td>0.678</td></tr><tr><td>Night</td><td>127*</td><td>262</td><td>0.227*</td><td>0.219*</td></tr><tr><td rowspan="2">Large</td><td>Day</td><td>332</td><td>341</td><td>0.281*</td><td>0.386</td></tr><tr><td>Night</td><td>339</td><td>332</td><td>0.494</td><td>0.440</td></tr><tr><td rowspan="2">An-oxic</td><td>Small</td><td>Day</td><td>341</td><td>174*</td><td>0.525</td><td>0.238</td></tr><tr><td>Night</td><td>126*</td><td>180*</td><td>0.215*</td><td>0.128*</td></tr><tr><td rowspan="6">June</td><td rowspan="2">Oxi c</td><td>Small</td><td>Day</td><td>287</td><td>239</td><td>0.455</td><td>0.273</td></tr><tr><td>Night</td><td>371</td><td>297</td><td>0.486</td><td>0.312</td></tr><tr><td rowspan="2">Large</td><td>Day</td><td>378</td><td>233*</td><td>0.408</td><td>0.419</td></tr><tr><td>Night</td><td>341</td><td>226*</td><td>0.478</td><td>0.267</td></tr><tr><td rowspan="2">An-oxic</td><td>Small</td><td>Day</td><td>350</td><td>304</td><td>0.441</td><td>0.333</td></tr><tr><td>Night</td><td>385</td><td>193*</td><td>0.397</td><td>0.443</td></tr><tr><td rowspan="6">Sep-tem-ber</td><td rowspan="2">Oxi c</td><td>Small</td><td>Day</td><td>453</td><td>325</td><td>0.406</td><td>0.306</td></tr><tr><td>Night</td><td>241</td><td>326</td><td>0.397</td><td>0.246</td></tr><tr><td rowspan="2">Large</td><td>Day</td><td>304</td><td>156*</td><td>0.496</td><td>0.248</td></tr><tr><td>Night</td><td>347</td><td>132*</td><td>0.452</td><td>0.234</td></tr><tr><td rowspan="2">An-oxic</td><td>Small</td><td>Day</td><td>235</td><td>242</td><td>0.440</td><td>0.342</td></tr><tr><td>Night</td><td>268</td><td>98*</td><td>0.405</td><td>0.543</td></tr><tr><td rowspan="6">No-vem-ber</td><td rowspan="2">Oxi c</td><td>Small</td><td>Day</td><td>164*</td><td>144*</td><td>0.297*</td><td>0.220</td></tr><tr><td>Night</td><td>262</td><td>120*</td><td>0.460</td><td>0.220</td></tr><tr><td rowspan="2">Large</td><td>Day</td><td>289</td><td>389</td><td>0.373</td><td>0.250</td></tr><tr><td>Night</td><td>307</td><td>371</td><td>0.452</td><td>0.652</td></tr><tr><td rowspan="2">An-oxic</td><td>Small</td><td>Day</td><td>344</td><td>367</td><td>0.479</td><td>0.433</td></tr><tr><td>Night</td><td>189</td><td>370</td><td>0.473</td><td>0.358</td></tr><tr><td rowspan="6">Novem-ber</td><td rowspan="2">Oxi c</td><td>Small</td><td>Day</td><td>280</td><td>317</td><td>0.435</td><td>0.404</td></tr><tr><td>Night</td><td>251</td><td>345</td><td>0.540</td><td>0.398</td></tr><tr><td rowspan="2">Large</td><td>Day</td><td>106*</td><td>184*</td><td>0.278*</td><td>0.142*</td></tr><tr><td>Night</td><td>100*</td><td>178*</td><td>0.248*</td><td>0.196*</td></tr></tbody></table>
<table><tr><td>Total read pairs</td><td>8,853</td><td>7,918</td><td>13.16</td><td>10.34</td></tr><tr><td>Average read pairs per biological sample</td><td>277</td><td>247</td><td>0.41</td><td>0.33</td></tr><tr><td>Number of samples sequenced</td><td>57</td><td>48</td><td>57</td><td>56</td></tr></table>

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>

## ILLUMINA SEQUENCING & BIOINFORMATICS ANALYSIS 3

<figure><img src="image_5.png" /></figure>

## PHYLOGENY, DISTRIBUTION & FUNCTIONAL ANNOTATION 4

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>

**Size fraction**

• Small

● Large

**Month**

● April

● June

● September

● November

**Depth**

o Upper layer (oxic)

● Lower layer (anoxic)
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /></figure>

② Sources

METDB

This project

Mendota

Tara Ocean

Public genomes
a

<figure><img src="image_11.png" /></figure>

Number of metagenome with ≥ 1 MAG/SAG detected (≥10%)

55

408

Total number of metagenomes

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>