Hugo Doré

Gregory K Farrant

Ulysse Guyet

Julie Haguait

Florian Humily

Morgane Ratin

Frances D Pitt

Martin Ostrowski

Christophe Six

Loraine Brillet-Guéguen

Mark Hoebeke

Antoine Bisch

Gildas Le Corguillé

Erwan Corre

Karine Labadie

Jean- Marc Aury

Patrick Wincker

Dong Han Choi

Jae Hoon Noh

Damien Eveillard

David J Scanlan

Frédéric Partensky

Laurence Garczarek
email: laurence.garczarek@sb-roscoff.fr

Evolutionary Mechanisms of Long-Term Genome Diversification Associated With Niche Partitioning in Marine Picocyanobacteria

Keywords: Marine cyanobacteria 1, Prochlorococcus2, Synechococcus 3, comparative genomics 4, niche adaptation 5, amino-acid substitutions 6, genomic islands 7, evolution 8 . Genome diversification in marine picocyanobacteria

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Understanding how phytoplankton species have adapted to the marine environment, a dynamic system through time and space, is a significant challenge, notably in the context of rapid global change [START_REF] Edwards | Ferredoxin and flavodoxin as biochemical indicators of iron limitation during open-ocean iron enrichment[END_REF][START_REF] Sears | Introduction to the symposium: Responses of organisms to climate change: A synthetic approach to the role of thermal adaptation[END_REF][START_REF] Irwin | Phytoplankton adapt to changing ocean environments[END_REF][START_REF] Doblin | Drift in ocean currents impacts intergenerational microbial exposure to temperature[END_REF]. Even though these microorganisms might adapt more rapidly than larger organisms to environmental change due to their short generation times and large population sizes, the underlying mechanisms and timescales required for such evolutionary processes to occur remain mostly unknown. One of the best ways to better understand these processes is by deciphering the links between current genomic diversity and niche occupancy of these organisms. Such an approach requires complete genomes with representatives of distinct ecological niches, a resource which remains limited even with the advent of high-throughput sequencing and the multiplication of partial single amplified genomes (SAGs; [START_REF] Stepanauskas | Matching phylogeny and metabolism in the uncultured marine bacteria, one cell at a time[END_REF][START_REF] Malmstrom | Ecology of uncultured Prochlorococcus clades revealed through single-cell genomics and biogeographic analysis[END_REF][START_REF] Kashtan | Single-cell genomics reveals hundreds of coexisting subpopulations in wild Prochlorococcus[END_REF][START_REF] Berube | Emergence of trait variability through the lens of nitrogen assimilation in Prochlorococcus[END_REF][START_REF] Moore | Single-cell genomics unveiled a cryptic cyanobacterial lineage with a worldwide distribution hidden by a dinoflagellate host[END_REF] or metagenomes assembled genomes (MAGs; [START_REF] Iverson | Untangling genomes from metagenomes: Revealing an uncultured class of marine Euryarchaeota[END_REF][START_REF] Haro-Moreno | Fine metagenomic profile of the Mediterranean stratified and mixed water columns revealed by assembly and recruitment[END_REF]. Due to their ubiquity, their natural abundance in situ, the occurrence of welldefined ecotypes and good knowledge of how environmental parameters influence their biogeography, marine picocyanobacteria constitute excellent model organisms to tackle evolutionary processes involved in niche partitioning.

Synechococcus and Prochlorococcus are the two most abundant photosynthetic organisms on Earth (Partensky et al., 1999a;[START_REF] Scanlan | Ecology of Cyanobacteria II: Their Diversity in Space and Time[END_REF]. As major primary producers, they have a pivotal role in CO2 fixation and carbon export and are key players in marine trophic networks [START_REF] Jardillier | Significant CO2 fixation by small prymnesiophytes in the subtropical and tropical northeast Atlantic Ocean[END_REF][START_REF] Flombaum | Present and future global distributions of the marine Cyanobacteria Prochlorococcus and Synechococcus[END_REF][START_REF] Grébert | Light color acclimation is a key process in the global ocean distribution of Synechococcus cyanobacteria[END_REF]. Although these organisms often co-occur in (sub)tropical and temperate waters, Synechococcus is present from the equator to sub-polar waters, while the Genome diversification in marine picocyanobacteria distribution of Prochlorococcus is restricted to the latitudinal band between 45°N and 40°S [START_REF] Johnson | Niche partitioning among Prochlorococcus ecotypes along ocean-scale environmental gradients[END_REF][START_REF] Flombaum | Present and future global distributions of the marine Cyanobacteria Prochlorococcus and Synechococcus[END_REF][START_REF] Paulsen | Synechococcus in the Atlantic Gateway to the Arctic Ocean[END_REF]. This broad distribution implies that these two microorganisms are able to survive in a large range of environmental niches along in situ gradients of temperature, light intensity as well as micro-and macro-nutrients [START_REF] Bouman | Oceanographic basis of the global surface distribution of Prochlorococcus ecotypes[END_REF][START_REF] Zwirglmaier | Global phylogeography of marine Synechococcus and Prochlorococcus reveals a distinct partitioning of lineages among oceanic biomes[END_REF][START_REF] Scanlan | Ecology of Cyanobacteria II: Their Diversity in Space and Time[END_REF][START_REF] Sohm | Cooccurring Synechococcus ecotypes occupy four major oceanic regimes defined by temperature, macronutrients and iron[END_REF][START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF].

The ability of marine picocyanobacteria to occupy various niches is likely related to the high intrinsic genetic diversity of these taxa. The Synechococcus/Cyanobium radiation has been split into three main groups, called Sub-Clusters (hereafter SC) 5.1 to 5.3 [START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF][START_REF] Huang | Novel lineages of Prochlorococcus and Synechococcus in the global oceans[END_REF]. While members of SC 5.2, currently encompassing strains assigned to both the Synechococcus and Cyanobium genera, are restricted to near coastal and estuarine areas, SC 5.1 and 5.3 are mainly marine, with SC 5.1 dominating in most oceanic waters and showing the highest genetic diversity currently comprising 18 distinct clades and 40 sub-clades so far described (Ahlgren and Rocap, 2012a;[START_REF] Mazard | Multi-locus sequence analysis, taxonomic resolution and biogeography of marine Synechococcus[END_REF]. The Prochlorococcus genus forms a branch at the base of the Synechococcus SC 5.1 radiation and although it includes seven major lineages, usually referred to as 'clades', the whole genus is actually equivalent to a single marine Synechococcus clade from a phylogenetic viewpoint [START_REF] Huang | Novel lineages of Prochlorococcus and Synechococcus in the global oceans[END_REF]Biller et al., 2014b;[START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF]. Lineages thriving in the upper mixed layer, socalled High Light-adapted (HL) clades, are genetically distinct from those occupying the bottom of the euphotic zone, so-called Low Light-adapted (LL) clades. Furthermore, while members of HLI were shown to colonize subtropical and temperate waters, HLII to IV are adapted to higher temperatures [START_REF] Johnson | Niche partitioning among Prochlorococcus ecotypes along ocean-scale environmental gradients[END_REF][START_REF] Zinser | Influence Genome diversification in marine picocyanobacteria of light and temperature on Prochlorococcus ecotype distributions in the Atlantic Ocean[END_REF]Martiny et al., 2009b), with HLII colonizing N-poor areas and HLIII and IV being restricted to iron(Fe)-limited environments [START_REF] Rusch | Characterization of Prochlorococcus clades from iron-depleted oceanic regions[END_REF][START_REF] West | A novel clade of Prochlorococcus found in high nutrient low chlorophyll waters in the South and Equatorial Pacific Ocean[END_REF][START_REF] Malmstrom | Ecology of uncultured Prochlorococcus clades revealed through single-cell genomics and biogeographic analysis[END_REF]. For Synechococcus, distribution and environmental preferences have only been well characterized for the five dominant clades in the field (clades I to IV and CRD1). Members of clades I and IV have been shown to be cold thermotypes that dominate in coastal, mixed and/or high latitude, nutrient-rich waters, while clades II and III are warm thermotypes, predominating in N-depleted areas and P-depleted regions, respectively [START_REF] Zwirglmaier | Global phylogeography of marine Synechococcus and Prochlorococcus reveals a distinct partitioning of lineages among oceanic biomes[END_REF][START_REF] Scanlan | Ecological genomics of marine picocyanobacteria[END_REF][START_REF] Pittera | Connecting thermal physiology and latitudinal niche partitioning in marine Synechococcus[END_REF][START_REF] Sohm | Cooccurring Synechococcus ecotypes occupy four major oceanic regimes defined by temperature, macronutrients and iron[END_REF][START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF]. Finally, members of clade CRD1 were recently found to be dominant in large Fe-depleted areas of the world Ocean [START_REF] Sohm | Cooccurring Synechococcus ecotypes occupy four major oceanic regimes defined by temperature, macronutrients and iron[END_REF][START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF]. Even though clades globally occupy distinct niches, it was also shown that distinct ecotypes within Prochlorococcus and Synechococcus clades can display specific distribution patterns [START_REF] Mazard | Multi-locus sequence analysis, taxonomic resolution and biogeography of marine Synechococcus[END_REF][START_REF] Kashtan | Single-cell genomics reveals hundreds of coexisting subpopulations in wild Prochlorococcus[END_REF]Larkin et al., 2016a), with for instance distinct genetic groups within clades II and CRD1 colonizing different thermal niches [START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF].

Despite good knowledge of both their genetic diversity and environmental preferences, little is known about how environmental factors influence genome diversity and shape the community structure of marine picocyanobacteria, especially for Synechococcus. However, the development of high throughput sequencing techniques now allows such questions to be addressed. In particular, comparative genomics approaches applied to bacteria have revealed the high variability of microbial gene content, even for closely related strains sometimes displaying identical 16S rRNA sequences Konstantinidis and Tiedje, 2005b). They notably led to the definition of i) the core genome, the conserved part of the genome that encompasses genes shared by all strains, and ii) the flexible genome, the content of which is much more variable and dependent on the local biotic and abiotic environment [START_REF] Lan | Intraspecies variation in bacterial genomes: The need for a species genome concept[END_REF][START_REF] Cordero | Explaining microbial genomic diversity in light of evolutionary ecology[END_REF]. In cyanobacteria, previous studies based on multiple genome comparisons have shown that these organisms still present a so-called 'open pan-genome' [START_REF] Tettelin | Genome analysis of multiple pathogenic isolates of Streptococcus agalactiae: implications for the microbial "pan-genome[END_REF][START_REF] Baumdicker | The infinitely many genes model for the distributed genome of bacteria[END_REF][START_REF] Simm | The composition of the global and feature specific cyanobacterial core-genomes[END_REF]. Indeed, each newly sequenced genome brings novel genes without diversity saturation, and this holds true for Prochlorococcus and Synechococcus, for which only 17 [START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF]Biller et al., 2014a) and 14 genomes, respectively [START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF][START_REF] Baumdicker | The infinitely many genes model for the distributed genome of bacteria[END_REF] have so far been compared. These studies Genome diversification in marine picocyanobacteria thus highlight that the genomic diversity of natural populations is still mostly under-sampled, which strongly limits the interpretation of comparative genomic analyses. Here, we use a dataset of 81 nonredundant genomes of marine or halotolerant picocyanobacteria, of which 34 are newly sequenced complete Synechococcus genomes, to further assess the genomic diversity within these genera and how occupancy of new realized niches has impacted the evolution of these genomes. Analysis of this unprecedented genome dataset with original bioinformatic tools allowed us to estimate the relative contribution of gene gains/losses and sequence divergence on the diversification of marine picocyanobacteria and to highlight key processes involved in their adaptation to various environmental niches.

