Mickael Le Gac
email: mickael.le.gac@ifremer.fr

Lou Mary

Gabriel Metegnier

Julien Quéré

Raffaele Siano

Francisco Rodríguez

Christophe Destombe

Marc Sourisseau

Mickael; Ifremer Le Gac

ODE/Dyneco/ Pelagos Mary

ODE/Dyneco/ Lou ; Ifremer

Pelagos Metegnier

ODE/ Gabriel ; Ifremer

Dyneco / Pelagos Quéré

Julien ; Ifremer

Rodríguez Hernández

Strong population genomic

Keywords: Adaptation, Population Genetics -Empirical, Protists, Metagenomics, Transcriptomics Dinoflagellate, Population Genomics, Harmful Algal Blooms, meta-transcriptomic

come

Introduction Microbial census population sizes are often huge and difficult to apprehend. This is especially the case for pelagic marine microbes. For instance, the smallest and most abundant photosynthetic organism on the planet, the cyanobacteria Prochlorococcus has a census population size estimated to be around 10 27 cells [START_REF] Biller | Prochlorococcus : The structure and function of collective diversity[END_REF]. Even without going to such extreme values, numerous marine protists (a polyphyletic group of organisms encompassing all unicellular eukaryotes) may easily reach cell densities of thousands of cells per liter, representing several hundreds of billions of cells in a small bay. With such population sizes, adequate sampling to investigate population genetic diversity and structure from individual cells is a non-trivial problem. A problem amplified by the need to isolate cells from natural samples and initiate clonal cultures to obtain sufficient genetic material for genotyping. With such constraints, even for species relatively easy to cultivate, genotyping tens of clones per natural population quickly becomes a daunting task and working on difficult to cultivate species is impossible. Strategies based on single cell genotyping might be an alternative but there are technical and financial hurdles preventing their widespread and high throughput development. However, community wide samples, containing the genetic material of several thousands and even millions of protist cells may easily be obtained by filtering a few liters of water. Such environmental DNA (eDNA) or RNA (eRNA) samples are classically used to infer community composition following targeted sequencing (metabarcoding) or to investigate functional genomic aspects of natural communities using non-targeted sequencing (metagenomic or metatranscriptomic). Here, we explore the use of metatranscriptomic (metaT) datasets to infer protist population genetics for species of interest. The metaT datasets are composed of mRNA sequences extracted from natural communities and give access to a wide range of Single Nucleotide Polymorphism (SNPs) markers. A priori, the main advantages of such approaches are the relative simplicity of the sample processing (filtering water instead of isolation and F o r R e v i e w O n l y cultivation steps) combined with the ability to infer population allelic frequencies from samples containing thousands of cells. Another non-negligible advantage is the absence of targeted sequencing, enabling reuse of previously obtained datasets. This may be especially interesting as metaT datasets obtained worldwide are quickly accumulating [START_REF] Alexander | Metatranscriptome analyses indicate resource partitioning between diatoms in the field[END_REF][START_REF] Alexander | Functional group-specific traits drive phytoplankton dynamics in the oligotrophic ocean[END_REF][START_REF] Carradec | A global ocean atlas of eukaryotic genes[END_REF][START_REF] Cohen | Diatom Transcriptional and Physiological Responses to Changes in Iron Bioavailability across Ocean Provinces[END_REF][START_REF] Gong | Eukaryotic phytoplankton community spatiotemporal dynamics as identified through gene expression within a eutrophic estuary[END_REF][START_REF] Hu | Shifting metabolic priorities among key protistan taxa within and below the euphotic zone[END_REF][START_REF] Ji | Metatranscriptome analysis reveals environmental and diel regulation of a Heterosigma akashiwo (raphidophyceae) bloom[END_REF][START_REF] Lambert | The dynamic trophic architecture of openocean protist communities revealed through machine-guided metatranscriptomics[END_REF][START_REF] Lampe | Divergent gene expression among phytoplankton taxa in response to upwelling[END_REF][START_REF] Marchetti | Comparative metatranscriptomics identifies molecular bases for the physiological responses of phytoplankton to varying iron availability[END_REF][START_REF] Metegnier | Species specific gene expression dynamics during harmful algal blooms[END_REF][START_REF] Wurch | Transcriptional Shifts Highlight the Role of Nutrients in Harmful Brown Tide Dynamics[END_REF][START_REF] Zhang | Metatranscriptomic Signatures Associated With Phytoplankton Regime Shift From Diatom Dominance to a Dinoflagellate Bloom[END_REF]. The main drawback of such approach, but also of all potential approaches based on eDNA/RNA, is that it relies on the direct estimation of allelic frequencies in a population, without obtaining individual genotypes, limiting the type of possible analyses. Moreover, population genetic analyses from metaT datasets may only be considered if the mapping of the environmental reads is specific and the estimated allelic frequencies are precise and unbiased. Mapping environmental reads to a metareference composed of the reference transcriptomes of several hundreds of species has previously been shown to be a good way to ensure specific mapping [START_REF] Metegnier | Species specific gene expression dynamics during harmful algal blooms[END_REF]. For marine protists, this reference corresponds to the concatenation of the species specific reference transcriptomes obtained from the Marine Micro-eukaryote Transcriptome Project [START_REF] Keeling | The Marine Microbial Eukaryote Transcriptome Sequencing Project (MMETSP) : Illuminating the Functional Diversity of Eukaryotic Life in the Oceans through Transcriptome Sequencing[END_REF]. Quantifying the precision of allelic frequency estimates in relation to sequencing depth is especially critical. Indeed, in a metaT sample, the sequencing depth of a given species strongly depends on the relative abundance of this species in a community at the time of sampling, and the sequencing depth of a given gene strongly depends on the relative expression of this gene. Using a resampling approach, we explore how allele frequency estimates are impacted by sequencing depth. Using metaT datasets, allelic frequencies are based on mRNA, not DNA. This may also be a critical point if differential gene expression tends to bias the estimated allelic frequencies, for instance in case of allele specific gene expression. This aspect is explored by comparing allelic frequencies inferred after either genotyping individual cells or pooling RNAseq reads obtained from these same strains. The former enabled to determine the true allelic frequency of the population of

F o r R e v i e w O n l y 5 99
cultivated strains, while the latter included any potential bias due to differential gene expression.

The dinoflagellate Alexandrium minutum is responsible for harmful algal blooms that may transiently dominate local micro-eukaryote coastal communities [START_REF] Lewis | A review of the global distribution of Alexandrium minutum (Dinophyceae) and comments on ecology and associated paralytic shellfish toxin profiles, with a focus on Northern Europe[END_REF]. During these events cell densities may reach several million of cells per liter [START_REF] Garcés | Relationship between vegetative cells and cyst production during Alexandrium minutum bloom in Arenys de Mar harbour (NW Mediterranean)[END_REF]. This species produces Paralytic Shellfish Toxins (PSTs) that tend to bioaccumulate in the trophic network and more precisely in shellfish with potential sanitary and socio-economic impacts [START_REF] Ben-Gigirey | First Report of Paralytic Shellfish Toxins in Marine Invertebrates and Fish in Spain[END_REF][START_REF] Nogueira | HABs in coastal upwelling systems : Insights from an exceptional red tide of the toxigenic dinoflagellate Alexandrium minutum[END_REF]. Using 18S and ITS markers, two main clades were identified, one cosmopolite clade and the other restricted to the Southern Pacific [START_REF] Lilly | Phylogeny, biogeography, and species boundaries within the Alexandrium minutum group[END_REF][START_REF] Mccauley | BIOGEOGRAPHIC ANALYSIS OF THE GLOBALLY DISTRIBUTED HARMFUL ALGAL BLOOM SPECIES ALEXANDRIUM MINUTUM (DINOPHYCEAE) BASED ON rRNA GENE SEQUENCES AND MICROSA℡LITE MARKERS1[END_REF]. Within the cosmopolite clade, a global geographic structure was suggested by microsatellites [START_REF] Mccauley | BIOGEOGRAPHIC ANALYSIS OF THE GLOBALLY DISTRIBUTED HARMFUL ALGAL BLOOM SPECIES ALEXANDRIUM MINUTUM (DINOPHYCEAE) BASED ON rRNA GENE SEQUENCES AND MICROSA℡LITE MARKERS1[END_REF]. Still using microsatellites but at a more local scale, a strong genetic structure was identified in the Mediterranean Sea. This structuration pattern was at least partly compatible with hydrodynamics [START_REF] Casabianca | Population genetic structure and connectivity of the harmful dinoflagellate Alexandrium minutum in the Mediterranean Sea[END_REF]. In the South-West of the English Channel a moderate spatio-temporal structure was identified across a couple of years among two estuaries less than 200 km apart [START_REF] Dia | Spatiotemporal changes in the genetic diversity of harmful algal blooms caused by the toxic dinoflagellate Alexandrium minutum[END_REF]. Finally, two highly divergent populations, that may be considered as cryptic species, were identified in Noth-Western Europe using a SNP based approach (Le Gac et al., 2016c). Here taking as a model system A. minutum populations from Western Europe, we: 1. Pooled and resampled RNAseq reads obtained from strains to quantify how differential gene expression may bias observed allelic frequencies and to estimate the relationship between coverage and allelic frequency precision. 2. Compared A. minutum genetic structure inferred from strains and metaT datasets. 3. Determined coding genome wide genetic diversity and divergence among A. minutum populations using metaT datasets.

F o r R e v i e w O n l y
Material and Methods 1. Strains Each strain corresponded to a clonal culture initiated after micropipetting a single cell into fresh culture medium under an inverted microscope. These clonal strains are haploid. Origin of strains is indicated in Supplementary Table 1 and Figure 1. Out of the 37 strains used in the present study, 25 were isolated from water samples. The twelve others were isolated after germination of resting cysts from dated sediment cores, as indicated in [START_REF] Delebecq | Revival of Ancient Marine Dinoflagellates Using Molecular Biostimulation[END_REF]. Strains were grown to late exponential phase in K medium. Cultures were centrifuged at 4500 g for 8 min. RNA extraction occurred either directly after centrifugation, or cell pellets were frozen into liquid nitrogen with RNA Later and stored at -80°C until RNA extraction 2. Metatranscriptomic samples MetaT samples were obtained during A. minutum blooms in Western Europe from 2013 to 2018 by filtering water on 20 µm polycarbonate filters using a peristaltic pump. The filters were frozen into liquid nitrogen with RNA Later and stored at -80°C until RNA extraction. Information regarding the 77 metaT samples are indicated in Supplementary Table 2 andFigure 1. 3. RNA extraction and library preparation Samples were ultra-sonicated on ice in extraction buffers. RNA was extracted using either the Qiagen Rneasy Plus Mini kit or NucleoSpin® RNA Plus kit (Macherey-Nagel) Kit (Supplementary Table 1 and2). Library prepared with either the Illumina Truseq mRNA V2 kit or Illumina mRNA TruSeq stranded kit. Samples were sequenced at Get-PlaGe France Genomics sequencing platform (Toulouse, France) on Illumina HiSeq 2000/2500 2*100 pb or HiSeq 3000 2*150 pb. Raw sequencing reads are available in public databases [START_REF] Metegnier | Metatranscriptomic sequences from Alexandrium minutum blooms sampled in situ in the bay of Brest (France) between 2013 and 2015[END_REF], Le Gac et al. 2016a [START_REF] Bolger | Trimmomatic : A flexible trimmer for Illumina sequence data[END_REF] was used to trim ambiguous, low quality reads and sequencing adapters with parameters ILLUMINACLIP: Adapt.fasta:2:30:10:8:TRUE LEADING:3 TRAILING:3 MAXINFO:40:0.5 MINLEN:80 for 2*150 or MINLEN:60 for 2*100 reads.

Generating simulated population datasets from strains

Using bash scripts, 5e+06 reads were subsampled for each trimmed strain fastq files. They were concatenated to obtain a forward and a reverse fastq files simulating a population of strains. These files were subsampled in order to obtain 10 replicate files for seven coverage levels corresponding to a total of 1e+04, 5e+04, 1e+05, 5e+05, 1e+06, 5e+06 and 1e+07 reads. This was done separately for the 18 strains sequenced using 2*100bp reads and the 19 strains sequenced using 2*150bp reads (Supplementary Table 1).

Aligning reads

The strain as well as the simulated population samples were aligned to the A. minutum reference transcriptome previously developed and corresponding to 153,222 contigs (Le Gac et al., 2016b,c). The metaT datasets were aligned to a metareference corresponding to the concatenation of 313 species specific reference transcriptomes, representing 213 unique genus. It corresponded to the resources developed during the Marine Micro-eukaryote Transcriptome Project [START_REF] Keeling | The Marine Microbial Eukaryote Transcriptome Sequencing Project (MMETSP) : Illuminating the Functional Diversity of Eukaryotic Life in the Oceans through Transcriptome Sequencing[END_REF] and also included the A. minutum reference transcriptome [START_REF] Metegnier | Species specific gene expression dynamics during harmful algal blooms[END_REF]. Alignments were performed using BWA-MEM [START_REF] Li | Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM[END_REF]. Only reads with a mapping score >10 were retained. Pairs for which the two reads did not map concordantly on the same transcript were removed. Samtools [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF] was used to sort and index bam files.