Results

Picocyanobacteria exhibit a wide intra-clade genomic diversity

In order to expand the coverage of Synechococcus in available marine picocyanobacterial genomes, 34 new strains were sequenced from cultured isolates, resulting in a quasi-doubling of the current number of complete or near-complete genomes publicly available for this genus. Strains were selected to cover the extent of the phylogenetic and pigment diversity of Synechococcus, as well as maximize their geographic origin and trophic regimes of their isolation site (Fig. 1, Supplementary Table S1). It should be noted though, that no cultured isolates are available yet for the EnvA and EnvB clades [START_REF] Mazard | Multi-locus sequence analysis, taxonomic resolution and biogeography of marine Synechococcus[END_REF][START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF]. The use of Wisescaffolder [START_REF] Farrant | WiseScaffolder: an algorithm for the semi-automatic scaffolding of Next Generation Sequencing data[END_REF] allowed us to close 28 out of the 31 genomes sequenced by the Genoscope and the Center for Genomic Research, with only one gap remaining in strains RS9915 and BOUM118 (both in the giant gene swmB [START_REF] Brahamsha | An abundant cell-surface polypeptide is required for swimming by the nonflagellated marine cyanobacterium Synechococcus[END_REF][START_REF] Mccarren | SwmB, a 1.12-megadalton protein that is required for nonflagellar swimming motility in Synechococcus[END_REF] and three gaps in strain BIOS-E4-1 (two in genes encoding a PQQ enzyme repeat family protein and one in an LVIVD repeat family protein). This high-quality genome dataset constitutes a key asset for comparative genomics analyses. Consistent with the genome streamlining that occurred in most Prochlorococcus lineages [START_REF] Dufresne | Accelerated evolution associated with genome reduction in a free-living prokaryote[END_REF][START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF][START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF] S1). Of note, members of the cold-adapted Synechococcus clades I and IV exhibited the lowest GC% values of all Synechococcus/Cyanobium strains (53.8 ± 0.73%) and this difference is even more marked using GC3%, i.e. the GC content at the third codon position (56.7 ± 1.25%; Fig. 2; p < 10 -8 Wilcoxon test for clades I and IV vs. all other Synechococcus/Cyanobium). By contrast, the warm-adapted clades II and III displayed significantly higher values (70.2 ± 1.5%; p < 10 -5 Wilcoxon test clades II and III vs. clades I and IV), while the highest GC3% was found for members of the brackish strains of Synechococcus clade VIII and SC 5.2 (81.1 ± 4.6%; p<10 -5 Wilcoxon test clade VIII and SC 5.2 vs. all other Synechococcus). Thus, although the strongest GC3% variation was associated to genome reduction in Prochlorococcus, some of the GC3% variations might be related to the ecological niches occupied by these organisms and notably to thermal and variable salinity niches [START_REF] Fuller | Clade-specific 16S ribosomal DNA oligonucleotides reveal the predominance of a single marine Synechococcus clade throughout a stratified water column in the Red Sea[END_REF].

Although they all belong to a monophyletic, long diverged branch within the cyanobacteria radiation [START_REF] Shih | Improving the coverage of the cyanobacterial phylum using diversity-driven genome sequencing[END_REF][START_REF] Sánchez-Baracaldo | Origin of marine planktonic cyanobacteria[END_REF], picocyanobacterial genomes show a tremendous diversity of both nucleotide sequences and gene content. Average nucleotide identity (ANI) and average amino acid identity (AAI) between pairs of picocyanobacterial genomes indeed ranged from

Genome diversification in marine picocyanobacteria

54.1 to 99.9% and 53.16 to 98.9%, respectively and intra-clade ANI and AAI were on average 91.8% and 91.04% (Fig. 3A, Supplementary Fig. S1). Thus, members of a given clade and even in most cases a given sub-clade, displayed ANI greater than 95%, classically used to define microbial species (Konstantinidis and Tiedje, 2005a;[START_REF] Goris | DNA-DNA hybridization values and their relationship to whole-genome sequence similarities[END_REF]. Interestingly, Synechococcus clades I and IV showed a particularly low ANI with other Synechococcus strains, while their ANIs with Prochlorococcus genomes were higher than for other Synechococcus-Prochlorococcus pairs. Since we did not observe this specificity with AAI, it is likely due to the low GC3% of Synechococcus clades I and IV (Fig. 2).

A plot of the relationship between 16S rRNA identity and AAI for the different pairs of genomes (Fig. 3B) additionally showed two major discontinuities. The first one at 80% AAI discriminated pairs of strains of the same clade from pairs of strains from different clades. Notable exceptions concerned the closely related and globally scarce Synechococcus clades V and VI as well as clades XX and UC-A, which fall within the intra-clade divergence level in terms of 16S rRNA identity and AAI, and Prochlorococcus clade LLII-III, which showed a divergence level similar to Synechococcus intra-SC divergence, suggesting that the gathering of these two clades into a single clade [START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF]Biller et al., 2014b) should be reconsidered, as suggested by Yan et al [START_REF] Yan | Genome rearrangement shapes Prochlorococcus ecological adaptation[END_REF]. The second discontinuity set apart Synechococcus strains of the same SC from strains of different SC (<98% 16S rRNA identity and <65% AAI), reflecting a very ancient genomic diversification between the three SC (see below). Because of this clear discontinuity, we propose to split the marine Synechococcus group into three distinct taxonomic groups: Ca. Marinosynechococcus (SC 5.1), Cyanobium (SC 5.2) and Ca. Juxtasynechococcus (SC 5.3). Prochlorococcus strains from different LL clades also fell below the 65% AAI discontinuity, highlighting the large divergence within this group. It is noteworthy that strains within SC 5.2 displayed a particularly low 16S rRNA identity compared to strains within SC 5.1, likely due to the low number of sequenced genomes relative to the wide diversity of this lineage, while in contrast the only two Synechococcus SC 5.3 genomes of our dataset were very closely related.

In order to manually refine the annotation of these genomes and ease comparative genomic analyses in terms of gene content, all genomes were included in the Cyanorak v2.1 information system (www.sb-roscoff.fr/cyanorak; Garczarek et al., submitted), in which predicted genes were grouped into clusters of likely orthologous genes (CLOGs) by all-against-all sequence similarity. This clustering allowed us to determine the core genome, i.e. CLOGs present in all strains, and the pan-genome, i.e. all CLOGs present in at least one strain, at various phylogenetic depths [START_REF] Tettelin | Genome analysis of multiple pathogenic isolates of Streptococcus agalactiae: implications for the microbial "pan-genome[END_REF]. When considering the whole dataset, the number of core CLOGs as a function of the number of genomes showed an asymptotic decline, tending toward a core set of 911 genes (Fig. 4B). In contrast, the pangenome of marine picocyanobacteria, containing 27,376 CLOGs, was still far from saturation, revealing that even with 81 genomes, every newly sequenced picocyanobacterial genome still brought about 192 new genes. This result held true when considering Prochlorococcus (7,537 CLOGs) and Synechococcus (20,986 CLOGs) independently, indicating that we still missed an essential part of the genetic diversity within both genera that is yet to be sequenced from the field. A major asset brought by the 34 newly sequenced Synechococcus genomes is the availability of several genomes per clade, which allowed us to estimate the relative sizes of the core set of CLOGs at different taxonomic levels (i.e. genus, SC and clades), the accessory genome, i.e. CLOGs shared by at least two strains but not core, and unique genes, i.e. CLOGs present in a single strain (Fig. 4A, Supplementary Table S2). While the proportion of accessory genes was pretty constant between genomes, constituting on average 13 ± 2.4% and 20.7 ± 6.3% of the Prochlorococcus and Synechococcus genomes respectively, unique genes constituted the most variable part of the genomes, ranging from 0.6% to 21.9% and 1.5% to 31.2% of the Prochlorococcus and Synechococcus genomes respectively, and were directly related to genome size. The newly sequenced strain BIOS-E4-1 (clade CRD1) contained by far the largest gene number of the genome dataset (4,426 genes), with a large proportion of unique genes (31.2%). Noteworthy, a significant proportion of CLOGs was present in all strains of a given clade (e.g. 335 genes for Synechococcus clade III, or 143 genes for Prochlorococcus HLI) and could thus potentially be involved in the adaptation of these taxa to specific environmental conditions. However, it should be noted that only a sub-set of these CLOGs were truly specific to each clade (e.g. 32 and 11 genes present in clades III and HLII, respectively; Supplementary Table S3) or ecologically significant taxonomic units (ESTU sensu [START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF]; see below and Supplementary Table S4), that is absent from all other clades or ESTUs.

Dynamics of the evolution of gene content in marine picocyanobacteria

To better understand the evolutionary processes that led to the diversification of gene content within marine picocyanobacterial genomes, we estimated by Maximum Likelihood the number of gene gain and loss events on each branch of a reference phylogenetic tree built from a concatenation of 821 single core proteins (Fig 5). As previously observed [START_REF] Dufresne | Accelerated evolution associated with genome reduction in a free-living prokaryote[END_REF][START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF], the gain and loss values obtained for Prochlorococcus were consistent with the scenario of a major genome streamlining process that occurred during the evolution of this genus, since an excess of gene loss was observed at the base of this radiation (Fig. 5). Globally, the number of genes gained and lost on each branch of the picocyanobacterial tree was quite variable. While on internal branches the number of gains and losses remained limited and balanced for both Prochlorococcus and Synechococcus SC 5.1 (gains ≤378, losses ≤479; not taking into account the genome streamlining at the base of the Prochlorococcus radiation), a higher number of events were generally observed on terminal branches as well as an excess of gains compared to losses, with up to 1,662 gained genes on the branch leading to Synechococcus BIOS-E4-1 (SC 5.1) and 831 on the one leading to Prochlorococcus MIT0701, for 105 and 108 lost genes, respectively.

By using calibration time points from a previous study [START_REF] Sánchez-Baracaldo | Origin of marine planktonic cyanobacteria[END_REF], we estimated that this corresponds to about 0.71 and 4.62 genes gained (1.67 and 1.80 genes lost) per million years (My) on internal and terminal branches of Synechococcus SC 5.1, respectively, while internal and terminal branches of Prochlorococcus HL gained 1.45 and 4.5 genes (0.87 and 3.72 lost; Table 1). The higher values observed for the terminal branches are related to the high number of strain-specific genes and reflect the fact that most of the variability in gene content occurs at the 'leaves' of the tree. If we assume the rate of gene gain to be constant over time, this suggests that most of the genes gained on internal branches have been secondarily lost and are therefore not represented in our genomic dataset.

As genomic islands have been shown to play a key role as repositories of laterally transferred genes potentially involved in niche adaptation in marine picocyanobacteria [START_REF] Coleman | Genomic islands and the ecology and evolution of Prochlorococcus[END_REF][START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF][START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF][START_REF] Delmont | Linking pangenomes and metagenomes : the Prochlorococcus metapangenome[END_REF][START_REF] Yan | Genome rearrangement shapes Prochlorococcus ecological adaptation[END_REF], we explored the contents of these islands in all analyzed genomes. Most genomes were too distant to compare genomic islands between strains by whole genome alignment as performed by [START_REF] Coleman | Genomic islands and the ecology and evolution of Prochlorococcus[END_REF] on Prochlorococcus, so here genomic islands were defined in each strain as regions of the genome enriched in gained genes using a similar approach as Kashtan et al. [START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF] but using a threshold to define the limits of the islands in each strain (see Methods; Supplementary Table S5). The number of gained genes located in genomic islands and shared by pairs of strains showed that closely related strains share many more island genes than distantly related ones and that only a few exchanges of genes occur between distantly related clades (Supplementary Fig. S2). These observations are