. Identifying SNPs

A single nucleotide polymorphism (SNP) database was developed using the strain data following the approach proposed previously (Le Gac et al., 2016c). Briefly, FreeBayes [START_REF] Garrison | Haplotype-based variant detection from short-read sequencing[END_REF] was run twice using the 37 strains by enforcing haploidy and diploidy. As in culture conditions, A. minutum cells are in a vegetative, haploid stage, only SNPs detected using the haploidy enforced run and identified as haploid using the diploidy enforced run were considered. This was done to exclude potential intragenomic variability due to multicopy genes. The genotypes of each of the strains were obtained using vcftools [START_REF] Danecek | The variant call format and VCFtools[END_REF], for SNPs (excluding indels) displaying two alleles, a quality criterion >40, and covered more than 10 times in each of the 37 strains (no missing data). This database was composed of 227,829 SNPs. For these same SNPs, the allelic frequencies of the simulated population samples and of the natural populations based on the metaT samples were obtained using FreeBayes [START_REF] Garrison | Haplotype-based variant detection from short-read sequencing[END_REF] enforcing diploidy, and extracting coverage for each of the two alleles at the various SNP positions. In addition, for the metaT samples, A. minutum allelic frequencies were also estimated using a de novo approach, i.e. without restricting the analysis to the SNPs identified using the strains. It was performed for A. minutum contigs using Freebayes [START_REF] Garrison | Haplotype-based variant detection from short-read sequencing[END_REF] and vcftools [START_REF] Danecek | The variant call format and VCFtools[END_REF] by enforcing diploidy, retaining positions with quality criterion >40, two alleles, and excluding indels. The MetaT samples were analyzed considering several coverage thresholds, representing a total of 18 datasets (Table 1). For each dataset, no missing value was allowed. They consisted of 6 main datasets obtained after filtering metaT samples based on A.minutum read coverage thresholds and SNPs based on site specific coverage thresholds. Each was analyzed in three ways: 1. by restricting the analyses to the SNPs previously identified using the strains, 2. by restricting the analyses after pruning (see below) SNPs previously identified using the strains

F o r R e v i e w O n l y
9 based on linkage disequilibrium 2. by using SNPs identified de novo from the metaT samples. Table 1 summarizes these datasets, indicating the minimum number of reads that must cover a given SNP (SNP coverage), the minimum number of reads aligning to A. minutum reference transcriptome (Min. reads), the number of metaT samples considered (Samples). The number of strain validated SNPs (Reference SNPs), the number strain validated SNPs after pruning (Pruned ref. SNPs), the number of SNPs identified directly from the metaT datasets (De novo SNPs) and the number of SNPs displaying one allele in strains from a given clade (NE_A, NE_B, Vigo) and the alternative allele in the two other clades (Diagnostic SNPs). In the Pooled_5P6 dataset, the metaT samples from the 5P6 dataset were pooled per geographic site. Only SNPs displaying a minimal allelic frequency > 0.05 in one sample were considered.

Population genomic analyses

5.1.

Nucleotide divergence between strains Nucleotide divergence among strains was calculated as the proportion of variable sites divided by the number of sites, only considering sites covered more than 10 times in all the strains (total of 15,103,704 sites). Strains were clustered using the function "hclust" as implemented in R. For each SNP, the allelic frequencies in the cultivated population of strains was calculated by dividing the number of strains with the reference allele by the total number of strains (Supplementary Figure 1, Strains).

Genetic differentiation

For each simulated population and metaT dataset, SNP coverage and reference allele read counts were extracted from VCF files. For each SNP, observed allelic frequencies were calculated by dividing the reference allele read counts by the coverage (Supplementary Figure 1, Simulated and metaT). The anova method implemented in the R package poolfstat [START_REF] Hivert | Measuring Genetic Differentiation from Pool-seq Data[END_REF], with poolsize parameter set to 10,000, was used to compute pairwise Fst estimates: 1. To investigate the precision of Fst estimates and how it is affected by the level

F o r R e v i e w O n l y
10 of coverage, pairwise Fst were calculated among the 10 simulated populations for each overall coverage level (1e+04, 5e+04, 1e+05, 5e+05, 1e+06, 5e+06 and 1e+07 reads) and considering seven SNP coverage levels (All SNPs, SNPs covered by more than 5, 10, 20, 30, 50, and 100 reads). 2. To investigate the potential bias that could result from an estimation of allelic frequencies using RNA and not DNA (for instance due to allele specific differential expression in populations), Fst between the simulated populations and the actual allelic frequencies of the cultivated population of strains (see above) was calculated at the various overall and SNP specific coverage levels (see above). 3. To determine the potential bias that may result from using several sequencing approaches, five metaT samples were sequenced using both 2x100bp and 2x150bp approaches. Fst was calculated within each of the five sample pairs considering SNPs covered by more than 5, 10, 20, 30, 50, 100 reads. 4. To determine the structure of the natural A. minutum populations, pairwise Fst were calculated between metaT samples for each dataset (Table 1). For points 1, 2 and 3 above, Fst=0 are expected in case of absolute precision and total absence of bias.

5.3.

Comparing genetic structure inferred using strains and metaT data The pattern of genetic variability inferred from the 37 strains and 77 metaT samples was explored using a PCA approach, starting from the vcf files, as implemented in PLINK [START_REF] Purcell | PLINK : A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses[END_REF]. Linkage disequilibrium pruning was performed on the strain dataset, using 50kb windows (meaning that each transcript is analyzed as a whole), a 10bp window step size, keeping SNPs displaying a R² < 0.1.

5.4.

Obtaining a folded joint allele frequency spectrum (JAFS) from metaT From the three versions of the 5P6 dataset (strain validated SNPs, pruned and de novo SNPs, Table 1), for each SNP, reference allele count and coverage were summed per geographic site, leading to a single pooled sample per site (pooled 5P6 dataset, Table 1). For all the SNPs

F o r R e v i e w O n l y
11 covered more than 30 times in each pooled sample, the rare allele frequency was computed. SNPs with a minimal allele frequency systematically < 0.05 per site were discarded. The distribution of the rare allele frequency was calculated per site and the folded joint allele frequency spectrum plotted for each of the three pairs of populations using the package hexbin in R.

Investigating coding regions genome wide divergence

The genetic linkage between contigs was established following an A. minutum linkage map previously developed (Mary et al., Submitted). Using the pooled 5P6 datasets, for each SNP in each geographic site, haplotype diversity was calculated as where p i is the 1 -∑ 𝑝 2 𝑖 frequency of each of the two alleles. The R package pcadapt with "pool" option and k=2 was used to identify SNPs displaying extreme allele frequency differences between populations [START_REF] Luu | pcadapt : An R package to perform genome scans for selection based on principal component analysis[END_REF]. The rollapply function from the zoo package implemented in R was used to determine moving average in terms of haplotype diversity and number of significant SNPs (from pcadapt) along the linkage groups. Values below the 0.5 th and above the 99.5 th percentiles were identified.