Genome diversification in marine picocyanobacteria

particularly striking for Prochlorococcus HL streamlined genomes that share only a low proportion of island genes with Synechococcus. A notable exception is Synechococcus clade VIII, which shares more island genes with strains of SC 5.2 than with most SC 5.1 strains, an expected pattern since these groups co-occur in coastal or estuarine waters of variable salinity [START_REF] Fuller | Clade-specific 16S ribosomal DNA oligonucleotides reveal the predominance of a single marine Synechococcus clade throughout a stratified water column in the Red Sea[END_REF][START_REF] Chen | Diverse and unique picocyanobacteria in Chesapeake Bay , revealed by 16S-23S rRNA internal transcribed spacer sequences[END_REF][START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF]. To further explore how strains share genomic islands, we used an innovative network method based on the partial similarity of gene contents between islands shared by pairs of strains. It allowed us to retrieve islands previously identified either by direct pairwise comparison of Prochlorococcus HLI MED4 and HLII MIT9312 strains [START_REF] Coleman | Genomic islands and the ecology and evolution of Prochlorococcus[END_REF] or by analyzing the deviation in tetranucleotide frequency in individual Prochlorococcus and Synechococcus genomes [START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF], demonstrating the validity of our automated approach (Supplementary Fig. S3 and Supplementary Tables S6 andS7). Interestingly, most islands identified by these authors in Prochlorococcus HL strains appeared to be shared by all HL strains, forming dense red, knot-shaped modules in the network (e.g., Pro_GI033 = MED4 ISL1; Pro_GI048 = MED ISL2; Pro_GI028 = MIT9312 ISL5; Pro_GI000 = MED SVR2; Pro_GI015 = MED4 SVR4; Pro_GI041 = MED4 ISL1.1; Pro_GI023 = MED4 ISL2.2; Fig. 6 and Supplementary Table S6). These red knots correspond to genomic regions prone to gene integration that have likely been acquired by the common ancestor of all HL strains, then vertically transferred to all descendants, much like the phycobilisome region that is shared by all Synechococcus strains [START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF]. In contrast, ISL4 island, initially identified in MED4 by [START_REF] Coleman | Genomic islands and the ecology and evolution of Prochlorococcus[END_REF] and later confirmed both by Dufresne and coworkers [START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF] and our automated island detection approach (Pro_GI004; Supplementary Fig. S3), does not form a red knot but only a fuzzy network of interconnected islands, each shared by only 2 to 4 strains (Fig. 6). So this island, whose gene content is highly variable, has seemingly been more recently acquired by a subset of the HL population. Our approach also unveiled previously undescribed islands specifically shared by sets of Prochlorococcus LL strains, including Pro_GI027, 039, 044 and 049 specific to LLI strains (several being enriched in hli genes, known to be amplified in LLI compared to other LL strains; (Partensky and Garczarek, 2010b), Pro_GI010, 018 and 025 specific to LLII/III strains, and Pro_GI002 as well as 13 other modules specific to LLIV strains, including several containing genes encoding lanthipeptides [START_REF] Tang | Structural characterization of four prochlorosins: A novel class of lantipeptides produced by planktonic marine cyanobacteria[END_REF]Fig. 6 and Supplementary Table S6).

In Synechococcus, the network included relatively few dense red knots compared to Prochlorococcus (Fig. 7). Among the most notable ones are three clade III-specific islands: the first one (Syn_GI013) gathers a gene cluster (cynA-B-D) involved in cyanate transport [START_REF] Kamennaya | Characterization of cyanate metabolism in marine Synechococcus and Prochlorococcus spp[END_REF]; Supplementary Table S7); the second one (Syn_GI087) includes a specific beta-glycosyltransferase and swmA, a protein involved in a special type of motility observed only in members of this clade [START_REF] Mccarren | Inactivation of swmA results in the loss of an outer Cell layer in a swimming Synechococcus strain[END_REF]; the third one (Syn_GI102) notably contains swmB, encoding a giant protein also involved in this motility process [START_REF] Mccarren | SwmB, a 1.12-megadalton protein that is required for nonflagellar swimming motility in Synechococcus[END_REF]. Another interesting example is Syn_GI100, which notably encompasses a 3-gene cluster composed of one nfeD homolog and two flotillin-like genes that both have similarity to the floT gene involved in the production of lipid rafts, whose deletion in Bacillus subtilis was found to strongly affect cell shape and motility [START_REF] Dempwolff | The deletion of bacterial dynamin and flotillin genes results in pleiotrophic effects on cell division, cell growth and in cell shape maintenance[END_REF]. Interestingly, this gene cluster was found in the only two clade III strains (A15-24 and A15-28) that lack swmA and swmB as well as in several distantly related strains. Conversely, no swmA-Bcontaining strain was found to possess the nfeD-floT1-floT2 gene cluster. The network approach also detected quite a few knots containing both red and blue edges. The latter color indicates that strains sharing these islands are distantly related to one another. Thus, knots that are mixing red and blue edges potentially emphasize relatively recent horizontal gene transfers between clades or longer phylogenetic distances. This includes i) Syn_GI022, a module found in many SC 5.1 strains with the notable exception of clade II strains, which encompasses a large gene cluster involved in glycine betaine synthesis (gbmt1-2) and transport (proV-W-X), located in some strains next to another gene cluster

Genome diversification in marine picocyanobacteria

involved in the biosynthesis of glucosylglycerate (gpgS-gmgG-gpgG; [START_REF] Scanlan | Ecological genomics of marine picocyanobacteria[END_REF] and ii) Syn_GI122, a module comprising strains from almost all lineages that encompasses genes encoding uncharacterized cell surface proteins, secreted CHAT domain-containing proteins and/or genes involved in the biosynthesis of cyclic AMP (cAMP), including adenylate cyclases located in the vicinity of cyclic nucleotide-binding proteins, such as the cAMP receptor protein (CRP) or a cAMPregulated small-conductance mechanosensitive ion channel. Altogether, this network approach nicely complements the detection of genomic islands in single genomes by providing insights about the evolutionary history of these genomic islands.

Relative contributions of variability at the sequence and gene content levels in the evolution of picocyanobacteria

The fairly low rate of gene acquisition evidenced in this study raises the question of the relative weight of gene content variations vs. substitutions in the nucleotide sequence in the long-term diversification and adaptation processes of these organisms. Figure 8 compares a phylogenetic tree built with a concatenation of 821 picocyanobacterial core protein sequences to a dendrogram based on the phyletic pattern (i.e. the pattern of presence/absence of each CLOG in each strain). Topologies of the two trees were globally similar, which reveals that fixation of genes and fixation of mutations occurred concomitantly during the evolutionary history of marine picocyanobacteria. Yet, Synechococcus clade VIII and SC 5.2 were found to be closely related in the dendrogram based on the phyletic pattern. Indeed, as previously reported in a study using 11 Synechococcus genomes [START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF], these taxa share a fair number of genes, potentially related to their co-occurrence in brackish environments. Interestingly, the closely related clades V and VI cluster together with these two taxa, indicating that they may also share with clade VIII and SC 5.2 some mechanisms of adaptation to low salinity niches (see below). Although the presence of SC 5.3 has been recently documented in freshwater environments [START_REF] Cabello-Yeves | Novel Synechococcus genomes reconstructed from freshwater reservoirs[END_REF], the presence of the two marine sequenced strains (RCC307 and MINOS11) at the base of this halotolerant group might instead be due to attraction by SC 5.2.

Among the Synechococcus SC 5.1 and Prochlorococcus radiations, we identified a few incongruences between the two trees within Synechococcus clades I, II, III and VI and Prochlorococcus HLII (Fig. 8) that are likely due to the relatively low number of specific genes within these clades. It is also worth noting that some clades were closer in terms of gene content than expected from the core phylogeny, in particular Synechococcus clades WPC1, XX and UC-A grouping with clade III in the tree based on the phyletic pattern. Finally, some clades lost their monophyly in the tree based on phyletic pattern, such as Synechococcus clades V and VI that were mixed together or Prochlorococcus HLI that was found to be mixed with HLII. This example is particularly interesting, since despite their clearly distinct phylogenetic clustering based on protein sequences and well-known ecological and physiological differences [START_REF] Johnson | Niche partitioning among Prochlorococcus ecotypes along ocean-scale environmental gradients[END_REF]Martiny et al., 2009b), these two clades have a quite similar gene content, with only a few genes (29) present in all HLII strains but not in all HL strains (Fig. 4A).

Similarly, Prochlorococcus clade LLI, which was previously shown to occupy an intermediate niche between HL and strict LL members (LLII-IV) and to share genes with both ecotypes [START_REF] Johnson | Niche partitioning among Prochlorococcus ecotypes along ocean-scale environmental gradients[END_REF]Partensky and Garczarek, 2010a), actually appeared to share more genes with the LLII-III clade (1,382 genes) than with HL (1,290 genes). Altogether, these two examples show that within Prochlorococcus, although HL and LL have different gene contents, differentiation within HL and to

Genome diversification in marine picocyanobacteria

a lesser extent within LLI-III rather relies on substitution accumulations than on variation in gene content.

Another major difference between these trees concerned branch lengths. By computing for each node at the base of a clade (blue dots in Fig. 8) the average length from the node to its descending leaves (terminal length), and the length from the node to its parent node (internal length), we showed that the ratio of terminal to internal branch lengths was significantly higher (Mann-Whitney paired test, p-value <0.0015) in the phyletic pattern tree than in the core tree (Supplementary Fig. S4). This suggests that there were more amino acid substitutions before the divergence of clades than after, whereas there was more gene content variation between strains of a clade than between clades. In other words, this comparison revealed that most of the changes that were fixed in the long term by evolution are substitutions and not changes in gene content.