Results

Nucleotide divergence between cultivated strains

Nucleotide divergence between the 37 monoclonal A. minutum strains was investigated at 227,829 SNP positions using mRNA sequences. Strain clustering based on the nucleotide divergence indicated the occurrence of three diverging clusters (Figure 2). The first one, hereafter named NE_A, was composed of 27 strains isolated from the Bay of Brest, Penzé and Rance Estuary in France, as well as from Cork harbor in Ireland. Strains from this cluster were isolated during algal blooms between 1989 and 2013, as well as after germination of

F o r R e v i e w O n l y
12 cysts preserved in sediment cores dated from 1947 to 2006. The second one, hereafter named NE_B, was composed of three strains isolated outside of blooming periods in 2010 and 2011 in the Bay of Concarneau and Brest. The third one was composed of seven strains isolated from the Bay of Vigo during a red tide in 2018 (Supplementary Table 1).

Fst precision and bias based on RNAseq data

The impact of coverage levels on the precision of allele frequency estimation resulting from population wide mRNA sequencing was quantified by pooling and subsampling the strain mRNA sequences at different coverage levels. Fst were quantified between replicate simulated populations. This was done separately for PE100 and PE150 datasets (Supplementary Table 1). Fst precision considerably improved with coverage levels. Below a coverage of 1e+05 reads, Fst estimation was imprecise, with Fst estimate between replicate simulated populations higher than 0.2. At a coverage of 5e+05, precision dropped to ~0.1, especially if the analysis was restricted to SNPs with a coverage > 5. At a coverage level > 5e+06, Fst precision was around 0.01 (Figure 3A andB). To quantify the potential bias of inferring allele frequencies from mRNA and not DNA (for instance due to potential allele specific expression patterns), allele frequencies from the simulated populations were compared to the ones calculated following independent genotyping of the strains. Fst bias was systematically lower than 0.05 and often as low as 0.015 (Figure 3C andD). To determine whether read length may influence Fst estimates, five metaT samples were sequenced using both 2x100bp and 2x150bp reads (Supplementary Figure 2). For each sample, Fst values between 2x100bp and 2x150bp read datasets were at time higher than 0.1. This indicated that the sequencing strategy may moderately influence Fst estimates. However, when calculated based on highly covered SNPs (>20), Fst were always below 0.05, indicating that such a potential issue may be solved by focusing on highly covered SNPs. The genetic variability determined using strain and metaT samples was compared using a PCA (Figure 4). The first axis separated strains from the three clades identified above (Figure 4A; NE_A, NE_B and Vigo). The second and third axes separated the strains belonging to the NE_B clade and highlighted the high genetic variability existing within this clade (Figure 4A,B). All the metaT samples from the Bay of Brest were grouped with NE_A strains, indicating that they are composed of NE_A cells. All the metaT samples from the Bay of Vigo were grouped with the Vigo strains, indicating that they are composed of cells belonging to the Vigo clade. The only metaT sample from the Penzé estuary was relatively close to the Vigo samples (strains and metaT) along the first axis, and close to the NE_A strains and Bay of Brest metaT samples along the third axis. It is worth mentioning that some strains from the NE_A clade were isolated from the Penzé estuary from 1989 to 2010 (Figure 2), while the Penzé metaT sample was sampled in 2015 (Supplementary Table 2). The Penzé estuary metaT sample may either be composed of cells belonging to a fourth distinct population or result from admixture from two or three of the previously identified populations. To investigate this, allelic frequencies at diagnostic SNP positions (SNPs displaying one allele in all strains from a given clade (NE_A, NE_B, Vigo) and the alternative allele in the two other clades) were analyzed in all metaT samples (Supplementary Figure 3). For a great majority of these SNPs, the Penzé metaT sample displayed a fixed allele. Depending on the SNP, the allele corresponded to the one identified in NE_A, Vigo, and more rarely NE_B clade. The absence of intermediate allelic frequencies indicated a new population, from which no strain had been isolated, and not an admixed one. genetic structure between samples within each site Pairwise Fst were calculated between all metaT samples and results were summarized in Figure 5. Within each site, all Fst values were below 0.03, i.e. about at the precision limit of the method, indicating an absence of intra-site genetic differentiation. It should be noted that the 36 metaT datasets from the Bay of Brest were sampled during A. minutum blooms that occurred over three consecutive years (Supplementary Table 2), indicating both intra-and inter-annual stability of the population genetic composition. In sharp contrast, genetic differentiation between geographic sites was extremely strong, with median Fst values of 0.55, 0.59, and 0.71 between Bays of Brest and Vigo, Bay of Brest and Penzé Estuary, and Bay of Vigo and Penzé Estuary, respectively.

Folded joint allele frequency spectrum (JAFS) is compatible with ancient divergence without gene flow

The metaT samples coming from the same geographic site were pooled and allelic frequencies at the three geographic sites compared using folded JAFS (Figure 6). The three JAFS highlighted extreme allele frequency differences in the three populations. Most of the SNPs were distributed along the axes, indicating that alleles segregating at intermediate frequencies at a given site are often absent of the two other sites (private alleles). Moreover, at numerous SNP positions, alleles appearing fixed, or almost fixed, at a given site were totally absent from the two other sites. The JAFS also showed a strong excess of SNP positions displaying alleles restricted to the Bay of Brest, indicating that this population was genetically more diverse than the two others. 6. Gene specific structure may be inferred from metaT samples For each population, haplotype diversity was investigated along A. minutum linkage groups (Figure 7A). As already identified using the JAFS (Figure 6), diversity was higher in the population from the Bay of Brest. Haplotype diversity fluctuated along the linkage groups, with transient increase or decrease. Haplotype diversity modifications were population specific rather than shared between the three populations. In agreement with the JAFS analysis, 10,099 SNPs were identified as displaying different allele frequencies in the three populations (pcadapt, adjusted p-value < 1e-10). These SNPs were widespread along the linkage groups, but a few genomic regions displayed higher proportions. Genomic regions enriched in significant SNPs often corresponded to regions displaying low haplotype diversity in one or two populations, but this was not systematically the case. Two genomic regions appeared of special interest in linkage groups L1 and L37 (Figure 7B,C). From 0 to 1.6 cM of L1, haplotype diversity was extremely low in the population from the Bay of Brest and displayed a slight increase in SNPs identified as displaying different allelic frequencies in the three populations (Figure 7B). This region contained 96 SNPs coming from 26 transcripts, including 12 annotated ones (Table 2). Haplotype diversity was extremely low in the Bay of Vigo population at 107.4 cM in L37 (Figure 7C). This region encompassed 327 SNPs in 43 transcripts, including 15 annotated (Table 3) and was characterized by a strong excess of SNPs displaying different allelic frequencies in the three populations.