In order to quantify more precisely this difference, we compared the estimated number of gene gains and losses per My (Supplementary Fig. S5) to the number of amino acid substitutions in core proteins per My (Supplementary Fig. S6) and results of these comparisons are shown in Table 1. It is important to note that the rates of gene gain/loss and amino acid substitutions calculated this way should only be considered as lower bound estimates for several reasons. First, since we only have access to the presentday genomes and not to ancestral ones, measurements of the rate of genes gained in fact refer to genes gained and successfully retained over time in at least one strain. Second, the amino acid substitution rates were measured on core proteins, whose genes likely undergo a strong purifying selection. This, together with the much longer generation time of picocyanobacteria compared to model bacteria and with their considerable population size (Partensky et al., 1999b;[START_REF] Dufresne | Accelerated evolution associated with genome reduction in a free-living prokaryote[END_REF][START_REF] Flombaum | Present and future global distributions of the marine Cyanobacteria Prochlorococcus and Synechococcus[END_REF], could explain why estimated rates were lower than for other bacterial lineages [START_REF] Lawrence | Molecular archaeology of the Escherichia coli genome[END_REF][START_REF] Mcdonald | Lateral gene transfer dynamics in the ancient bacterial genus Streptomyces[END_REF]. With this caveat in mind, in Prochlorococcus HL, 356x more amino acid substitutions than gene gains were estimated for internal branches per My, and 69.6x for terminal branches, primarily due to a higher rate of gene gain in the latter branches. In Synechococcus SC 5.1, a ratio of 164 and 20 was obtained for internal and terminal branches, respectively, the difference between the two genera likely being due to the higher rate of protein sequence evolution observed in Prochlorococcus [START_REF] Dufresne | Accelerated evolution associated with genome reduction in a free-living prokaryote[END_REF].

We also compared at each node the fixation rate of amino acid substitutions in core proteins (i.e. amino acids in the alignment that are identical in all descending strains and different in all other strains) to the fixation rate of genes (i.e. present in all descending strains and in no other strain). 201x more amino acid variants than genes were fixed per My in Prochlorococcus HL (and 116x more for Synechococcus SC 5.1). This corresponds to a fixation rate of 78 and 18 amino acid changes in core proteins per My for Prochlorococcus HL and Synechococcus SC 5.1, respectively, while one gene is fixed once every 2.6 My for Prochlorococcus HL and once every 6.3 My for Synechococcus SC 5.1. While these numbers show that substitutions played a major role in genomic diversification, the question remains as what part of this diversification is related to an adaptive process.

Role of gene content in the adaptation of Synechococcus to specific niches

In contrast to Prochlorococcus [START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF]Partensky and Garczarek, 2010a;Biller et al., 2014b;[START_REF] Delmont | Linking pangenomes and metagenomes : the Prochlorococcus metapangenome[END_REF][START_REF] Yan | Genome rearrangement shapes Prochlorococcus ecological adaptation[END_REF], few genomic diversity studies have been conducted so far in Synechococcus. In order to reveal whether the presence or absence of genes might be related to Synechococcus adaptation to specific niches, we defined sets of clades co-occurring in the field and

Genome diversification in marine picocyanobacteria

occupying similar niches, based on assemblages of ESTUs as defined in [START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF]. We then searched for genes occurring in strains within a given set and absent from other picocyanobacterial strains using a relaxed, niche-related definition of specificity (Supplementary Table S4). These analyses revealed only 18 CLOGs specific to members of both cold thermotypes, clades I and IV, among which 6 had a putative function, though with seemingly no direct relationship with adaptation to low temperature. However, the set of 19 CLOGs specific to clade I includes a particular isoform of the chaperone protein DnaK (DnaK4, CK_00056929; Supplementary Table S3) in addition to the three gene copies present in most Synechococcus SC 5.1 strains. This additional copy might be involved in protein folding in cold conditions [START_REF] Genevaux | The Hsp70 chaperone machines of Escherichia coli: A paradigm for the repartition of chaperone functions[END_REF].

Members of clades III, WPC1 and SC 5.3, co-occurring in warm, P-depleted oligotrophic waters, were found to share a much higher number of genes (85; Supplementary Table S4), among which 2 were previously reported to be related to phosphate availability: a yet uncharacterized gene (CK_00002088) found to be downregulated in early phosphate stress [START_REF] Tetu | Microarray analysis of phosphate regulation in the marine cyanobacterium Synechococcus sp[END_REF] and a chromate transporter (ChrA), which was recently suggested to be involved in phosphate acquisition in Prochlorococcus, based on its enrichment in P-poor oligotrophic areas [START_REF] Kent | Global biogeography of Prochlorococcus genome diversity in the surface ocean[END_REF]. Clades III and WPC1 also share a cluster of 12 consecutive genes potentially involved in capsular polysaccharide synthesis and export (including genes similar to kps genes in Escherichia coli K1, responsible for the formation of a polysialic acid extracellular capsule; see Kps cluster in Supplementary Table S4) and another cluster of 7 genes that might be involved in the use of organic nitrogen sources since it encompasses a putative nitrilase (CK_00002256). Additionally, 32 genes were found to be specific to the 8 clade III strains, including the above-mentioned cyanate transporter genes (cynABD; [START_REF] Kamennaya | Characterization of cyanate metabolism in marine Synechococcus and Prochlorococcus spp[END_REF] as well as a phosphate starvation-induced protein (PsiP1; [START_REF] Scanlan | Molecular ecology of the marine cyanobacterial genera Prochlorococcus and Synechococcus[END_REF] and a specific alkaline phosphatase (CK_00052500) that potentially hydrolyses extracellular organic phosphates (Supplementary Table S3). Similarly, the two members of SC 5.3 also share a large number of strictly specific genes (215), including a regulator of phosphate uptake (PhoU; CK_00005756; diCenzo et al., 2017) as well as two putative phosphatases (CK_00005504, CK_00005619) and a putative pyrophosphatase (CK_00005811), in addition to the 4 potential pyrophosphatases present in most picocyanobacterial genomes (CK_00000642, CK_00000654, CK_00000805 and CK_00008108; Supplementary Table S3). Altogether, these results suggest that the occurrence of these genes might contribute to the success of clade III, WPC1 and SC 5.3 cells in oligotrophic, P-depleted environments such as the Mediterranean Sea in summer [START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF], and indicates that members of these three taxa have adopted partially different strategies to cope with P depletion. To further explore the adaptive strategies of these clades to cope with low inorganic P concentrations, we compiled a Table displaying the number of copies of each CLOG related to P transport and metabolism in all Synechococcus strains (Supplementary Table S8). All clade III strains share at least three copies encoding the PstS transporter and one copy of sphX, in addition to the ChrA transporter mentioned above. The number of transporters is also high in clades VIII, WPC1 and members of SC5.2, while it is systematically lower in clades I, II, IV, VII and CRD1. Interestingly, clade I, II and IV strains virtually all sphX, with only one clade II strain (A15-44) possessing this gene. All members of clades I and IV also lack the genes phoB and phoR coding for the two-component system involved in P sensing and the regulation of P metabolism, as previously observed on fewer strains in [START_REF] Scanlan | Ecological genomics of marine picocyanobacteria[END_REF]. While all clades have the genetic potential for phosphonate utilization, only some clade II strains and a single strain from clade III (A15-28) possess the genes for phosphite assimilation. This trait is however not conserved at the clade level. Finally, this detailed analysis revealed the particularly high number of shared phosphatase genes in clades III (8 to 12 genes) and WPC1 (8 and 11 genes, median = 9), in contrast to the lower number observed in clades I, II, IV and VII (3 to 6, median = 5). This suggests an adaptive strategy to diversify sources of organic phosphate available to members of these clades, likely as an adaptation to environments depleted in inorganic P. Clade VIII seems to have

Genome diversification in marine picocyanobacteria

specialized in a specific organic source with 3 or 4 copies of the same phosphatase while clades V, VI, CRD1 and SC5.2 have more variable numbers of phosphatases, reflecting strain-level variation rather than clade-level strategies.

Genes potentially involved in niche adaptation were also found in all three strains of the CRD1 clade, known to dominate in iron-depleted oceanic regions, which share a quite high number of specific CLOGs (81, Supplementary Tables S3 andS4), though most of them have no known function. Among the characterized ones were a second copy of the flavodoxin IsiB, a Cu-containing protein known to replace ferredoxin in iron-depleted conditions (Erdner and Anderson, 1999), the ferrous iron transport protein FeoA, an iron-sulfur cluster biosynthesis family protein (CK_00008433) as well as 3 specific high light-induced proteins (HLIPs) that might provide protection from oxidative stress to photosystems [START_REF] He | The high light-inducible polypeptides in Synechocystis PCC6803[END_REF].

Finally, in agreement with their clustering in the dendrogram based on phyletic pattern (Fig. 8), clades VIII and SC 5.2 share 28 genes including a few strictly specific genes (Supplementary Table S4), such as a fatty acid hydroxylase (CK_00002851) involved in lipid biosynthesis, and one or two copies of a P-type ATPase (CK_00045881), a family of ATP-driven pumps known to transport a variety of different ions and phospholipids across membranes [START_REF] Axelsen | Evolution of substrate specificities in the P-type ATPase superfamily[END_REF]. It is also noteworthy that SC 5.2 and clade VIII share a fair number of genes potentially involved in the adaptation to low salinity with members of clades V, VI and sometimes VII, whose ecological niches are still poorly known [START_REF] Zwirglmaier | Global phylogeography of marine Synechococcus and Prochlorococcus reveals a distinct partitioning of lineages among oceanic biomes[END_REF][START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF]Xia et al., 2017a) and possibly encompass environments with variable salinity (Supplementary Table S4). This includes a specific small-conductance mechanosensitive ion channel (MscS family) that might be involved in the response to osmotic stress (CK_00056919; [START_REF] Haswell | Mechanosensitive channels: What can they do and how do they do it?[END_REF] and a bacterial regulatory protein of the ArsR family that besides regulating the efflux of arsenic and arsenite was suggested to participate in salt tolerance in Staphylococcus aureus through a Na + efflux activity [START_REF] Scybert | NaCl-sensitive mutant of Staphylococcus aureus has a Tn917-lacZ insertion in its ars operon[END_REF]. In addition, members of clade VIII share 22 specific genes, including a second potential mechanosensitive ion channel (MscS; CK_00056915), while members of SC 5.2 share 31 specific genes, including another mscS gene copy (CK_00003081) as well as genes encoding a putative chloride channel (CK_00042275) and a NAD-dependent malic enzyme, a protein known to be enhanced under salt stress in plants [START_REF] Liu | Expression of an NADPmalic enzyme gene in rice (Oryza sativa L) is induced by environmental stresses: overexpression of the gene in Arabidopsis confers salt and osmotic stress tolerance[END_REF]; Supplementary Table S3). Despite these few examples, it seems that the number of genes potentially related to the ecological niche occupied by each clade or assemblage of clades is fairly limited and varies depending on the considered niche, with for instance few genes related to thermal niche adaptation. Most of the diversity in gene content therefore relies on differences between individual strains rather than between phylogenetic groups or ESTUs, a large proportion of the sparsely distributed genes having yet unknown functions, some potentially being involved in niche adaptation.

Role of substitutions in adaptation

Given our observation that a high number of amino acid substitutions have been fixed in the long term, we also searched for those potentially involved in niche adaptation. We identified "specific variants" as positions in core protein alignments for which a particular amino acid is found in all strains of a given clade, ESTU or set of ESTUs and a different amino acid is found in other strains. In order to reduce the noise due to the accumulation of clade-specific substitutions and to better identify the niche

Genome diversification in marine picocyanobacteria

adaptation signal, we focused on variants shared by clades I and IV, which do not form a monophyletic group (Fig. 8, left) but usually co-occur in cold, temperate waters [START_REF] Zwirglmaier | Basin-scale distribution patterns of picocyanobacterial lineages in the Atlantic Ocean[END_REF][START_REF] Zwirglmaier | Global phylogeography of marine Synechococcus and Prochlorococcus reveals a distinct partitioning of lineages among oceanic biomes[END_REF]Martiny et al., 2009b;[START_REF] Sohm | Cooccurring Synechococcus ecotypes occupy four major oceanic regimes defined by temperature, macronutrients and iron[END_REF][START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF][START_REF] Kent | Parallel phylogeography of Prochlorococcus and Synechococcus[END_REF]. We identified 180 proteins mainly involved in i) energy metabolism, ii) biosynthesis of cofactors, prosthetic groups and carriers, such as pigments and vitamins, iii) protein synthesis and protein fate, and to a lesser extent iv) transport and DNA metabolism (Supplementary Table S9). The first category encompassed proteins responsible for carbon fixation (RuBisCO subunits RbcS and RbcL, carbonic anhydrase CsoSCA, carboxysome proteins CsoS1E and CsoS2, and Calvin cycle enzyme Fbp-Sbp), two photosystem II subunits (the extrinsic PsbU protein and the manganese cluster assembly protein, Psb27) and a number of proteins involved in electron transport for photosynthesis and/or respiration (CtcAI, CtcEI, NdhA and two ATP synthase subunits: AtpA and AtpD). Furthermore, this set includes six proteins potentially involved in the response to light or oxidative stress: two High Light Inducible Proteins (HLIPs; CK_00001609 and CK_00001414), two peroxiredoxins (PrxQ), a glutaredoxin (CK_00000445) and a flavoprotein involved in the Mehler reaction (Flv1). We also identified a few enzymes involved in sugar metabolism and in particular in the pentose phosphate pathway (Pgl, TalA and Zwf). As concerns the 'protein synthesis' and 'protein fate' categories, this includes six ribosomal proteins and nine amino acid biosynthesis proteins, several tRNA/rRNA modification enzymes and tRNA aminoacyltransferases as well as seven proteins responsible for folding and stabilization of polypeptides. Of particular interest are the proteins belonging to the 'biosynthesis of cofactors, prosthetic groups and carriers' category, including enzymes involved in chlorophyll (HemC, ChlN, ChlB), cobalamin (CobO, cobQ, CobU-CobP) and carotenoid biosynthesis. The latter includes CrtE and GpcE, two enzymes involved in the phytoene biosynthesis pathway and CrtP, CrtQ and CrtL-b, the three enzymes catalyzing all the steps required to transform phytoene into β-carotene. It is also interesting to note that the five proteins displaying the largest number of specific substitutions relative to protein length are a putative ABC multidrug efflux transporter (CK_00008042; 19 positions specific to clades I and IV out of 607 amino acids), lycopene β-cyclase, responsible for the last step of βcarotene synthesis (CrtL-b; 7/347), the bifunctional enzyme fructose-1,6biphosphatase/sedoheptulose-1,7-biphosphate phosphatase involved in both Calvin cycle and glycolysis (7/347), the photosystem II manganese cluster assembly protein Psb27 (3/160) and the ribosomal protein RpmB (1/78). Even though the number of substitutions is not directly correlated to the level of selection pressure, the high proportion of specific substitutions in these proteins suggests that they have been subjected to positive selection and therefore have potentially a role in adaptation to cold environments.