7. Genetic structure may be inferred from metaT samples even in absence of pre-existing strain validated SNP database Results presented above in terms of Fst, JAFS and genome wide divergence analyses from metaT datasets were restricted to strain validated SNPs. These same analyses were also performed using SNPs identified de novo using the metaT datasets. Overall results were extremely similar using strain validated and de novo SNPs, with a strong structure between

F o r R e v i e w
O n l y geographic sites and absence of structure within sites (Supplementary Figure 4), with the occurrence of numerous private alleles and of fixed polymorphism (Supplementary Figure 5). Fluctuating haplotype diversity along linkage groups and low diversity in the Bay of Brest at the beginning of L1 and in the Bay of Vigo at the end of L37 (Supplementary Figure 6) was also detected. A few differences may nevertheless be noticed. First, the number of SNPs considered was several times higher when using de novo SNPs (Table 1). Second, interpopulation Fst values were lower (0.10<Fst<0.20 with de novo SNPs compared to 0.5<Fst<0.7 with strain validated SNPs). Third, a high proportion of SNPs displayed similar allelic frequencies in the three populations (Supplementary Figure 5). Finally, haplotype diversity was higher in the Bay of Vigo and Penzé Estuary populations when using de novo SNPs (but still slightly lower than in the Bay of Brest).

Discussion

In the present study, the genetic structure of the dinoflagellate A. minutum populations in Western Europe was investigated using strains and metaT samples. MetaT datasets were shown to be extremely insightful to decipher the complex genetic structure of such a microbial organism. Results highlighted very strong genetic structure probably resulting from an ancient divergence without gene flow, numerous SNP markers displaying private alleles spread-out in the genomes in the different populations, as well as a few genomic regions displaying very low genetic diversity in one population. These results are discussed below.

Using metaT samples to infer protist genetic structure

The estimation of allelic frequencies is not biased by gene expression. This could have been a major issue in case of widespread and systematic allele specific expression. To illustrate the potential issue, one may imagine two populations living in two distinct environments, A and B, each displaying two alleles, x and y, at a 50/50 allele frequency. However, if allele x is more expressed in environment A and allele y in environment B, estimated allele frequency from F o r R e v i e w O n l y mRNA sequences would be biased compared to actual allelic frequencies. Allele specific expression (ASE) is mostly studied in humans, by monitoring the relative expression of the two gene copies in heterozygotes across cell types and tissues. Results suggest unequal expression of gene copies may be widespread, but mostly for gene copies displaying genetic variation in cis-regulatory regions. Moreover the observed bias tends to vary from cell to cell or tissue to tissue [START_REF] Cleary | Perspectives on Allele-Specific Expression[END_REF][START_REF] Montgomery | Rare and Common Regulatory Variation in Population-Scale Sequenced Human Genomes[END_REF][START_REF] Wagner | Computational Analysis of Whole-Genome Differential Allelic Expression Data in Human[END_REF]. Our pooling and resampling approach clearly shows that population wide ASE is not a global issue. Nevertheless, we cannot rule out that for specific SNP markers, especially in strong linkage disequilibrium with a cis-regulatory variant, the observed allelic frequency inferred from metaT datasets would be influenced by population wide ASE. In metaT samples, SNP coverage is difficult to control a priori, because it strongly depends upon the relative frequency of the species of interest in the sampled community, as well as upon the relative expression of the carrying gene relative to all other genes expressed by this species. However, it is a critical factor to properly estimate allele frequencies. Indeed, as SNP coverage decreases, observed allelic frequency would be strongly affected by sampling error. The pooling and resampling approach helped identify the relationship between species coverage and expected Fst precision. For A. minutum, the resampling approach suggested that a minimum coverage of 5e+05 reads is sufficient to detect genetic structure corresponding to Fst > 0.1 and coverage higher than 5e+06 reads enabled extremely precise Fst estimates (Fst ~ 0.01). The main consequence is that, depending on the sample sequencing depth, a high relative abundance of the species of interest in sampled communities may be required to accurately identify very low levels of genetic differentiation (for instance, Fst ~ 0.01 may be detected with a relative abundance of 0.25 if the sample sequencing depth is 2e+07, but only 0.05 if the sample sequencing depth is 1e+08). However, high levels of differentiation (Fst > 0.1) could be identified, even in very moderately abundant species (for instance a relative abundance of 0.025 or 0.005 if the sample sequencing depth is 2e+07 or 1e+08, respectively). As a matter of comparison, the average sequencing depth in the samples used in the present

F o r R e v i e w O n l y
18 study was 2e+07 [START_REF] Metegnier | Species specific gene expression dynamics during harmful algal blooms[END_REF] and 1.6e+08 in Tara Ocean samples [START_REF] Carradec | A global ocean atlas of eukaryotic genes[END_REF]. The comparison between strain and metaT samples was extremely insightful. Population structure inferred from the two approaches was extremely similar. This has several implications. At the same time, it encourages the use of metaT samples to determine genetic structure of protist populations, but also confirms that strain based population samples are not necessarily biased compared to natural populations. This could have been the case if the number of strains was too small to capture the genetic diversity of natural populations, or if single cell isolation and subsequent culturing steps to obtain clonal strains selected an unrepresentative subset of natural populations. We should note that, in the present system, genetic structure was extremely strong, with numerous private alleles, making it easy to detect even with a small strain sample size. Another insightful result of the metaT sample analyses was that replicate samples from the same population displayed very low temporal variability, indicating that the population genetic diversity may be captured using a very low number of samples. As an illustration, we may refer to the extreme similarity of the samples from the Bay of Vigo, sampled over a distance of ten kilometers, or of the samples from the Bay of Brest obtained at the same locality during three consecutive years. This has profound implications for future investigation of protist population structure. Indeed, using strains, a strong sampling effort has to be taken to capture the genetic diversity of the population of interest at a given locality, limiting the feasibility of genetic structure analyses to a handful of sites. Using metaT samples, as the population genetic diversity may be captured using a very limited number of easy to obtain samples, a finer spatial sampling resolution may be considered. Increasing the spatial coverage with a relatively low resolution can be an essential step towards a better understanding of the extremely complex spatio-temporal protist population genetic structure (see below).

F o r R e v i e w O n l y
The present study highlighted that population genetic analyses can be performed from metaT datasets even in absence of pre-validated SNPs. Indeed, results obtained from strain validated SNPs and SNPs detected de novo from the metaT datasets were extremely similar. This was true for the overall genetic structure, as well as for the genomic regions, genes and even SNPs (identified as displaying different allelic frequencies across populations). We noted that Fst values were much lower using de novo SNPs and that numerous SNPs identified de novo appeared to display similar allelic frequencies across populations, while this was rarely the case when using strain validated SNPs. Such differences may be explained by the peculiar organization of dinoflagellate genes. Indeed, these are organized in tandem repeats [START_REF] Stephens | Genomes of the dinoflagellate Polarella glacialis encode tandemly repeated single-exon genes with adaptive functions[END_REF][START_REF] Wisecaver | Dinoflagellate Genome Evolution[END_REF]. During mapping, reads belonging to different gene copies align to the same reference transcript. However, the different gene copies within a single cell may display SNPs, reflecting the genetic divergence of gene copies following gene duplication in a given genome. As dinoflagellate vegetative cells are haploid, it is possible to exclude these markers by restricting the analyses to SNPs displaying monomorphism within a given genome but polymorphism across genomes when genotyping individual strains (see methods). Using SNPs identified de novo from metaT datasets, such restriction is impossible and the population genetic analyses include both markers of genetic diversity between and within genomes. The observed difference between strain-validated and de novo SNPs (lower Fst and shared polymorphism with de novo SNPs) thus probably resulted from ancestral gene duplications and mutations of the various gene copies that preceded the population splits. Despite the confounding effect of the genetically variable gene copies within a given genome, the signal of genetic divergence between populations is still captured from metaT samples. Using a de novo SNP approach is not only a possibility, but it may also be a better choice than using strain validated SNPs. Indeed, in the present study, the haplotype diversity of the Penzé Estuary populations was much lower using the strain validated SNPs compared to the de novo approach, probably because no strains from this population were isolated and genotyped, preventing the identification of the Penzé Estuary population private polymorphism using the strain validated approach. These results are

F o r R e v i e w
O n l y 20 especially promising for applying metaT based population genetic approaches to a wide range of species, including species difficult to grow in the lab but with a reference transcriptome or genome available, or even species not cultivated but for which reference transcriptomes could be computed directly from metaT or metaG datasets [START_REF] Delmont | Functional repertoire convergence of distantly related eukaryotic plankton lineages revealed by genome-resolved metagenomics[END_REF][START_REF] Vorobev | Transcriptome reconstruction and functional analysis of eukaryotic marine plankton communities via high-throughput metagenomics and metatranscriptomics[END_REF].

Strong divergence between A. minutum populations in Western

Europe Within A. minutum, a strong genetic structure was inferred both from strains and metaT samples in Western Europe. A total of four highly divergent groups were identified. The first group (NE_B) was identified only using strains sampled outside of massive A. minutum developments, i.e. bloom events. The metaT samples were exclusively obtained from blooms and there was no sign of the presence of individuals belonging to the NE_B group in these samples. Strains from NE_B were phenotypically distinct from the ones from the other three populations, as they did not produce PSTs [START_REF] Geffroy | From the sxtA4 Gene to Saxitoxin Production : What Controls the Variability Among Alexandrium minutum and Alexandrium pacificum Strains[END_REF]. Interestingly, the genetic diversity of this population was much higher than the other three populations. Theoretically, genetic diversity is directly related to the effective population size [START_REF] Ellegren | Determinants of genetic diversity[END_REF], which would imply that the NE_B effective population size is much higher than the three others. Considering the hypothesis that the NE_B population is found at low cell density, this might seem counterintuitive. However, dense A. minutum blooms tend to occur transiently, in restricted geographic areas, like the one in Bay of Vigo in 2018, to which the strains used in this study belong [START_REF] Nogueira | HABs in coastal upwelling systems : Insights from an exceptional red tide of the toxigenic dinoflagellate Alexandrium minutum[END_REF]. Moreover, massive development probably involves numerous rounds of mitotic division, contributing to decreasing the effective population size despite huge census population sizes [START_REF] Ellegren | Determinants of genetic diversity[END_REF]. Given the volumes at play at sea, very low density populations across broad areas could potentially outnumber both census and effective population sizes of bloom forming populations. The role of low density