Discussion

The availability of 81 complete and closed picocyanobacterial genomes with extensive manually refined annotations, including 34 novel Synechococcus, constitutes a key asset for comparative genomics analyses. With regard to previous studies (see e.g. [START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF][START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF][START_REF] Scanlan | Ecological genomics of marine picocyanobacteria[END_REF], sequencing of several strains for most major Synechococcus clades revealed that the extent of genomic diversity is tremendous, at all taxonomic levels including within clades and most sub-clades. As previously observed for SAR11 [START_REF] Nayfach | An integrated metagenomics pipeline for strain profiling reveals novel patterns of bacterial transmission and biogeography[END_REF][START_REF] Tsementzi | SAR11 bacteria linked to ocean anoxia and nitrogen loss[END_REF], ANI and AAI were indeed most often well below the cut-off of 95% (Fig. 3), usually considered to be the limit between bacterial species (Konstantinidis andTiedje, 2005b, 2005a;[START_REF] Jain | High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries[END_REF]. Thus, based on this cut-off, most clades within cluster 5 sensu [START_REF] Herdman | Form-genus XIII. Synechococcus[END_REF]) would correspond to one or even several species, as suggested by one research group [START_REF] Thompson | Genomic taxonomy of the genus Prochlorococcus[END_REF][START_REF] Coutinho | Comparative genomics of Synechococcus and proposal of the new genus Parasynechococcus[END_REF]Coutinho et al., , 2016a)). However, the delineation of so many species in a radiation that mostly exhibits a continuum