F o r R e v i e w O n l y
21 populations is rarely considered in protist population genetics although we might speculate that they could play an important role in terms of metapopulation dynamics. The other three populations were exclusively sampled during blooms. The first one, NE_A, was sampled thanks to strains isolated in the Atlantic Ocean and English Channel and corresponded to all metaT sampled from the Bay of Brest during three consecutive years. This population was the most genetically diverse bloom forming population and displayed an extremely stable genetic structure in the Bay of Brest, during the three years sampled using metaT, but also probably during several decades, as two strains isolated from sediment cores and dated to 1947 belonged to this same population. The second one corresponded to strains and metaT samples obtained from the Bay of Vigo during a huge red tide. Finally, the last population was sampled from a single metaT sample from the Penzé Estuary in 2015. Private alleles tend to indicate that this last population was truly distinct from the others and did not correspond to the admixture of other populations. Surprisingly, strains isolated from the Penzé Estuary from 1989 to 2013 all belonged to the NE_A clade. With such a limited sampling, it is difficult to determine whether this new population was transient or if it replaced the NE_A population. Altogether, our study highlighted an extremely complex A. minutum population structure, with extremely high levels of divergence supported by the presence of numerous private alleles, including fixed polymorphism. The divergence between NE_A and NE_B had been characterized previously using a few strains (Le Gac et al., 2016c) and the present study showed that several A. minitum populations exhibit extensive divergence in Western Europe. The observed divergence pattern is compatible with ancestral divergence and an absence of contemporary gene flow between populations [START_REF] Almeida | A population genomics insight into the Mediterranean origins of wine yeast domestication[END_REF][START_REF] Ellegren | Genome sequencing and population genomics in non-model organisms[END_REF][START_REF] Gutenkunst | Inferring the Joint Demographic History of Multiple Populations from Multidimensional SNP Frequency Data[END_REF]. Due to their extensive dispersal abilities in the marine environment, protist species are commonly expected to be cosmopolite and display homogeneous populations worldwide (De [START_REF] De Wit | Everything is everywhere, but, the environment selects'; what did Baas Becking and Beijerinck really say[END_REF][START_REF] Finlay | Global Dispersal of Free-Living Microbial Eukaryote Species[END_REF]. However, strong genetic structure is extremely common in marine protist species [START_REF] Casabianca | Population genetic structure and connectivity of the harmful dinoflagellate Alexandrium minutum in the Mediterranean Sea[END_REF][START_REF] Craig | Patterns of population structure and complex haplotype sharing among field isolates of the green alga Chlamydomonas reinhardtii[END_REF][START_REF] Gao | Spatiotemporal genetic structure of regional-scale Alexandrium catenella dinoflagellate blooms explained by extensive dispersal and environmental selection[END_REF][START_REF] Godhe | The role of intraspecific variation in the ecological and evolutionary success of diatoms in changing environments[END_REF][START_REF] Paredes | Population Genetic Structure at the Northern Edge of the Distribution of Alexandrium catenella in the Patagonian Fjords and Its Expansion Along the Open Pacific Ocean Coast[END_REF][START_REF] Rengefors | Genetic diversity and evolution in eukaryotic phytoplankton : Revelations from population genetic studies[END_REF]. This pattern is in sharp contrast with marine macroorganism populations often displaying extremely low levels of divergence between them [START_REF] Gagnaire | Using neutral, selected, and hitchhiker loci to assess connectivity of marine populations in the genomic era[END_REF][START_REF] Waples | Separating the wheat from the chaff : Patterns of genetic differentiation in high gene flow species[END_REF]. Protist genetic structure is sometimes correlated with geographic or hydrodynamic features [START_REF] Casabianca | Population genetic structure and connectivity of the harmful dinoflagellate Alexandrium minutum in the Mediterranean Sea[END_REF][START_REF] Casteleyn | Limits to gene flow in a cosmopolitan marine planktonic diatom[END_REF][START_REF] Godhe | Physical barriers and environmental gradients cause spatial and temporal genetic differentiation of an extensive algal bloom[END_REF][START_REF] Postel | Adaptive divergence across Southern Ocean gradients in the pelagic diatom Fragilariopsis kerguelensis[END_REF] as well as with various environmental parameters [START_REF] Sassenhagen | Comparison of Spatial and Temporal Genetic Differentiation in a Harmful Dinoflagellate Species Emphasizes Impact of Local Processes[END_REF][START_REF] Sjöqvist | Local adaptation and oceanographic connectivity patterns explain genetic differentiation of a marine diatom across the North Sea-Baltic Sea salinity gradient[END_REF][START_REF] Whittaker | Evidence for environmental and ecological selection in a microbe with no geographic limits to gene flow[END_REF]. Moreover divergent populations sometimes co-occur for extensive periods of time [START_REF] Härnström | Hundred years of genetic structure in a sediment revived diatom population[END_REF][START_REF] Lundholm | Exploring the impact of multidecadal environmental changes on the population genetic structure of a marine primary producer[END_REF]. The level of divergence identified in the present study and in other species is on the order of what is expected for sibling species. One possibility to explain the observed level of divergence would be the existence of numerous cryptic species, evolving without exchanging genes, but displaying virtually no morphological difference, within taxonomically recognized protist species. Hence numerous cryptic species have been repeatedly identified from what were previously thought to be single cosmopolitan protist species [START_REF] De Luca | Resolving cryptic species complexes in marine protists : Phylogenetic haplotype networks meet global DNA metabarcoding datasets[END_REF][START_REF] Smayda | Cryptic planktonic diatom challenges phytoplankton ecologists[END_REF]. Even using molecular tools, identifying cryptic species may be especially difficult for protists. Indeed, the coalescence time (the time since the most common ancestor of two gene copies) is expected to be directly related to the effective population size [START_REF] Hudson | Gene genealogies and the coalescent process[END_REF]. With huge population sizes, the coalescence time may easily be much higher than the speciation time. As a consequence, ancestral polymorphism may remain in contemporary cryptic species, blurring the phylogenetic signal of species split and preventing the identification of cryptic species [START_REF] Rannala | The Multi-species Coalescent Model and Species Tree Inference[END_REF]. Given their census population size, intra-population genetic diversity would be expected to be extremely high. However, the present study and previous ones [START_REF] Blanc-Mathieu | Population genomics of picophytoplankton unveils novel chromosome hypervariability[END_REF][START_REF] Filatov | Extreme Lewontin's Paradox in Ubiquitous Marine Phytoplankton Species[END_REF], indicated that the observed levels of diversity in protists tend to be on the order of what was observed for moderately diverse animal species displaying much lower census [START_REF] Romiguier | Comparative population genomics in animals uncovers the determinants of genetic diversity[END_REF]. Previous studies calculated that given the observed genetic diversity, the effective population size should be on the order of tens of millions of cells, a number of cells that may be found in a few liters of sea water [START_REF] Filatov | Extreme Lewontin's Paradox in Ubiquitous Marine Phytoplankton Species[END_REF]. Such discrepancy between the census population size and the level of diversity is known as the Lewontin's paradox [START_REF] Buffalo | Quantifying the relationship between genetic diversity and population size suggests natural selection cannot explain Lewontin's Paradox[END_REF][START_REF] Ellegren | Determinants of genetic diversity[END_REF][START_REF] Filatov | Extreme Lewontin's Paradox in Ubiquitous Marine Phytoplankton Species[END_REF][START_REF] Galtier | How Much Does Ne Vary Among Species?[END_REF]. Three main solutions have been proposed to solve this paradox. The first one emphasizes the possibility of reduced mutation rate in species with large populations. Mutation rate estimation are rather limited in protist species, but tend to indicate that it does not deviate from rates of species displaying much lower census population sizes [START_REF] Krasovec | Evolution of Mutation Rate in Astronomically Large Phytoplankton Populations[END_REF]. The second one emphasizes the importance of demographic fluctuations and asexuality on reducing the effective population size. The third one, the reduction of genetic diversity at the neighborhood of sites increasing in frequency in a population due to selection (genetic hitchhiking or draft). A thorough investigation in the haptophyte Emiliania huxleyi appeared compatible with a major importance of recurrent selective events in shaping the genetic diversity of this species [START_REF] Filatov | Extreme Lewontin's Paradox in Ubiquitous Marine Phytoplankton Species[END_REF]. In the present study, the numerous SNPs displaying extremely divergent allelic frequencies could result from adaptive evolution occurring independently in the different populations. These SNPs appeared spread-out in the genome and most of the time they were not surrounded by genomic regions of low haplotype diversity. Such a signature would be compatible with selection from standing genetic variation during which neutral alleles segregating in the population become adaptive following a modification of the selective pressures [START_REF] Barrett | Adaptation from standing genetic variation[END_REF]. As they appear in various genomic contexts, such mutations could increase in frequency in the population with minimal effect on the genetic diversity at neighboring sites. Such events would correspond to soft selective sweeps [START_REF] Ellegren | Genome sequencing and population genomics in non-model organisms[END_REF]. Nevertheless, a few genomic regions tend to display a very low level of diversity in one population compared to the others, but also compared to the other genomic regions of the same population. Such events are compatible with hard selective sweeps during which a de novo adaptive mutation appears in a population and rises to fixation with the genomic context in which it appeared [START_REF] Ellegren | Genome sequencing and population genomics in non-model organisms[END_REF].

F o r R e v i e w O n l y 22 Wit &
F o r R e v i e w O n l y population sizes
F o r R e v i e w O n l y 24
As a conclusion metaT datasets may be great to investigate the population genetics of protist species relatively abundant in given communities. By circumventing the time-consuming culturing steps required to genotype individual cells from natural populations, metaT but also any other population or community wide dataset (metaG, multiplex amplicon sequencing…) should be of primary importance to improve the spatio-temporal sampling of protist populations. Such improvements are mandatory if we want to better understand the population genetics of these species that display moderate diversity and complex spatio-temporal genetic structure despite theoretical expectations of extremely high diversity and low genetic structure. One of the major constraints is that most of the population genomic tools (linkage disequilibrium analyses, simulations…) were developed considering individual genotypes and not allelic frequencies estimated from population wide samples. A situation that may change given the growing interest in using eDNA/RNA to determine the population genetics of a wide range of organisms [START_REF] Sigsgaard | Population-level inferences from environmental DNA-Current status and future perspectives[END_REF] 1).The three strain populations NE_A, NE_B and Vigo are indicated. Figure 3: Precision and bias of Fst estimates based on mRNA sequences. A and B correspond to pairwise Fst among pairs of replicate simulated populations based on PE100 (A) and PE150 (B) strain samples for various coverage levels. C and D correspond to Fst between the population of strains (based on individual strain genotypes) and each of the ten simulated populations, based on PE100 (A) and PE150 (B) strain samples for various coverage levels. R indicates the number of reads subsampled for the simulated populations, C is coverage level used to filter out SNPs (All indicates no filtering, all SNPs are considered), S is the number of SNP analyzed. Figure 4: Population structure for strain and metaT samples. Genetic variation along four axes following PCA analysis. The results correspond to the pruned 5P6 dataset, but similar results were obtained for all the other datasets (not shown). PC1 and 2 are represented on A and PC3 and PC4 on B. The inset indicates the percentage of variance explained by each principal component. Colors indicate the geographic origin of the strains and symbols indicate strains (circles) and metaT datasets (triangles). Figure 5: Boxplot indicating the pairwise Fst between and within populations based on metaT samples. The results correspond to the 5P6 dataset, but similar results were obtained for all the other datasets (not shown). Figure 6: Folded Joint Allele Frequency Spectrum (JAFS) comparing the rare allele frequencies in the three populations using metaT samples. Colors indicate the number of SNPs falling in each bin defined by a unique combination of allele frequency in the two populations considered. For each of the three population pairs, The results correspond to the 5P6 dataset, but similar results were obtained for all the other datasets (not shown). Figure 7: Genome wide diversity inferred from metatT. A. Haplotype diversity (H) moving average (width=100SNPs, steps=20 SNPs) along A. minutum linkage groups in the Bay of Brest (purple), Bay of Vigo (red) and Penze Estuary (green) populations. The black line represents the moving average proportion of SNPs displaying different allelic frequencies across the three populations (adjusted p-values < 1e-10). Dots indicate regions falling above or below of the 95th and 5th quartiles, respectively. Background colors correspond to the various linkage groups. B. and C. Haplotype diversity (purple, red and green, color code as in A.) and proportion of significant SNPs (black broken lines) moving average (width=50SNPs, steps=10 SNPs) in linkage groups L1 and L37, respectively. The number of SNPs in each linkage group is indicated