Genome diversification in marine picocyanobacteria in terms of within clade sequence identity (ID% range: 84 to 100%; Fig. 3B) would create more confusion than clarification as it would result in most cases into single-strain species, which cannot be clearly differentiated based on their fundamental (see e.g. [START_REF] Moore | Single-cell genomics unveiled a cryptic cyanobacterial lineage with a worldwide distribution hidden by a dinoflagellate host[END_REF][START_REF] Pittera | Connecting thermal physiology and latitudinal niche partitioning in marine Synechococcus[END_REF] and/or environmental realized niches [START_REF] Huang | Novel lineages of Prochlorococcus and Synechococcus in the global oceans[END_REF][START_REF] Sohm | Cooccurring Synechococcus ecotypes occupy four major oceanic regimes defined by temperature, macronutrients and iron[END_REF][START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF][START_REF] Kent | Global biogeography of Prochlorococcus genome diversity in the surface ocean[END_REF]. With this caveat in mind, it is clear that besides the Prochlorococcus lineage, there are three extremely divergent monophyletic groups within the marine Synechococcus/Cyanobium radiation (Sánchez- [START_REF] Sánchez-Baracaldo | Insights into the evolution of picocyanobacteria and phycoerythrin penes (mpeBA and cpeBA)[END_REF], which furthermore can be clearly discriminated based on 16S similarity vs. AAI plots (Fig. 3B), with an AAI divergence below the 65% limit that has been proposed to discriminate distinct genera [START_REF] Konstantinidis | Prokaryotic taxonomy and phylogeny in the genomic era: advancements and challenges ahead[END_REF]. Based on these criteria, our proposition to split the marine Synechococcus group into three distinct taxonomic groups: Ca. Marinosynechococcus (SC 5.1), Cyanobium (SC 5.2) and Ca. Juxtasynechococcus (SC 5.3). This proposal notably solves the inconsistency to have a mix of strains named Cyanobium spp. and Synechococcus spp. within SC 5.2, which should clearly all be called Cyanobium spp. For the universal acceptance of the revised taxonomy of this group and cyanobacteria at large [START_REF] Komárek | A polyphasic approach for the taxonomy of cyanobacteria: principles and applications[END_REF], both temporary names proposed for SC 5.1 and 5.3 as well as the potential definition of species within each of these radiations await validation by a large panel of cyanobacterial community members. In any case, any creation of new species within this group should likely take into account previously defined monophyletic clades and subclades as these phylogenetic groups have been used in most previous laboratory and environmental studies, whatever the genetic marker used [START_REF] Palenik | Cyanobacterial diversity in marine ecosystems as seen by RNA polymerase (rpoC1) gene sequences[END_REF][START_REF] Penno | Diversity of Synechococcus and Prochlorococcus populations determined from DNA sequences of the N-regulatory gene ntcA[END_REF]Ahlgren and Rocap, 2012b;[START_REF] Huang | Novel lineages of Prochlorococcus and Synechococcus in the global oceans[END_REF][START_REF] Mazard | Multi-locus sequence analysis, taxonomic resolution and biogeography of marine Synechococcus[END_REF][START_REF] Scanlan | Ecology of Cyanobacteria II: Their Diversity in Space and Time[END_REF].

The particularly high degree of genomic divergence occurring within Cyanobacteria Cluster 5 needs to be taken into account when putting results from comparative genomics of marine picocyanobacteria in the context of other highly sequenced bacterial groups such as pathogens and commensals [START_REF] Harris | Evolution of MRSA during hospital transmission and intercontinental spread[END_REF][START_REF] Kennemann | Helicobacter pylori genome evolution during human infection[END_REF][START_REF] Mather | Distinguishable epidemics of multidrug-resistant Salmonella typhimurium DT104 in different hosts[END_REF]. While high divergence and associated low level of synteny somehow limit the application of classical population genetics approaches, such as calculation of recombination rates (McDonald and Currie, 2017), our dataset is in contrast well suited to study the long-term evolutionary processes that have shaped the genomes of these abundant and widespread organisms in relation to their ecological niche occupancy. Comparative genomic analyses on marine picocyanobacteria have so far mainly focused on comparing gene repertoires from strains isolated from distinct niches, with the idea that niche adaptation largely relies on differential gene content [START_REF] Rocap | Genome divergence in two Prochlorococcus ecotypes reflects oceanic niche differentiation[END_REF][START_REF] Palenik | Genome sequence of Synechococcus CC9311: Insights into adaptation to a coastal environment[END_REF][START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF][START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF]. Here, a comparison of several strains per clade led in most cases to the identification of relatively few specific genes of known function that may confer a trait necessary for niche adaptation, even using relaxed stringency criteria (e.g. by selecting genes present in >80 or 90% of strains within a clade/ESTU assemblage and in <20 or 10% of others; Supplementary Tables S3 andS4). This may be due to the existence of an extended within-taxa microdiversity (Martiny et al., 2009b;[START_REF] Kashtan | Single-cell genomics reveals hundreds of coexisting subpopulations in wild Prochlorococcus[END_REF][START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF]Larkin et al., 2016b), where the more genomes in a taxon, the lower the number of genes found in all strains of this taxon. This fairly low number of niche-specific genes might also suggest that gene gain/loss, and fixation of these events during evolution, is a less prominent process to explain niche adaptation of marine picocyanobacteria than previously thought. Although lateral gene transfer is often considered to "commonly" occur between cells, and was notably shown to be involved in adaptation to nitrogen-or phosphorus-poor conditions in Prochlorococcus, no previous study explicitly stated the evolutionary time scale at which these adaptations took place [START_REF] Martiny | Phosphate acquisition genes in Prochlorococcus ecotypes: Evidence for genome-wide adaptation[END_REF](Martiny et al., , 2009a;;[START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF][START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF][START_REF] Scanlan | Ecological genomics of marine picocyanobacteria[END_REF][START_REF] Berube | Physiology and evolution of nitrate acquisition in Prochlorococcus[END_REF][START_REF] Yan | Genome rearrangement shapes Prochlorococcus ecological adaptation[END_REF]. Here, although the higher estimated rate of gene gains on the terminal branches of the phylogenetic tree indicates that most detectable events occurred fairly recently with regard to the long evolutionary history of both genera (Fig. 5, Table 1), adding time calibration to the tree led to an estimation of only 4.5 and 5.6 genes gained per My on terminal branches in Prochlorococcus HL and

Genome diversification in marine picocyanobacteria

Synechococcus SC 5.1 strains, respectively. Thus, gene gains appear to be rather rare events. Even though these rates are approximate due to uncertainties in time calibration and probably underestimated, they are entirely in line with those estimated for Prochlorococcus HLII populations, thought to have diverged a few million years ago but only possessing a dozen unique genes [START_REF] Kashtan | Single-cell genomics reveals hundreds of coexisting subpopulations in wild Prochlorococcus[END_REF]. Furthermore, in accordance with previous studies on other bacterial groups [START_REF] Lerat | Evolutionary origins of genomic repertoires in bacteria[END_REF][START_REF] Ochman | Examining bacterial species under the specter of gene transfer and exchange[END_REF][START_REF] Nowell | The extent of genome flux and its role in the differentiation of bacterial lineages[END_REF][START_REF] Mcdonald | Lateral gene transfer dynamics in the ancient bacterial genus Streptomyces[END_REF], the fact that rates of gene gain/loss are estimated to be higher on terminal branches of the tree (Supplementary Fig. S4), together with the high number of unique genes in every sequenced strain (Fig. 4A), clearly suggests that most recently acquired genes will not be kept in the long term in both genera. Our calculation indeed gives an approximate value of 1.45 and 0.71 genes gained and subsequently kept per My in Prochlorococcus HL and Synechococcus SC 5.1, respectively (Table 1). This low fixation rate suggests that most of the recently gained genes have no or little beneficial effect on fitness in the long term and that we observe them in genomes because purging selection has not deleted them yet [START_REF] Hao | The fate of laterally transferred genes: life in the fast lane to adaptation or death[END_REF][START_REF] Abby | Comparative genomics and the evolution of prokaryotes[END_REF][START_REF] Rocha | Evolutionary patterns in prokaryotic genomes[END_REF]. Still, these recently gained genes could be involved in more transient adaptation processes at the evolutionary scale such as biotic interactions (e.g. resistance to viral attacks or grazing pressure).

Such a result also has important implications for interpreting the role of flexible genomes in the context of adaptation to distinct niches. Indeed, genes conferring adaptation to a specific niche are mixed in the genomes with genes with no or little beneficial effect and are thus difficult to identifyin particular when they have only a putative function. The relatively low gene fixation rate that we observed (Table 1) also implies that flexible genes that are fixed within a clade (i.e. clade-specific genes) were gained tens of millions of years ago, and thus might be more reflective of past selective forces than of recent adaptation to newly colonized niches. In this context, genes specifically shared by Synechococcus clade VIII and SC 5.2 suggest that adaptation to low salinity environments was a critical factor in their differentiation from other taxa and the most parsimonious evolutionary scenario would be a lateral transfer of these genes from a SC5.2-like strain to the common ancestor of clade VIII, which might date back to 51.6 My (confidence interval 0-141 My). Similarly, adaptation to phosphorus-depleted oligotrophic areas might have driven the differentiation of Synechococcus clade III, as revealed by the occurrence of P-and other nutrient-uptake genes specific to this clade. Interestingly, co-occurring ESTUs IIIA, WPC1A and SC 5.3A only share a few common genes potentially involved in the adaptation to this limitation. Instead, these ESTUs seem to have independently acquired different sets of genes to improve P-uptake and/or assimilation and potentially use different sources of organic phosphate (see Results and Supplementary Tables S4 andS8). It is notable that some clade II strains have also potentially adapted to inorganic P depletion by acquiring or conserving the ability to use phosphite. It is also noteworthy in this context that in Prochlorococcus, P metabolism is not claderelated but dependent on within-clade variability in the gene content of specific genomic islands [START_REF] Martiny | Phosphate acquisition genes in Prochlorococcus ecotypes: Evidence for genome-wide adaptation[END_REF](Martiny et al., , 2009a)), further highlighting the variety of evolutionary paths that led to adaptation to low-P environments in these different lineages.

As proposed recently for other bacterial model organisms [START_REF] Thrash | Singlecell enabled comparative genomics of a deep ocean SAR11 bathytype[END_REF][START_REF] Mcdonald | Lateral gene transfer dynamics in the ancient bacterial genus Streptomyces[END_REF], natural selection of specific substitutions also appears to play a crucial role in genome diversification of marine picocyanobacteria and to have driven their adaptation to specific environments. Indeed, in the time necessary for one gene to be gained, we found that 20 to 60 amino acid substitutions accumulate in any picocyanobacterial genome (as estimated based on terminal branches of the phylogeny, Table 1). This finding brings new evidence to support the "Maestro Microbe" model of bacterial genome evolution recently proposed by Larkin and Martiny [START_REF] Larkin | Microdiversity shapes the traits, niche space, and biogeography of microbial taxa[END_REF], which posits that some phenotypic traits, such as thermal preferences, evolve by progressive fitness optimization of protein sequences rather than gene gains and losses. This theory is

Genome diversification in marine picocyanobacteria

mainly based on the lack of specific genes that may explain trait differences between closely related organisms inhabiting distinct niches, and one of the best examples concerns Prochlorococcus clades colonizing temperate (HLI) and warm (HLII) environments [START_REF] Coleman | Genomic islands and the ecology and evolution of Prochlorococcus[END_REF][START_REF] Martiny | Phosphate acquisition genes in Prochlorococcus ecotypes: Evidence for genome-wide adaptation[END_REF][START_REF] Kettler | Genome diversification in marine picocyanobacteria Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF][START_REF] Larkin | Microdiversity shapes the traits, niche space, and biogeography of microbial taxa[END_REF], which were partly mixed on our tree based on gene content despite a clear phylogenetic separation based on core marker genes (Fig. 8). The sequencing of new Synechococcus genomes also allowed us to extend the Maestro Microbe hypothesis to Synechococcus thermotypes [START_REF] Zwirglmaier | Global phylogeography of marine Synechococcus and Prochlorococcus reveals a distinct partitioning of lineages among oceanic biomes[END_REF][START_REF] Pittera | Connecting thermal physiology and latitudinal niche partitioning in marine Synechococcus[END_REF], since particularly few genes were found to be specific to the cold-adapted clades I and IV (Supplementary Table S4). In contrast, our analysis of Synechococcus core proteins containing amino acid variants shared exclusively by all members of these cold thermotypes revealed potential candidates for adaptation to cold waters (Supplementary Table S9). A number of these core proteins target essential cell functions such as protein metabolism or carbon fixation and metabolism, suggesting that sequence variations of these proteins have an impact on their efficiency at different temperatures. We also identified proteins involved in carotene biosynthesis and the oxidative stress response, suggesting that these pathways are impacted by cold temperature in marine picocyanobacteria. Overall, while experimental testing is needed to validate the role of these substitutions in adaptation to cold environments, this analysis provides numerous strong candidates for such validation (Supplementary Table S9). The fact that all members of clades I and IV share specific variants of the three proteins involved in the β-carotene synthesis pathway (with e.g. >2% of the protein sequence comprising residues specific to these clades in CrtL-b) is particularly striking, since physiological experiments have shown that members of clades I and IV were able to maintain or increase their β-carotene:chlorophyll a ratio in response to cold stress, while this ratio decreased in strains representative of warm thermotypes [START_REF] Pittera | Connecting thermal physiology and latitudinal niche partitioning in marine Synechococcus[END_REF]. Thus, these substitutions might allow cells of the former clades to maintain β-carotene synthesis in cold conditions, resulting in a reduction of the cold-induced oxidative stress. Additionally, four proteins potentially involved in the response to oxidative stress were found to display variants specific to clades I and IV (Supplementary Table S9). In much the same way, a recent study identified two substitutions in genes encoding the two subunits of phycocyanin in Synechococcus between these cold-adapted clades and the warm-adapted clades II and III, which were also thought to be involved in adaptation to distinct thermal niches: RpcA G-43 and RpcB S-42 in the former clades and RpcA A-43 and RpcB N-42 in the latter [START_REF] Pittera | Adaptive thermostability of light-harvesting complexes in marine picocyanobacteria[END_REF]. It is worth noting that these genes were not detected by the stringent approach used here either because of the absence of the multi-copy cpcA gene in the CB0101 genome due to assembly issues or to a single exception among the newly sequenced genomes, the clade I strain PROS-9-1 having an RpcB S-42. Given that clades I and IV have diverged about 425 My ago (confidence interval 308-468 My), the most parsimonious explanations for these many shared substitutions would be either an adaptive convergence or an ancient homologous recombination between ancestors of these clades. In this context, it is interesting to note that mutations were found to arise in just a few generations in a clonal Prochlorococcus strain as an adaptation to selective conditions such as UV radiation [START_REF] Osburne | UV hyper-resistance in Prochlorococcus MED4 results from a single base pair deletion just upstream of an operon encoding nudix hydrolase and photolyase[END_REF], antibiotics [START_REF] Osburne | The spontaneous mutation frequencies of Prochlorococcus strains are commensurate with those of other bacteria[END_REF] or phage pressure (Avrani et al., 2011), emphasizing the role of such substitutions in short-term adaptation, although only a subset of these are kept in the long term.

Conclusions

Current clades of marine picocyanobacteria might be considered as survivors of a former set of "backbone" populations (as defined by [START_REF] Kashtan | Single-cell genomics reveals hundreds of coexisting subpopulations in wild Prochlorococcus[END_REF] that appeared hundreds of millions years ago, and then optimized their sequence, while eventually losing most of the genes that initially allowed niche colonization [START_REF] Lawrence | Gene transfer in bacteria: speciation without species?[END_REF][START_REF] Cohan | The origins of ecological diversity in prokaryotes[END_REF][START_REF] Polz | Horizontal gene transfer and the evolution of bacterial and archaeal population structure[END_REF][START_REF] Kashtan | Single-cell genomics reveals hundreds of coexisting subpopulations in wild Prochlorococcus[END_REF].

Genome diversification in marine picocyanobacteria

More recently, each of these clades further diversified into a number of new backbone populations, which correspond to the within-clade microdiversity recently described in Prochlorococcus and Synechococcus (see e.g. Martiny et al., 2009b;[START_REF] Kashtan | Single-cell genomics reveals hundreds of coexisting subpopulations in wild Prochlorococcus[END_REF][START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF]Larkin et al., 2016b). One explanation for the topology of the phylogenetic tree based on core proteins (short branches at the leaves of the tree and long branches at the base of clades, Fig. 8) would be the occurrence of periods of rapid diversification, as previously suggested for the occurrence of the different Synechococcus clades within SC 5.1 and of the Prochlorococcus radiation [START_REF] Urbach | Rapid diversification of marine picophytoplankton with dissimilar light-harvesting structures inferred from sequences of Prochlorococcus and Synechococcus (cyanobacteria)[END_REF][START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF]) and more extended periods during which each population stays relatively genetically homogeneous (e.g. by homologous recombination or by frequent genomic sweeps).

Alternatively, and perhaps more likely, picocyanobacterial populations might undergo continuous diversification at a fairly constant rate, with diversity purged during rare but severe extinction events, leaving traces only of the surviving ones. While it is tempting to relate these events (diversification or purge) to past geological and climatic shifts, this would need a more thorough examination with an improved time calibration.

One of the next challenges will be to more precisely relate variants (genes or substitutions) to a particular niche. We could advocate achieving this via comparative genomics, but this usually necessitates hundreds to thousands of closely related genomes (for review see [START_REF] Read | Characterizing the genetic basis of bacterial phenotypes using genome-wide association studies : a new direction for bacteriology[END_REF][START_REF] Chen | Diverse and unique picocyanobacteria in Chesapeake Bay , revealed by 16S-23S rRNA internal transcribed spacer sequences[END_REF], as well as a refined phenotypic characterization of the sequenced strains. Alternatively, one could search in situ data for genes or substitutions related to a particular niche or environmental parameter (see e.g. [START_REF] Kent | Global biogeography of Prochlorococcus genome diversity in the surface ocean[END_REF][START_REF] Grébert | Light color acclimation is a key process in the global ocean distribution of Synechococcus cyanobacteria[END_REF][START_REF] Ahlgren | Genomic mosaicism underlies the adaptation of marine Synechococcus ecotypes to distinct oceanic iron niches[END_REF][START_REF] Garcia | Linking regional shifts in microbial genome adaptation with surface ocean biogeochemistry[END_REF]. Given the wealth of marine metagenomes that are becoming available for a large variety of environmental niches, such an approach should be particularly powerful to unveil niche adaptation processes in the forthcoming years.

Methods

Genome sequencing and assembly

Thirty-four Synechococcus strains were chosen for genome sequencing based on their phylogenetic position, pigment content and isolation sites (Fig. 1, Supplementary Table S1). All but the three KORDI strains were retrieved from the Roscoff Culture Collection (RCC; http://roscoff-culturecollection.org/) and transferred three times on 0.3% SeaPlaque Agarose (Lonza, Switzerland) to clone them and reduce contamination by heterotrophic bacteria. A first set of 25 Synechococcus genomes (including WH8103) were generated at the Genoscope (CEA, Paris-Saclay, France) by shotgun sequencing of two libraries: a short-insert forward-reverse pair-end (PE) library (50-150 bp) and a long-insert reverse-forward mate-pair library (4-10 kb), both sequenced by Illumina™ technology. Additionally, seven other genomes were sequenced at the Center for Genomic Research (University of Liverpool, UK) by shotgun sequencing of 250 bp reads. Single or PE reads were first assembled into contigs using the CLC AssemblyCell © 4.10 (CLC Bio, Aarhus, Denmark). Synechococcus contigs were identified based on their different coverage compared to heterotrophic bacteria, scaffolded using WiseScaffolder and 28 out of 31 genomes were closed by manual finishing as described in [START_REF] Farrant | WiseScaffolder: an algorithm for the semi-automatic scaffolding of Next Generation Sequencing data[END_REF]. Three genomes (BIOS-E4-1, BOUM118 and RS9915), had only one to three gaps in highly repeated genomic regions. The base numbering of the circularized genomes was set at 174 bp before the dnaN start, corresponding approximately to the origin of replication. Automatic structural and functional annotation of the genomes was then realized using the Institute of Genome Science (IGS) Annotation Engine (http://ae.igs.umaryland.edu/cgi/index.cgi; [START_REF] Galens | The IGS standard operating procedure for automated prokaryotic annotation[END_REF]. As concerns Genome diversification in marine picocyanobacteria KORDI-49, KORDI-52 and KORDI-100 strains, genomes were sequenced from axenic cultures using a 454 GS-FLX Titanium sequencing system (Roche) at Macrogen (Seoul, Korea). The obtained reads were assembled using the Newbler assembler (version 2.3, Roche). To fill contig gaps, additional PCR and primer walking was conducted. Sequences of all new Synechococcus genomes were deposited in GenBank under accession numbers CP047931-CP047961 (BioProject PRJNA596899), except Synechococcus sp. WH8103 that was previously deposited to illustrate the performance of the pipeline used to assemble, scaffold and manually finish these genomes (Supplementary Table S1).

Clustering of orthologous genes

Protein and RNA sequences retrieved from new genomes were clustered with genomes previously available (Supplementary Table S1) into CLOGs using the OrthoMCL algorithm [START_REF] Li | OrthoMCL: Identification of ortholog groups for eukaryotic genomes[END_REF] and were then imported into the custom-designed Cyanorak v2.1 information system (www.sbroscoff.fr/cyanorak/) for further manual curation and functional annotation. Clustering into CLOGs allowed us to build phyletic patterns (i.e. the number of copies of each gene in each genome per CLOG), which was used to extract lists of genes shared at different taxonomic levels. Core genomes were defined at the genus, sub-cluster and clade levels when more than three genomes were available for a given taxonomic level (see Supplementary Table S2).

The phyletic pattern was also used to estimate the size of the pan-genome and core genome. The sampling of genome combinations necessary to draw pan-genome curves was performed with the software PanGP [START_REF] Zhao | PanGP: A tool for quickly analyzing bacterial pan-genome profile[END_REF] using as parameters 'Totally Random', SR=100 and SS=1000. Pangenome curves were then drawn with R custom designed scripts (v3.3.1.; [START_REF] Core | R: A language and environment for statistical computing[END_REF]. The results of PanGP exponential fits were used as estimates of the asymptotic number of core genes.

ANI/AAI calculation

Whole-genome ANI and percentage of conserved DNA between pairs of genomes (percentage of the genome length aligned by Blast with more than 90% ID) were calculated following the method described in [START_REF] Goris | DNA-DNA hybridization values and their relationship to whole-genome sequence similarities[END_REF]. AAI was calculated following the method described by (Konstantinidis and Tiedje, 2005b). When AAI values differed for a given pair of strains depending on which strain was used as a query for BLAST, the highest value was kept.

Phylogeny and tree comparisons

The petB phylogenetic tree was built using PhyML 3.1 [START_REF] Guindon | A simple, fast, and accurate algorithm to estimate large phylogenies by Maximum Likelihood[END_REF] with the HKY model and by estimating gamma parameters and the proportion of invariant sites, based on a database of 230 petB sequences [START_REF] Mazard | Multi-locus sequence analysis, taxonomic resolution and biogeography of marine Synechococcus[END_REF][START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF]. The confidence of branch points was determined by performing bootstrap analyses, including 1000 replicate data sets. Phylogenetic trees were edited using the Archaeopteryx v0.9901 beta program [START_REF] Han | PhyloXML: XML for evolutionary biology and comparative genomics[END_REF]. The tree was drawn using iTOL (http://itol.embl.de; [START_REF] Letunic | Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees[END_REF]. Additionally, a set of 821 single-copy core proteins were aligned with MAFFT v7.164b [START_REF] Katoh | MAFFT: Iterative refinement and additional methods[END_REF] and concatenated into a single alignment, resulting in a total of 226,778 amino acids. A phylogenetic tree was built with PhyML 3.1 with the WAG model and estimation of parameters of the gamma distribution and of the proportion Genome diversification in marine picocyanobacteria of invariant sites. The phylogeny based on gene content was performed as described in [START_REF] Wolf | Genome trees and the Tree of Life[END_REF]: a Jaccard distance matrix was computed from the phyletic pattern with the package vegan [START_REF] Oksanen | Vegan: Community ecology package. R package version 1[END_REF] and the matrix was then used by the Neighbor-Joining algorithm implemented in the R package ape [START_REF] Paradis | APE: Analyses of phylogenetics and evolution in R language[END_REF] to generate a tree with 100 bootstraps.

The phylogenetic tree based on core proteins was then compared to the tree based on the phyletic pattern using the R package dendextend v.1.3.0 [START_REF] Galili | dendextend: An R package for visualizing, adjusting and comparing trees of hierarchical clustering[END_REF]. Branch lengths were compared using custom python scripts based on the ete2 toolkit [START_REF] Huerta-Cepas | ETE: a python environment for tree exploration[END_REF]. Briefly, for each node at the base of a clade (highlighted by blue dots in Fig. 8), the average distance from the node to the descending leaves ('external' length) and the distance to the parent node ('internal' length) were calculated. Boxplots of the distribution of ratios of external to internal branch lengths were drawn in R for both trees and a paired Mann-Whitney-Wilcoxon test assessed the difference between the mean ratios.

Estimation of gene gains and losses

The number of gene gains and losses were assessed from phyletic patterns using the software Count [START_REF] Csurös | Count: Evolutionary analysis of phylogenetic profiles with parsimony and likelihood[END_REF] that implements a Maximum Likelihood method for estimating the ancestral states (presence, absence or multiple copies) of every CLOG in the dataset using the phylogenetic core protein tree as reference and allowing four categories for the gamma distribution of duplications and branch lengths (options -transfer_k 1 -length_k 4 -loss_k 1 -duplication_k 4). Cut-off on posterior probability was set at 90%, which allowed us to obtain 2,921 CLOGs at the root of the tree, a number similar to the average number of CLOGs in present-day Synechococcus strains. The state of presenceabsence of each gene family was then extracted at each node of the tree, and used to calculate the number of gene gains and losses on every branch. These estimations of gained genes were also used to predict genomic islands in each strain. A genomic island, starting and finishing with full-length gained genes, was defined from consecutive sliding windows (size 10,000 bp, interval 100 bp) with a ratio of nucleotides from gained CDS to total coding nucleotides higher than 50%. A network approach was then applied on all predicted islands to compare the gene content of these islands between all strains. Jaccard distances based on shared gene content were calculated between islands and an edge was drawn to connect two islands if their distance was higher than 0.1 (i.e., when two islands shared at least 10% of their pooled gene content). Network modules detection was then performed using the modularity algorithm [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]; resolution = 0.2) implemented in Gephi version 0.9.2 [START_REF] Bastian | Gephi: An open source software for exploring and manipulating networks[END_REF]. Furthermore, in order to take into account the phylogenetic relatedness between strains sharing genomic islands, a distance matrix based on core protein sequences was computed and used to color edges between nodes. Networks were then represented following the "Atlas 2" spatialization implemented in Gephi.

Time calibration of the tree

The core protein phylogeny was used as input for the reltime algorithm [START_REF] Tamura | Estimating divergence times in large molecular phylogenies[END_REF] and the JTT matrix-based model [START_REF] Jones | The rapid generation of mutation data matrices from protein sequences[END_REF], as implemented in MEGA7 [START_REF] Kumar | MEGA7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets[END_REF], with default parameters and SC 5.3 designated as an outgroup. Two calibration points were used, based on [START_REF] Sánchez-Baracaldo | Origin of marine planktonic cyanobacteria[END_REF] and TimeTree [START_REF] Kumar | TimeTree: A resource for timelines, timetrees, and divergence times[END_REF]: the first calibration point was set on node n2 (Supplementary Fig. S7), i.e. the common ancestor of strains WH5701 and WH8102 estimated Genome diversification in marine picocyanobacteria to have occurred between 582 and 878 My ago, and the second on node n4 (i.e. the common ancestor of strains CC9311 and WH8102; Supplementary Fig. S7), set between 252 and 486 My. This method allowed us to relate gain/loss events with the time elapsed on each branch of the tree, taking into account the higher evolution rate of protein-coding genes in Prochlorococcus than in Synechococcus [START_REF] Dufresne | Accelerated evolution associated with genome reduction in a free-living prokaryote[END_REF]. We also calculated the number of substitutions for each branch of the tree by multiplying branch length by the total number of residues in the alignment, and divided it by the time elapsed and the branch to obtain a substitution rate per My.

Estimation of the number of fixed genes and fixed substitutions specific to a taxon or shared between taxa

At a given node of the tree, genes that were found in all descending leaves and no other strain in the dataset were considered as fixed genes specific to this node. Similarly, every position that showed the same amino acid variant in all leaves below a node and another amino-acid in every other strain were considered as fixed variants specific to this node. Terminal branches were not taken into account in these calculations since, by definition, strain-specific amino acids or genes occurring in these branches cannot be considered as fixed.

Additionally, we also looked in Synechococcus-Cyanobium core genes for amino acid variants specific to a set of strains corresponding to clades (Supplementary Table S9). A variant was considered as specific to a set of strains if it showed the same amino acid in every strain within the set and any other amino acid in every other strain. To allow comparison between proteins of different lengths, the number of specific variants was normalized by gene length. Given that older clades are expected to have accumulated more substitutions, each set of strains proteins were ranked according to their proportion of specific variants. To identify candidate proteins potentially involved in adaptation to cold conditions in clades I and IV, we took the ratio of the protein rank for the "clades I and IV" set of strains to the median rank for other clades (excluding the clades containing a single strain). We kept only proteins for which this ratio was below 0.33, i.e. proteins with a rank 3 times higher in the "clades I and IV" set than in other clades (Supplementary Table S9).

Genome diversification in marine picocyanobacteria strain_subclade (or higher taxonomic level when no sub-clade has been defined). (B) Evolution of the pan and core genomes for an increasing number of picocyanobacterial genomes (red, 81 genomes), Synechococcus (orange, 53 genomes) and Prochlorococcus (green, 28 genomes). The grey zone around each curve represents the first and third quartiles around the median of 1,000 samplings by randomly modifying the order of genome integration. The ancestral state of presence/absence of every cluster of likely orthologous genes (CLOGs) was assessed using Count [START_REF] Csurös | Count: Evolutionary analysis of phylogenetic profiles with parsimony and likelihood[END_REF] and used to infer the number of gains and losses of gene families on each branch of the tree using the phylogenetic core protein tree as reference.

The number of gained and lost genes is labeled in blue and red, respectively. Nodes highlighted in red correspond to the major genome streamlining events that have occurred in the Prochlorococcus radiation.

Figure 6: Network of shared gene islands between all Prochlorococcus strains analyzed in this study. Each node corresponds to a genomic island in a given strain, the gene content of which is listed in Supplementary Table S5. Edges were colored according to the phylogenetic distance between strains, with red indicating closely-related strains and blue more distantly related strains, as indicated in the color bar. Edge width corresponds to the Jaccard distance between islands based on gene content. Nodes were colored based on Prochlorococcus clade. Modules cited in the text are surrounded with a grey line for those containing islands already described in the literature (subtitled with their names in [START_REF] Coleman | Genomic islands and the ecology and evolution of Prochlorococcus[END_REF][START_REF] Dufresne | Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria[END_REF]) and a black line for new modules described in the present study. The gene and genomic island composition of each module is described in Supplementary Table S6.

Figure 7: Same as Fig. 6 but for marine Synechococcus/Cyanobium strains. The gene and genomic island composition of each module is described in Supplementary Table S7. 15 Tables Table 1: Estimation of the number of gained, lost and/or fixed genes per million years (My) as well as total and fixed number of substitutions on internal branches (int. b.) or terminal branches (ter. b.) for Prochlorococcus HL and Synechococcus (Syn) SC 5.1. SE: standard error, adj. R 2 : adjusted R 2 . [20] and the background colors correspond to the finest possible taxonomic resolution obtained for each strain using the petB marker gene (left hand side legend). Colored circles surrounding the tree indicate newly sequenced genomes, while squares indicate previously available ones. Note that the WH8020 genome indicated by a diamond was not used in this study due to its poor quality. Symbols are colored according to their pigment type as defined previously ([137-139]; right hand side legend). Distribution of clusters of likely orthologous genes (CLOGs) in picocyanobacterial genomes. A CLOG is considered as core in a taxonomic group if it is present in ≥ 90% of the strains within this group. Sets of core CLOGS are inferred only for taxonomic groups with more than 3 genomes. Strains are labeled as strain_subclade (or higher taxonomic level when no sub-clade has been defined). B. Evolution of the pan and core genomes for an increasing number of picocyanobacterial genomes (red, 81 genomes), Synechococcus (orange, 53 genomes) and Prochlorococcus (green, 28 genomes). The grey zone around each curve represents the first and third quartiles around the median of 1,000 samplings by randomly modifying the order of genome integration.

Rate

Figure 5: Estimation of the gene gains and losses during the evolution of marine picocyanobacteria. The ancestral state of presence/absence of every cluster of likely orthologous genes (CLOGs) was assessed using Count [130] and used to infer the number of gains and losses of gene families on each branch of the tree using the phylogenetic core protein tree as reference. The number of gained and lost genes is labeled in blue and red, respectively. Nodes highlighted in red correspond to the major genome streamlining events that have occurred in the Prochlorococcus radiation. S5. Edges were colored according to the phylogenetic distance between strains, with red indicating closely-related strains and blue more distantly related strains. Edge width corresponds to the Jaccard distance between islands based on gene content. Nodes were colored based on Prochlorococcus clade. Modules cited in the text are surrounded with a grey line for those containing islands already described in the literature (subtitled with their names in [51] and [22]) and a black line for new modules described in the present study. The gene and genomic island composition of each module is described in Additional file 1: Table S6. S7.

Figure 5 :

 5 Figure5: Estimation of the gene gains and losses during the evolution of marine picocyanobacteria. The ancestral state of presence/absence of every cluster of likely orthologous genes (CLOGs) was assessed using Count[START_REF] Csurös | Count: Evolutionary analysis of phylogenetic profiles with parsimony and likelihood[END_REF] and used to infer the number of gains and losses of gene families on each branch of the tree using the phylogenetic core protein tree as reference. The number of gained and lost genes is labeled in blue and red, respectively. Nodes highlighted in red correspond to the major genome streamlining events that have occurred in the Prochlorococcus radiation.

Figure 8 :

 8 Figure 8: Comparison of phylogenies based on core protein sequences and phyletic patterns of non-core genes. Left, Maximum Likelihood tree based on the alignment of 821 concatenated core proteins. Right, Neighbor-Joining tree based on the Jaccard distance between the phyletic patterns of 27,376 accessory gene families found in the 81 picocyanobacterial genomes. Labels are colored according to the strain sub-clade. Red branches indicate discrepancies between the topology of the two trees. Nodes located at the base of a clade and highlighted by blue dots were used for branch length comparisons in Supplementary Fig. S4.

Figure 1 :

 1 Figure 1: Phylogenetic position of the 53 marine Synechococcus genomes used in this study. A maximum-Likelihood tree was generated based on 231 petB marine Synechococcus sequences from both cultured and environmental samples. Diamonds indicate bootstrap support over 70%. Sequences were named after strain name_sub-cluster_clade_subclade (sub-clade assignments as in Farrant et al.[20] and the background colors correspond to the finest possible taxonomic resolution obtained for each strain using the petB marker gene (left hand side legend). Colored circles surrounding the tree indicate newly sequenced genomes, while squares indicate previously available ones. Note that the WH8020 genome indicated by a diamond was not used in this study due to its poor quality. Symbols are colored according to their pigment type as defined previously([137-139]; right hand side legend).

Figure 2 :

 2 Figure 2: Relationship between genome size and GC3% (GC content at the third codon position). Each symbol corresponds to a different genome, with Prochlorococcus indicated by circles and Synechococcus by triangles. The color of each symbol indicates the clade or SC.

Figure 3 :

 3 Figure3: Genomic diversity of marine picocyanobacteria. A. Heatmap of average nucleotide identity (ANI, bottom left triangle) and average amino acid identity (AAI, upper right triangle) between pairs of genomes. Each lane corresponds to a strain, and strains are ordered according to their phylogenetic relatedness. Strains are as labeled as strain_subclade (or higher taxonomic level when no sub-clade has been defined). B. Relationships between 16S rRNA identity, AAI, and taxonomic information for Synechococcus (left panel) and Prochlorococcus (right panel) genomes. Dots correspond to comparisons between pairs of genomes belonging to the same clade, triangles between pairs of genomes belonging to the same SC but different clades and squares between pairs of genomes belonging to different SC.

Figure 4 :

 4 Figure 4: Core, accessory and pan genomes of marine picocyanobacteria. A. Distribution of clusters of likely orthologous genes(CLOGs) in picocyanobacterial genomes. A CLOG is considered as core in a taxonomic group if it is present in ≥ 90% of the strains within this group. Sets of core CLOGS are inferred only for taxonomic groups with more than 3 genomes. Strains are labeled as strain_subclade (or higher taxonomic level when no sub-clade has been defined). B. Evolution of the pan and core genomes for an increasing number of picocyanobacterial genomes (red, 81 genomes), Synechococcus (orange, 53 genomes) and Prochlorococcus (green, 28 genomes). The grey zone around each curve represents the first and third quartiles around the median of 1,000 samplings by randomly modifying the order of genome integration.

Figure 6 :

 6 Figure6: Network of shared gene islands between all Prochlorococcus strains analyzed in this study. Each node corresponds to a genomic island in a given strain, the gene content of which is listed in Additional file 1: TableS5. Edges were colored according to the phylogenetic distance between strains, with red indicating closely-related strains and blue more distantly related strains. Edge width corresponds to the Jaccard distance between islands based on gene content. Nodes were colored based on Prochlorococcus clade. Modules cited in the text are surrounded with a grey line for those containing islands already described in the literature (subtitled with their names in[51] and [22]) and a black line for new modules described in the present study. The gene and genomic island composition of each module is described in Additional file 1: TableS6.

Figure 7 :

 7 Figure 7: Same as Fig. 6 but for marine Synechococcus/Cyanobium strains. The gene and genomic island composition of each module is described in Additional file 1: TableS7.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Author Contributions FH, MR, FDP, MO, DHC, JHN and CS purified newly sequenced Synechococcus strains, FH, MR, FDP, DHC, JHN and MO extracted the DNA, KL, JMA, PW, FDP, DHC, JHN, FP, LG, MH and GKF participated in sequencing and/or assembly of the genomes, MH, GLC, EC, AB, LBG and GFK developed and run the automatic clustering and annotation pipelines. HD, UG, FP and LG participated in the expert manual annotation of the genomes. HD, GKF, UG and LG generated and processed the data. UG, HD, JH and DE produced the genomic island networks. HD, GKF, UG, DE, DJS, FP and LG analyzed the results. HD, UG, JH, GKF and LG made the figures. All authors contributed to the preparation of the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by the French "Agence Nationale de la Recherche" Programs SAMOSA (ANR-13-ADAP-0010) and CINNAMON (ANR-17-CE2-0014-01), the Genoscope project METASYN and Natural Environment Research Council grant NE/I00985X/1.

Acknowledgments

We would like to thank the Institute for Genome Sciences Annotation Engine service at the University of Maryland School of Medicine, and in particular Michelle Giglio and Suvarna Nadendla, for providing automatic structural and functional annotation of the sequences, Brian Palenik and Tanja Woike for authorizing us to use the two unpublished Synechococcus genomes WH8016 and CC9616 as well as Garance Monier and Théo Sciandra for participating in the curation of the Cyanorak v2.1 database. We warmly thank the Roscoff Culture Collection and Sophie Mazard for maintaining and isolating some of the Synechococcus strains used in this study as well as the ABiMS platform for providing computational support for this work. This work is dedicated to our esteemed colleague Christophe Caron, who deceased on May 5th 2018.

Genome diversification in marine picocyanobacteria

Supplementary Material 13 Data Availability Statement

The datasets generated and analyzed in this study can be found in the Genbank repository (https://www.ncbi.nlm.nih.gov/genbank/). NCBI accession number of each genome is available in Supplementary Table S1. Sequences were named after strain name_sub-cluster_clade_subclade (sub-clade assignments as in [START_REF] Farrant | Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria[END_REF] and the background colors correspond to the finest possible taxonomic resolution obtained for each strain using the petB marker gene (left hand side legend). Colored circles surrounding the tree indicate newly sequenced genomes, while squares indicate previously available ones. Note that the WH8020 genome indicated by a diamond was not used in this study due to its poor quality. Symbols are colored according to their pigment type as defined previously (Humily et al., 2014;Xia et al., 2017b;[START_REF] Grébert | Light color acclimation is a key process in the global ocean distribution of Synechococcus cyanobacteria[END_REF]; right hand side legend). S3 of [54] and Supplementary Material 5 of [22] are indicated in blue and red, respectively. Abbreviations: ISL and SVR correspond to 'islands' and 'smaller variable regions', respectively as defined in previous work; GI, genomic islands, as defined in the present work.

Figure Legends