F o r R e v i
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Prop. Signif.

  structure inferred from strains and metaT samples are similar

  genetic structure between the Bay of Brest and Vigo, but no

Table 3 :

 3 . Benefits from this research accrue from the sharing of our data on public databases as described above. The data that supports the findings of this study are available in the supplementary material of this article, as well as from the public database cited above. Annotated genes in Linkage group L37 at 107.4 cM Geographic origin of the strains and metaT samples. Dots and stars indicate the geographic origin of metaT samples and strains, respectively. Several strains and metaT samples may have the same geographic origin, see Supplementary Table1 and 2. Strain clustering. Strains are clustered based on nucleotide divergence. colored dots indicate the geographic origin of the strains. Isolation date for strains isolated from water samples or sediment core datations for strains isolated from sediment cores (Supplementary Table

	F o
	r
	R
	e v
	i e
	w
	O n
	l y

Acknowledgements

The study was funded by PRIMROSE (EC Interreg Atlantic Area EAPA 182 /2016) and the Brittany Region as part of the Paleoecology of Alexandrium minutum dans la Rade de Brest-Marché n°2017-90292 project PALMIRA. We thank all the participants in the crew of the RV Ramón Margalef from the Remedios cruise (Research project -grant number CTM2016-75451-C2-1-R), particularly B. Mourino-Carballido, for their support to the sample collection. We thank the Ifremer Sebimer team for bioinformatic support as well as the INRAE GeT-PlaGe sequencing platform for sequencing.

References:

Author Contributions

The research was designed by MLG, CD, and MS. The sampling was performed by MLG, MS, GM, JQ, RS, FR. The molecular biology was performed by MLG, JQ, GM. The bioinformatics and population genomic analyses were done by MLG, GM, LM. Writing of the article was carried out by MLG, CD, LM, MS, RS, FR.

Tables Table 1: Summary of the metaT datasets considered. 4 : Boxplot indicating the pairwise Fst between and within populations based on metaT samples using de novo SNPs. The results correspond to the 5P6 dataset, but similar results were obtained for all the other datasets (not shown). Supplementary Figure 5 : Folded Joint Allele Frequency Spectrum (JAFS) comparing the rare allele frequencies in the three populations using de novo SNPs identified from metaT samples.

Colors indicate the number of SNPs falling in each bin defined by a unique combination of allele frequency in the two populations considered. For each of the three population pairs, The results correspond to the 5P6 dataset, but similar results were obtained for all the other datasets (not shown). Supplementary Figure 6: Genome wide diversity using de novo SNPs identified from metaT samples. A. Haplotype diversity (H) moving average (width=100SNPs, steps=20 SNPs) along A. minutum linkage groups in the Bay of Brest (purple), Bay of Vigo (red) and Penze Estuary (green) populations. The black line represents the moving average proportion of SNPs displaying different allelic frequencies across the three populations (adjusted p-values < 0.05). Dots indicate regions falling above or below of the 95th and 5th quartiles, respectively. Background colors correspond to the various linkage groups. B. and C. Haplotype diversity (purple, red and green, color code as in A.) and proportion of significant SNPs (black broken lines) moving average (width=50SNPs, steps=10 SNPs) in linkage groups L1 and L37, respectively. The number of SNPs in each linkage group is indicated