David Demory
email: david.demory@biosci.gatech.edu

Anne-Claire Baudoux

Adam Monier

Nathalie Simon

Christophe Six

Pei Ge

Fabienne Rigaut-Jalabert

Dominique Marie

Antoine Sciandra

Olivier Bernard
email: olivier.bernard@inria.fr

Sophie Rabouille
email: srabouille@obs-vlfr.fr

Picoeukaryotes of the Micromonas genus: sentinels of a warming ocean

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The Intergovernmental Panel for Climate Change (IPCC) stressed unequivocal warming of the climate system. Their Fifth Report anticipates rises in the global mean surface temperature by the end of 21st century ranging from 0.3-1.7 • C (RCP2.6) to 2.6-4.8 • C (RCP8.5) [1]. Oceans participate in buffering the increasing emissions of greenhouse gases, thus modulating the warming; in addition to the chemical equilibration of gas species between the atmosphere and dissolved phases, phytoplankton is an important contributor of carbon remediation through CO 2 sequestration in the ocean [2]. Should dramatic shifts occur in species biodiversity and distribution following temperature increases [3,4], the resilience of ecosystems could severely be impaired. The likely responses of ecosystems to such rapid temperature changes are at the core of debates, with worrisome consequent impacts on oceanic biogeochemical cycles and feedbacks on the climate system [5].

Phytoplankton live in a thermally fluctuating environment that constrains growth capacity [3,6,7]. The temperature growth response of phytoplankton varies widely, both between and within taxa. Phenotypic plasticity determines the ability to acclimate to short-term environmental variations while genetic adaptations characterize evolutionary processes under long-term changes. These features will provide, or not, each species with the capacity to survive in a given biotope and to evolve by modifying their thermal niche. Since temperature depends on latitude [8,9,10,11], it is therefore a probable driver of niche partition in the oceans, creating large-scale biogeographic patterns [12]. Hence, the structure and diversity of phytoplankton communities could partly reflect observed trends in the global temperature [6,7].

Temperature-related interspecific distributions have been studied for the whole phytoplankton community [3] but few studies explored intragenus diversity [START_REF] Pittera | Connecting thermal physiology and latitudinal niche partitioning in marine Synechococcus[END_REF][START_REF] Martiny | Interactions between Thermal Acclimatation, Growth Rate, and Phylogeny Influence Prochlorococcus Elemental Stoichiometry[END_REF]. Micromonas species have emerged as emblematic representative of the 2 Submitted to the ISME Journal eukaryotic pico-phytoplankton communities, thriving in a variety of ecosystems from polar to tropical waters [START_REF] Lovejoy | Distribution, phylogeny, and growth of cold-adapted picoprasinophytes in Arctic seas[END_REF][START_REF] Worden | Ecology and diversity of picoeukaryotes. Microbial Ecology of the Oceans[END_REF][START_REF] Monier | Gene invasion in distant eukaryotic lineages: discovery of mutually exclusive genetic elements reveals marine biodiversity[END_REF][START_REF] Monier | Oceanographic structure drives the assembly processes of microbial eukaryotic communities[END_REF]. They often dominate phytoplankton in coastal environments [START_REF] Not | A Single Species, Micromonas pusilla (Prasinophyceae), Dominates the Eukaryotic Picoplankton in the Western English Channel[END_REF], where their major contribution to primary production influences the biogeochemical cycles [START_REF] Massana | Eukaryotic Picoplankton in Surface Oceans[END_REF]. In the past decade, phylogenetic analyses identified several distinct genetic lineages within Micromonas and have suggested that this genus was composed of cryptic species [START_REF] Guillou | Diversity of Picoplanktonic Prasinophytes Assessed by Direct Nuclear SSU RDNA Sequencing of Environmental Samples and Novel Isolates Retrieved From Oceanic and Coastal Marine Ecosystems[END_REF][START_REF] Slapeta | Global Dispersal and Ancient Cryptic Species in the Smallest Marine Eukaryotes[END_REF][START_REF] Foulon | Ecological niche partitioning in the picoplanktonic green alga Micromonas pusilla: evidence from environmental surveys using phylogenetic probes[END_REF][START_REF] Marin | Molecular phylogeny and classification of the Mamiellophyceae class. nov. (Chlorophyta) based on sequence comparisons of the nuclear-and plastid-encoded rRNA operons[END_REF]. Four species have now been formally described [START_REF] Simon | Revision of the genus Micromonas Manton et Parke (Chlorophyta, Mamiellophyceae), of the type species M. pusilla (Butcher) Manton & Parke and of the species M. commoda van Baren, Bachy and Worden and description of two new species based on the genetic and phenotypic characterization of cultured isolates[END_REF]. Micromonas spp. may co-occur at various latitudes, but were found to occupy different temporal or depth niches within their sympatric ranges [START_REF] Foulon | Ecological niche partitioning in the picoplanktonic green alga Micromonas pusilla: evidence from environmental surveys using phylogenetic probes[END_REF].

As observed for picocyanobacteria [START_REF] Johnson | Niche partitioning among Prochlorococcus ecotypes along ocean-scale environmental gradients[END_REF][START_REF] Pittera | Connecting thermal physiology and latitudinal niche partitioning in marine Synechococcus[END_REF], the temperature response of such a widely distributed and phylogenetically diverse eukaryote is expected to vary between Micromonas species. The interspecific diversity within the genus Micromonas, the number of characterized strains, and abundant omics data make it a relevant model organism to both explore the impact of temperature on latitudinal distribution and diversity of phytoplankton, and to shed light on the mechanisms that drive phytoplankton thermal responses in the ocean. We therefore studied the thermotolerance and thermal growth response of eleven Micromonas strains in the laboratory under controlled conditions (hereafter referred as experimental strains) and we derived a mathematical model that describes the impact of temperature on growth rate. With this model, we uncover the logic that lies behind the observed distribution of species and their co-occurrence; we also reveal the existence of thermotypes within the genus. We extrapolated the thermal response to a set of 46 additional strains from the Roscoff Culture Collection (hereafter referred as collection strains), observed in various oceanic regions, showing that temperature is the main driver of diversity and distribution in this genus. Then, we developed a predictive model of niche partition to characterize Micromonas interspecific diversity, which we successfully validated against the Tara Oceans dataset [START_REF] Vargas | Eukaryotic plankton diversity in the sunlit ocean[END_REF], making it a plausible prediction tool. We demonstrated that Micromonas distribution is a relevant and accurate proxy of the whole phytoplankton community distribution.

More than a sentinel of the ocean biogeochemistry as previously suggested by Worden and colleagues [START_REF] Worden | Green evolution and dynamic adaptations revealed by genomes of the marine picoeukaryotes Micromonas[END_REF], Micromonas is a probe for global warming. To explore how phytoplankton communities may respond to a future, warmer ocean, we ran the niche partition model under IPCC Sea Surface Temperature (SST) projections, adding an evolutionary model that accounts for the potential adaptation of growth to temperature changes.

Results and discussion

Micromonas strains feature distinct physiological responses to temperature To estimate the temperature tolerance and growth responses of the four described Micromonas species, we selected three strains of M. commoda, M. bravo and M. pusilla as well as two strains of M. polaris. We measured their exponential growth rate after being grown for two months between 4 2) depending on the strain origin.

To increase the accuracy in the temperature response estimation, the experimental protocols followed the recommendations given in [START_REF] Boyd | Marine phytoplankton temperature versus growth responses from polar to tropical waters-outcome of a scientific community wide study[END_REF][START_REF] Grimaud | Modeling the temperature effect on the specific growth rate of phytoplankton: a review[END_REF] (see Methods). The chosen strains, obtained from the Roscoff Culture Collection (RCC), were originally isolated from contrasted thermal niches of the Atlantic, Pacific and Artic basins (Figure 1a, Supplementary Figure 3 and Supplementary Table 1). All showed a typical [START_REF] Raven | Temperature and Algal Growth[END_REF][START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF] asymmetric growth response to temperature, which we characterized by four cardinal growth parameters: T min and T max , respectively the minimum and maximum temperatures for growth; µ opt , the maximum specific growth rate obtained at the optimum temperature T opt (Figure 1b). Overall, the Micromonas genus was able to grow over the thermal range tested, but with diverse and specific responses for each strain, depicted by distinct cardinal parameters (Supplementary Table 6). Temperature stimulates enzymatic processes and metabolic rates, but also accelerates cell mortality [START_REF] Serra-Maia | Influence of temperature on Chlorella vulgaris growth and mortality rates in a photobioreactor[END_REF]. In the suboptimal range (T < T opt ), enzymatic activity increases more than mortality in response to increasing temperatures. At T opt this balance between metabolic activity and mortality is optimized and yields the highest observed net growth rate. At supra optimal temperatures (T > T opt ), the denaturation of key metabolic enzymes, like rubisco [START_REF] Salvucci | Inhibition of photosynthesis by heat stress: The activation state of Rubisco as a limiting factor in photosynthesis[END_REF] and the thermolability of Photosystem II [START_REF] Rokka | Dephosphorylation of photosystem II reaction center proteins in plant photosynthetic membranes as an immediate response to abrupt elevation of temperature[END_REF] are exacerbated, along with an increase of the membrane damages [START_REF] Da | Membrane fluidity and its roles in the perception of environmental signals[END_REF] ; as a consequence, the net growth rate sharply decreases with temperature up to the maximal growth temperature the strain can withstand (T max at which µ is null).

Several patterns appeared when comparing the growth response to the annual average SST (T S ) at the site where each strain was isolated. Strains isolated in locations where T S was above 19.7 • C (RCC 299 and RCC 829) were able to grow up to high temperatures (T max = 32.6±0.02 and 37.0±0.12 • C, respectively); they showed a high µ opt (1.1±0.05 to 1.3±0.07d -1 , respectively) at an elevated opti- In Summary, the four formerly described Micromonas species exhibited specific temperature tolerance and growth optima in vitro and their according response parameters were related to the thermal environment from which the strains were isolated. Model parameters T min , and to a lesser extent T max , are difficult to accurately estimate [START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF]. Since measurements for temperatures close to T min (but slightly higher) and close to T max (but slightly lower) are generally rare, they must be extrapolated from a mathematical model. These parameters also bracket the thermal niche, i.e. the breadth of the thermal response. For instance, it appears that Arctic strains showed a much narrower niche: they were more stenotherm compared to the other strains.

The Micromonas genus includes six thermotypes: evidence from the most recent phylogeny.

The phylogenetic analysis of the 57 Micromonas 18S DNA sequences from the eleven experimental and 46 collection strains highlighted the existence of six distinct phylogenetic groups (see Methods and Supplementary Figure 2). To identify whether they were associated with specific thermal conditions in the ocean, we analyzed available data of average SST in areas where Micromonas spp.

were sampled. We computed a non-metric dimensional scaling (NMDS) of the thermal environment dataset (see Methods). The significant ordination (stress = 0.004) identified six different distributions in the thermal environment, from warmer, low latitudes to colder, high latitudes, that showed a good match with the phylogenetic tree (Figure 2a and Supplementary Figure 3), demonstrating that the thermal niche of Micromonas was related to its phylogenetic affiliation. M. polaris and M. pusilla strains occupied respectively a narrow and wide thermal niche while M. bravo and M. commoda each included two distinct groups. One isolated from a warmer (lower latitude; warm group) and one isolated in a colder (higher latitudes; cold group) environment (Figure 2a and Supplementary Figure 3).

There are few examples in the literature of latitudinal segregation within eukaryotic phytoplankton genera [START_REF] Cuvelier | Targeted metagenomics and ecology of globally important uncultured eukaryotic phytoplankton[END_REF][START_REF] Monier | Phylogenetic diversity and biogeography of the Mamiellophyceae lineage of eukaryotic phytoplankton across the oceans[END_REF]. For example, the global distribution of Ostreococus clades, a picoeukaryote close to Micromonas is related to temperature but first seems to discriminate rather coastal, high-light adapted clades from more oceanic, low-light adapted clades [START_REF] Demir-Hilton | Global distribution patterns of distinct clades of the photosynthetic picoeukaryote Ostreococcus[END_REF]. In agreement with the hypothesis of Submitted to the ISME Journal Foulon et al. [START_REF] Foulon | Ecological niche partitioning in the picoplanktonic green alga Micromonas pusilla: evidence from environmental surveys using phylogenetic probes[END_REF], our experimental and phylogenetic results showed that a niche segregation within Micromonas did occur that is consequent to thermal, groupspecificities and which compels with the recently identified, four known species.

The present analysis further revealed the existence of two thermotypes within both M. commoda and M. bravo species, making a total of six distinct Micromonas thermotypes.

Establishing a thermal response model for Micromonas thermotypes

To obtain a better appraisal of the thermal response of strains, we looked for possible correlations between cardinal growth parameters and environmental features where strains had been isolated. Among the tested descriptors of the SST dynamics, the average surface temperature at the isolation site (T S ) best correlated with the cardinal temperature. For T min , the latitude was also included in the regression (Table 1 and Supplementary Figure 6a). The optimal growth rate (µ opt ) increased with T S , following the Eppley's hypothesis of a faster growth rate at warmer temperatures [START_REF] Eppley | Temperature and phytoplankton growth in the sea[END_REF]. The maximal growth temperature (T max ) and the optimal growth temperature (T opt ) were also both positively correlated with T S , suggesting that environmental temperature featured the upper tolerance window of strains. The minimal temperature of growth (T min ) had the lowest correlation with the environmental temperature (Supplementary Figure 6b), as also reported by [START_REF] Chen | Patterns of thermal limits of phytoplankton[END_REF] for different phytoplankton species. We found that the minimal growth temperature T min best correlated (negatively) with a combination of the yearly average temperature T S and latitude (Lat, Supplementary Figure 7). In the end, the growth response (µ opt , T min , T opt and T max ) of cultured strains can thus accurately be predicted from the thermal environments (T S ) and latitude from which they were isolated, using the relations defined in Table 1. Last, statistically significant correlations were also found between cardinal parameters (Supplementary Figure 8). In particular, the optimal temperature of growth (T opt ) linearly correlated with the maximal temperature of growth (T max ) by a factor close to 1, as previously highlighted for a wide range of bacterial species [START_REF] Rosso | An unexpected correlation between cardinal temperatures of microbial growth highlighted by a new model[END_REF].

The relationships between cardinal growth parameters and environmental temperatures deduced from the culture experiments (Table 1) were used to extrapolate the cardinal parameters of 46 additional Micromonas collection strains, using the latitude and average annual temperature of their isolation site (Table 1 and Supplementary Table 9). This data set confirmed a segregation of the four species into six different thermotypes. To deduce a representative thermal response for each thermotype, we randomly chose 100,000 values within the confidence interval of the cardinal parameters of each group and ran Monte Carlo simulations of the related thermal responses (see Methods). The Bernard and Rémond (BR) model was then fitted to each bundle of simulated responses [START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF] to obtain the average thermal response curve representative of each thermotype (Figure 2b and Supplementary Figures 9,10 and 11). Last, we calibrated the envelope curve, inspired from [START_REF] Grimaud | Modelling the temperature effect on phytoplankton: from acclimation to adaptation[END_REF], on the Micromonas genus, by fitting the BR model [START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF] to the set of (T opt ,µ opt ) obtained for each thermotype (see Methods and Figure 2b ).

With the narrowest thermal niche (23.04 ± 2.42 Tara Oceans dataset validates the global segregation of thermotypes.

To validate our hypothesis that temperature is a key factor that greatly influences Micromonas biogeography over a yearly period, we retrieved the 18S V9 metabarcodes dataset obtained in the frame of Tara Oceans [START_REF] Vargas | Eukaryotic plankton diversity in the sunlit ocean[END_REF](Figure 3). Read abundance data assigned to each of the Micromonas thermotypes were identified across 47 stations, spanning 6 marine regions with different thermal environments: Mediterranean Sea, Red Sea, Indian Ocean, South Atlantic Ocean, Southern Ocean and South Pacific Ocean (Figure 3a). Using an NMDS ordination method, we first compared the relative abundance of each Micromonas thermotype at sampling stations (see Methods and Figure 3b) to the physicochemical environmental conditions observed along the Tara Oceans circumnavigation. The presence of Micromonas species was better explained by temperature (R 2 = 0.48, p-value < 0.001) than by nutrient availability, mixing, or geographical location. To a lower extent, nutrients (NO 2 +NO 3 , PO 4 and NO 2 ; R 2 < 0.23, p-value <0.032), Chl a concentration (R 2 = 0.1710, p-value = 0.003) and mixed layer depth (MLD;R 2 = 0.13, p-value = 0.03) also explained significantly the Micromonas assemblages along the transect. Temperature is thus the strongest descriptor of the change in diversity between Tara Oceans stations.

Submitted to the ISME Journal

We then compared the relative abundance of thermotypes at all stations in relation to yearly SST (Figure 3c). A very clear thermal separation appeared between the two M. commoda thermotypes, further supporting our identification of two distinct thermotypes. M. commoda cold was most abundant in waters with temperature below 20 • C and rarely found beyond 25 • C, while M. commoda warm mostly occurred between 25 • C and 30 • C and was completely absent at stations where temperatures were below 15 • C. Species M. bravo was less often observed than M. commoda and showed overlapping distributions of its two warm and cold thermotypes, which we believe was due to the large thermal niche of the warm thermotype spreading over that of the more restrained, cold thermotype (Figure 2).

A non-distinct distribution (Figure 3c) in the Tara Oceans data could also suggest that the evolution of the two M. bravo thermotypes was more recent. Species M. polaris was observed only at stations with T<10 • C with highest abundances near 0 • C, validating the psychrophilic characteristics of this thermotype. Species M. pusilla was only found at a few stations compared to M. commoda and M. bravo; it was observed from 12 • C to 30 • C with a maximum abundance above 25 • C. This distribution may well be related to the fact that its thermal response is close to the barycenter of the whole Micromonas thermal response (average parameters:

T opt = 21.26, µ opt = 0.84 and (T max -T min ) = 28.34). The reported occurrences of this species at low concentrations all around the globe [START_REF] Foulon | Ecological niche partitioning in the picoplanktonic green alga Micromonas pusilla: evidence from environmental surveys using phylogenetic probes[END_REF][START_REF] Baudoux | Interplay between the genetic clades of Micromonas and their viruses in the Western English Channel[END_REF] could support the idea that it plays a "seed bank" role, acting as a dormancy stage of Micromonas compared to other species [START_REF] Lennon | Microbial seed banks: the ecological and evolutionary implications of dormancy[END_REF]. Interestingly, Foulon et al. [START_REF] Foulon | Ecological niche partitioning in the picoplanktonic green alga Micromonas pusilla: evidence from environmental surveys using phylogenetic probes[END_REF] also suggested a possible niche partition over depth, along a light gradient that may explain the low concentration of M. pusilla in the Tara Oceans dataset. In the end, temperature is a sufficient parameter to describe the latitudinal segregation of Micromonas between Tara Oceans stations. The current typology of Tara Oceans (they mainly are open ocean areas), does not allow to fully assess a possible effect of nutrients [START_REF] Botebol | Acclimatation of a low iron adapted Ostreococcus strain to iron limitation through cell biomass lowering[END_REF].

Influence of temperature on the intragenus diversity of Micromonas assemblages

To further understand the thermal niche partition of Micromonas at the global scale, we proposed a simple index to relate Micromonas intragenus diversity to the global average SST (Figure 4 and Supplementary Figure 12). We computed an interspecific Micromonas diversity index (Shannon derived/based) from the growth response of a given thermotype i to a considered local temperature T Submitted to the ISME Journal according to the equation:

H(T) = n ∑ i=1 D i ln(D i ) with D i = µ i (T) ∑ n i=1 µ opt,i (1)
Where D i is the distribution index, µ i (T) is the growth rate at the temperature T, and µ opt,i is the optimal growth for the thermotype i. We compared H(T) to a Shannon-like index for the Micromonas genus at each Tara Oceans sampling station using the proportion of each Micromonas thermotype OTU in the total counted

Micromonas OTU and the local SST annual average (Figure 4a andb). Based on the calculated diversity index H(T), we were able to qualitatively predict the Micromonas intragenus diversity estimated from the Tara Oceans V9-18S dataset (Spearman test: ρ = 0.417, p-value = 0.0035), thereby validating our theoretical developments. The diversity index followed a fluctuating trend through the cruise path characterized by different thermal environments (Figure 3a).

When running the Micromonas diversity model at the global scale (Figure 4c and Supplementary Figure 13), the predicted diversity was minimal at the poles (Lat > 60 • N and >50 • S) and at the equator (between 20 • N and 20 • S), especially in the Indian Ocean and the Pacific Ocean (Figure 4c). Maximum diversity levels were found from 20 to 60 • N and from 20 to 40 • S. We used the relationship between the phytoplankton diversity as calculated by Thomas and colleagues [3] and our

Micromonas diversity to normalize our diversity index within Thomas's scale (see Methods). Our simulated global Micromonas diversity was point by point compared to the whole phytoplankton potential diversity calculated by Thomas and collaborators [3] (Figure 4d). We found a very strong relationship between the two diversity patterns (R 2 adj = 0.97, p-value < 0.05; see Methods and Supplementary Figure 15). This result strongly suggests that the diversity between Micromonas thermotypes, at mesoscale and on a yearly basis, is representative of the whole phytoplankton community. It likely explains the overall success of the genus to colonize very contrasted biotopes [START_REF] Not | A Single Species, Micromonas pusilla (Prasinophyceae), Dominates the Eukaryotic Picoplankton in the Western English Channel[END_REF][START_REF] Foulon | Ecological niche partitioning in the picoplanktonic green alga Micromonas pusilla: evidence from environmental surveys using phylogenetic probes[END_REF]. Micromonas could thus serve as a relevant marker of the biodiversity of phytoplankton communities. The term "sentinel", originally proposed by [START_REF] Worden | Green evolution and dynamic adaptations revealed by genomes of the marine picoeukaryotes Micromonas[END_REF] to depict the role of Micromonas on ocean biogeochemistry is all the more relevant considering this genus reflects the pattern of the whole phytoplanktonic community and can help to better anticipate the impact of ocean warming.

Submitted to the ISME Journal

Diversity evolution in a warmer ocean: a matter of the adaptation time scale

To explore the impact of future temperature changes on phytoplankton diversity, we investigated its evolution using SST projections over the period 2001-2100.

To account for the adaptation capability [4,[START_REF] Grimaud | Modelling the effect of temperature on phytoplankton growth across the global ocean[END_REF], we proposed a very simple adaptive model. This model assumes that the evolution time scale is related to the local doubling time ln(2)

µ i (T) of each thermotype i. Adaptation is thus faster for the warm thermotypes in warm environments. The adaptation dynamics describes the evolution of the cardinal temperatures (T min , T opt and T max ) from their present value to their value at the end of the century. The evolution rate is estimated according to the characteristic number of generations Na required to adapt to a different temperature, i.e. to shift each cardinal (i.e. represented by the character "c") temperature T c to its asymptotic value T * c , defined as the evolutionary equilibrium given the changes of the surface temperature T at each time step (Supplementary Figure 16). The evolution dynamics of each cardinal parameter T c,i is described by a simple first order equation:

dT c,i dt = N i (T(t)) Na (T * c,i (T(t)) -T c,i (t)) (2)
Where N i (T(t)) =

µ i (T(t),T c,i ) ln (2)
, with µ i (T(t), T c,i (t)) the growth rate at the temperature T, calculated using the set of cardinal parameters T c,i for the thermotype i.

We ran this model for different Na, from fast adaptation scales (Na < 100 generations) to slow adaptation scales (Na = 10 6 generations) and calculated the evolution of thermotypes diversity between the present period (2001 to 2010) and future period (2091 to 2100, Figure 5). We considered two realistic evolution hypotheses to describe the dichotomy between specialist and generalist species:

the Specialist-generalist hypothesis with constant thermal niche width (Figure 5a and b) and the Specialist-generalist hypothesis with dynamical thermal niche [START_REF] Stinchcombe | Genetics and evolution of function-valued traits: understanding environmentally responsive phenotypes[END_REF] (Figure 5c and d -see Methods). Over the 21st century, SST will globally increase by 2 to 3 • C over the whole ocean surface and up to 5 • C around 45 • N, with the exception of the highest latitudes, which may see a slight decrease in their average temperature (Supplementary Figure 17).

Similar erosion patterns were found for both specialist-generalist hypotheses that showed diversity losses between 40 • S to 40 • N. At latitudes higher than 40 • , we found possible gains in biodiversity, regardless of the adaptation scenario and the evolution hypothesis. At these latitudes, for most phytoplankton species, the optimum temperature (T opt ) is higher than the average environmental temperature Submitted to the ISME Journal (T S ). With a fast adaptation scenario, thermal traits follow the thermal environment and T opt remains above T S , each thermotype keeps its thermal niche and diversity is not affected. In contrast, thermal traits will not change fast enough in a slow adaptation scenario; T S gets closer to T opt , and each thermotype ends up with a fitness that is out of phase with the thermal environment (Supplementary Figure 18). While these conditions are still favorable for growth, they typically increase the diversity. Finally, for the adaptation scenario where the thermal niche can increase, it gives more chance for a species to adapt faster even for a higher change in the thermal environment (Supplementary Figure 19). At latitudes lower than 40 • , ocean warming will drive a decrease in phytoplankton diversity, with a mitigation of diversity losses tightly dependent on the adaptation time scale and similar for both hypotheses (Figure 5a andc). Slow adaptation scenarios lead to an important diversity erosion compared to fast adaptation scenarios, suggesting that the adaptation time scale is a key parameter in the mitigation of diversity loss and matters far more than the strategy of adaptation itself. In areas most vulnerable to diversity erosion (Supplementary Figure 20

Methods

A graphic abstract of the overall, scientific approach is provided in Supplementary Figure 1.

Growth measurements and thermal response model

Culture conditions. Eleven Micromonas spp. strains were selected from the RCC for the laboratory experiments. We chose strains representative of all the currently known species and according to their isolation site, to consider a range of organisms found along a latitudinal gradient (Supplementary Table 1). Cells were grown in batch cultures in ventilated polystyrene flasks (Nalgene, Rochester, NY, USA) in K-Si medium [START_REF] Keller | Media for the culture of oceanic ultraphytoplankton[END_REF]. Cultures were maintained in temperature-controlled chambers (Aqualytic, Dortmund, Germany) at different temperatures (4, 7.5, 9.5, 12.5, 20, 25, 27.5, 30 and 32.5 • C) for two months (see Supplementary Table 2 for the number of generations) under a 12h:12h light-dark cycle with 100 µmol photons m -2 s -1 provided by fluorescent tubes (Mazda 18WJr/865).

Growth response curves.

Cell concentration was determined on fresh samples using flow cytometry according to [START_REF] Marie | Enumeration of Marine Viruses in Culture and Natural Samples by Flow Cytometry[END_REF]. The maximum cell growth rate (µ max ) was calculated as the slope of the linear regression relating cell concentration logarithm vs. time observed during the exponential phase of growth. The Cardinal Temperature Model with Inflection (BR model) from [START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF] was used to estimate the optimal temperature of growth (T opt ) at which the growth rate is optimal (µ opt ), and the minimal and maximal temperatures of growth (T min and T max ) at which µ = 0. The growth µ(T) at temperature T is described as follows:

µ(T) =      0 for T < T min µ opt .φ(T) for T min < T < T max 0 for T > T max (3)
where

φ(T) = (T -T max )(T -T min ) 2 (T opt -T min )[(T opt -T min )(T -T opt ) -(T opt -T max )(T opt + T min -2T)
Selection of the thermal growth response model. Number of models exist that represent the response of phytoplankton strains to temperature; we selected the one we believe to be the most relevant for the purpose of the present study. We first short-listed the most appropriate models after the two recent reviews of [START_REF] Grimaud | Modeling the temperature effect on the specific growth rate of phytoplankton: a review[END_REF] and [START_REF] Low-Decarie | Predictions of response to temperature are contingent on model choice and data quality[END_REF]. Grimaud and colleagues [START_REF] Grimaud | Modeling the temperature effect on the specific growth rate of phytoplankton: a review[END_REF] discussed the strengths and limitations of several thermal response models in regard to four criteria: the fit quality, the easiness of calibration, the biological interpretation of parameters, and the applicability to phytoplankton growth. They convincingly argued that the BR (Eq. 3, [START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF]) and Eppley-Norberg (Eq. 4, [START_REF] Norberg | Biodiversity and ecosystem functioning: a complex adaptive systems approach[END_REF]) models presented the overall best performances. Following the analysis from [START_REF] Low-Decarie | Predictions of response to temperature are contingent on model choice and data quality[END_REF], we also considered the Boatman model (Eq. 5, [START_REF] Boatman | A key marine diazotroph in a changing ocean: The interacting effects of temperature, CO2 and light on the growth of Trichodesmium erythraeum IMS101[END_REF]) and we calibrated all three models to our growth measurements (Supplementary Figure 4 and5).

µ(T) = 1 - T -T opt w 2
ae bT where w = abs(T max -T min ) (4)

µ(T) = µ max sin π T -T min T max -T min a b (5)
We then computed an Akaike Information Criterion (AIC) and a Bayesian Information Criterion (BIC) for each model (Supplementary Table 5) according to the following equations:

AIC = 2k -2ln(MSE) (6) BIC = -2ln(AIC) + kln(n) ( 7
)
Where k is the number of model parameters to be estimated, MSE the Mean Square Error between measured and predicted growth rates and n the number of data points. These two criteria provide with an estimation of the relative quality of the models tested. Being an increasing function of MSE and k, the BIC is a selection criterion between models. The BR model yielded the smallest criteria and, in this regard, represented the best model tested to represent the growth response to temperature in Micromonas, in agreement with the findings of [START_REF] Grimaud | Modeling the temperature effect on the specific growth rate of phytoplankton: a review[END_REF] for other phytoplankton species.

Phylogenetic tree reconstruction and evolutionary placements.

Sequence alignment. 18S amplicon sequences from Micromonas RCC strains were aligned to a reference Mamiellophyceae sequence alignment. This reference alignment spans the rDNA operon and was originally used to describe the phylogenetic relationships amongst Mamiellophyceae genera (Marin and Melkonian, 2010). The reference alignment was trimmed to represent only the 18S rDNA region; long Micromonas RCC 18S amplicons (> 1000 nt; n = 35) were added to this alignment using MAFFT v7 [START_REF] Katoh | MAFFT multiple sequence alignment software version 7: improvements in performance and usability[END_REF]. The resulting alignment was then edited using the mask from the original alignment annotation [START_REF] Marin | Molecular phylogeny and classification of the Mamiellophyceae class. nov. (Chlorophyta) based on sequence comparisons of the nuclear-and plastid-encoded rRNA operons[END_REF] and was composed of a total of 2158 sites.

Phylogenetic tree reconstruction. The edited alignment was used for maximumlikelihood (ML) tree reconstructions. The best ML tree was identified from 100 independent tree reconstructions. All ML reconstructions were run using RAxML v8 [START_REF] Stamatakis | RaxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF] with the HKY85+G+I model, which was determined as the best-fit model of nucleotide substitution with jModelTest v2 [START_REF] Darriba | jModelTest 2: more models, new heuristics and parallel computing[END_REF] and by both the Akaike and Bayesian information criteria. Node supports of the resulting phylogenetic tree were determined using 1000 non-parametric bootstrap replicates. Bayesian inferences were conducted using BEAST v2 [START_REF] Bouckaert | BEAST 2: A Software Platform for Bayesian Evolutionary Analysis[END_REF] using the HKY85+I+G with a log-normal, relaxed molecular clock and default priors. A total of 4 MCMC chains of 10 6 generations were conducted, and a 25% 'burnin' value was applied on the resulting tree set. The iTol web-server [START_REF] Letunic | Interactive Tree Of Life v2: online annotation and display of phylogenetic trees made easy[END_REF] was used to generate vector scalable graphic rendering. Thermal environment analysis. To identify possible correlation of isolated strains to temperature, a non-metric dimensional scaling (NMDS) was realized on a Euclidean distance matrix computed on the thermal environment dataset (T - S , T S ,

T + S , T + S -T - S )
using the R package vegan [START_REF] Oksanen | vegan: Community Ecology Package. R package version 2.4-0[END_REF]. The stress value is the measure of how well the NMDS configuration represents the dissimilarities and is referred as the Kruskal stress [START_REF] Kruskal | Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis[END_REF].

Relation between strains and environmental temperatures. Relationships between environmental temperatures (T - S , T S , T + S , T + S -T - S ), latitude of the isolation site (Lat) and the species cardinal parameters (T min ,T opt ,T max and µ max ) were calculated for the eleven experimental strains that were grown in the laboratory.

We tested simple and multiple linear regression models and chose the best relationship according to a high R 2 adjusted and p-value < 0.05. Best relationships were obtained with T S and were used to determine cardinal parameters of all other 46 collection strains that were not experimentally tested but referenced in the RCC.

Thermotypes construction.

For each thermotype, we computed 100,000 growth vs. temperature curves through a Monte Carlo procedure with the BR model [START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF] and cardinal parameters of the i-th thermotype randomly taken from the parameter distributions (assuming a gaussian repartition of the parameters in the interval [p * -2σ, p * + 2σ] where p * are the parameters value). In order to ensure a biological coherence in the random samples of the cardinal parameters, the µ opt parameter is generated slightly differently. An Eppley model is used to link µ opt and T opt [START_REF] Eppley | Temperature and phytoplankton growth in the sea[END_REF]:

µ opt = a.e b.T opt (8)
Where parameters a and b are obtained from the best fit with all the strains of the thermotype (Supplementary Table 7). The values of µ opt for a random strain are then directly deduced from random values of T opt using this model.

Finally, we used the BR model to get the average thermal response and its standard deviation for each thermotype.

The optimal growth response envelope [START_REF] Grimaud | Modelling the temperature effect on phytoplankton: from acclimation to adaptation[END_REF] for the whole Micromonas genus was calculated with a BR curve calibrated on a data set consisting in 57 couples (T opt ,µ opt ) from the eleven experimental strains and the 46 collection strains.

Moreover, the decreasing part of the curve was constrained with 8 couples (T,µ(T))

simulated from the M. commoda Warm thermotype model for temperatures equally distributed in the (T opt ,T max ) interval for this thermotype. The increasing part of the curve was also constrained with eleven couples (T,µ(T)) simulated from the M. polaris model for temperatures equally distributed in the (T min ,T opt ) interval for this species.

Tara Oceans

Tara Oceans V9 dataset analysis. Molecular and contextual data from the Tara Oceans project were retrieved from PANGAEA [START_REF] Pesant | Open science resources for the discovery and analysis of Tara Oceans data[END_REF]. The Tara Oceans V9-18S

dataset [START_REF] Vargas | Eukaryotic plankton diversity in the sunlit ocean[END_REF] is available both at the barcode level (non-redundant sequences) and clustered at the Swarm/operational taxonomic unit level [START_REF] Mahé | Swarm: robust and fast clustering method for amplicon-based studies[END_REF]. Thermotypes inside the Tara Oceans V9 dataset. To explore the impact of temperature on species occurrence, we computed an NMDS on a Bray-Curtis distance matrix calculated from a community matrix of Micromonas species abundance per station (expressed in percentage of barcodes) with the R package "vegan" [START_REF] Oksanen | vegan: Community Ecology Package. R package version 2.4-0[END_REF].

Submitted to the ISME Journal

Results display a cloud of sampling stations from the different oceanic basins, discriminating surface and deep chlorophyll maximum (DCM); the closer proximity of stations, in terms of Bray-Curtis distances, expresses their similarities in their 18S diversity. We then fitted environmental variables (nutrients, temperature and mixed layer depth) and total chlorophyll a abundance on the ordination space with the vegan function env f it in vegan package [START_REF] Oksanen | vegan: Community Ecology Package. R package version 2.4-0[END_REF] with p-value based on 999 permutations was used to assess the significance of the fit.

The Micromonas distribution for each thermotype was computed against yearly SST (from NOAA) for each Tara Oceans station. We then computed Loess regressions with polynomial fitting to illustrate the temperature patterns with the R package "ggplot2" [START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF].

Global temperature response and diversity index

Global SST dataset. We used global SST data from the Copernicus Marine Environment Monitoring Service (product: GLOBAL_REP_PHYS_001_013) to calculate monthly averages SST in the period 1993 -2012 at the global scale.

Species distribution as a function of temperature. Cardinal parameters (T min ,

T opt , T max ) and optimum growth rate µ opt for each thermotype i were used to calculate the growth rate µ i (T) for each temperature T using the BR model [START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF].

Then, normalized distribution D i (T) of each thermotype was calculated following the equation:

D i (T) = µ i (T) ∑ n i=1 µ opt,i
for each temperature T of the global ocean surface.

Remark that this normalisation removes the effect of other factors which also influence net growth at the same location (nutrients, light, predations, etc.).

Diversity index.

To get a diversity index, we computed 10,000 thermal distribution via a Monte Carlo procedure for each species (Supplementary Figure 13).

We then computed an averaged and standard deviation of a Shannon-like based interspecific diversity index within the Micromonas genus according to Eq. 1 (Supplementary Figure 14) and compared it with a Shannon diversity index based on Tara Oceans V9 dataset thermotypes relative abundance:

H TARA (s) = n ∑ i=1 E(s, i)ln(E(s, i)) (9)
Where E(s, i) is the number of barcodes for the Micromonas thermotype i at the station s. The Tara Oceans dataset was used along the transect from station 4 to 125 [START_REF] Vargas | Eukaryotic plankton diversity in the sunlit ocean[END_REF]. The spatial distance between stations was calculated as a distance as the crow flies. In addition, we compare the Shannon-like base interspecific diversity Submitted to the ISME Journal index (Eq. 1) calculated for Micromonas (H M ) to the diversity index calculated by Thomas and colleagues for the phytoplankton (H P ) with a linear regression model (R 2 adj = 0.95 and p-value < 0.5):

H P = 83.21H M + 65.05 (10)
Then, we used Eq. 10 to quantify the diversity in the same index as the Thomas et al. study [3] (Supplementary Figure 15).

Cardinal parameters adaptation model

Cardinal Parameters Evolution. We studied the evolution of diversity in a warmer ocean with a dynamical model of the thermal growth response over the period 2001 to 2100. Projections of future, global temperature regimes were obtained from the NOAA GFDL CM2.1 [START_REF] Griffies | Formulation of an ocean model for global climate simulations[END_REF][START_REF] Delworth | GFDL's CM2 global coupled climate models. Part I: Formulation and simulation characteristics[END_REF] driven with the SRES A2 emissions scenario [START_REF] Nakicenovic | Special report on emissions scenarios: A special report of Working Group III of the Intergovernmental Panel on Climate Change[END_REF]. This dataset spans from 2001 to 2100 and was also used by Thomas and colleagues [3].

First, we computed the evolution of cardinal parameters T c,i (T min ,T opt and T max )

for each thermotype i depending on the temperature T(t, l, L) with t the year, l the latitude and L the longitude. The evolution of cardinal parameters follows Eq.

2, which is parameterized by the number of generations Na required to adapt to a different temperature (Supplementary Figure 16):

dT opt,i dt = N i (T(t)) Na (T * opt,i (T(t)) -T opt,i (t)) (11)
dT max,i dt = N i (T(t)) Na (T * max,i (T(t)) -T max,i (t)) (12)
Where T * opt,i and T * max,i are computed from the derivative of the relationships in Table 1 depending on the local temperature T(t, l, L):

dT * opt,i dt = 0.84 dT(t, l, L) dt ( 13
)
dT * max,i dt = 0.77 dT(t, l, L) dt (14)
The evolutive minimal temperature of growth was computed contingent to the Submitted to the ISME Journal evolution hypothesis:

T * min,i =      T ini min,i + T * max,i -T ini max,i
Constant thermal niche dT min,i dt

= N i (T(t)) Na (T * min,i (T(t)) -T min,i (t)) Dynamical model (Eq. 2) (15)
Where T ini min,i and T ini max,i are the initial value of T min,i and T max,i respectively at time t = 2001 and T * min,i is computed from the derivative of the relationships in Table 1 depending on the local temperature T(t, l, L):

dT * min,i dt = -0.92 dT(t, l, L) dt (16)
We constrained T * min,i and T * max,i by the envelope curve [START_REF] Grimaud | Modelling the temperature effect on phytoplankton: from acclimation to adaptation[END_REF] of the Micromonas genus (Figure 2b) that represents its evolution boundaries.

Second, we calculated µ opt,i at T opt,i with the BR model calibrated with the cardinal parameters of the envelope curve.

Third, we calculated the related growth rate µ i (T) of each thermotype i depending on its cardinal parameters T c,i at temperature T(t, l, L) following the BR model [START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF].

Third, we calculated the diversity for the present (2001 to 2010 -H now ) and future (2091 to 2100 -H f uture ) periods following the Eq. 1 averaged on 10 years and expressed as the diversity index used by Thomas and colleagues [3] with the Eq. 10.

Diversity erosion. We performed this cardinal parameter evolution framework for different values of Na, from fast (Na < 100 generations) as highlighted by [START_REF] Padfield | Rapid evolution of metabolic traits explains thermal adaptation in phytoplankton[END_REF][START_REF] Bonnefond | Continuous selection pressure to improve temperature acclimation of Tisochrysis lutea[END_REF] to slow (Na = 10 9 ) adaptation kinetics. This slow time scale corresponds to two to six months in the lab, which means a time scale in the range of years in the natural environment (assuming µ = 0.2 day -1 as a typical growth rate in the sea). For long-term evolution, we refer to a time scale slower than climate change. We call slow evolution an evolution with a typical adaptation kinetics with a millennium, which means Na = 10 6 generations for an average growth rate of 0.2 day -1 . We then calculated a diversity erosion index representing the loss of diversity along the latitude gradient with the equation:

H erosion (l) = h L L max L max ∑ L=0 (H now (l, L) -H f uture (l, L)) (17)
With L the longitude and l the latitude, L max the maximal longitude of the dataset (n = 359.7) and h the longitude resolution (h L = 0.1).

Submitted to the ISME Journal

The averaged latitudinal erosion (H erosion ) per latitude was calculated as follows:

H erosion = h l n l max ∑ l=l min (H erosion (l)) (18)
With l the latitude, l min and l max the minimum and maximum latitude of the dataset (l min = -82 and l max = 90), h l the latitude resolution (h l = 0.1) and n the H erosion vector's length. A negative erosion signifies a diversity gain.

The tipping point (p) of the H erosion vs. Na curve was calculated as the inflection point following the equation: growth response to temperature for each phylogenetic group computed from 100,000 possible response curves simulated within the ranges observed in each phylogenetic group. The black line represents the overall, optimal growth response envelope [START_REF] Grimaud | Modelling the temperature effect on phytoplankton: from acclimation to adaptation[END_REF] of Micromonas computed as µ opt vs. T opt , where µ opt and T opt are given by the average response of each thermotype. The grey shaded area is the standard deviation around µ opt .

p = max dH erosion dNa (19)

Choice of the light intensity during the experiments

The experiments were conducted at 100 µmol photons m 2 s 1 to find an optimal trade-off between non-photolimiting and non-photoinhibiting conditions. As supported by the work of [1], temperature growth response of phytoplankton is weakly coupled with light intensity for moderate light, whereas it could become more dependent to light at higher light intensities for which photoinhibition occurs. [2] showed that marked light-limitation can reduce the optimal growth temperature of phytoplankton by about 5 C. [3] studied the light growth response of Micromonas commoda and observed photoinhibition for light intensities higher than 300 µmol photons m 2 s 1 . In order to develop a model that accounts for the response to temperature, it was critical to experiment on the potential response, i.e. to assess the maximum growth capacity of strains at each temperature, not to introduce any bias (such as photoinhibition) in the experiments that would have led to an inaccurate estimation of the sole impact of temperature. The intensity of 100 µmol photons m 2 s 1 was therefore a reasonable trade-off. The BR model being tailored to account for light limitation on growth as well, it would then be possible to describe the coupled limitation of light and temperature, should it appear necessary. However, in the present study, the model proved to accurately compare to in situ data sets with the sole response to temperature, which indicated that additional model complexity through the inclusion of a light response was not necessary Supplementary Figure 4: Growth thermal response model fits for the 11 experimental strains. BR model (blue, solid line) with its 95% confidence interval (blue, shaded area), Eppley-Norberg model (red, solid line; [5]), Boatman model (green, solid line, [6]) and average experimental data (black circles) with their standard deviation (n at least = 3). Lat Supplementary Figure 6: Linear relationships between cardinal parameters and environmental parameters for the eleven strains tested experimentally. a) Relationships between T opt , T max and µ opt vs. the average surface temperature at the isolation site TS . b) Relationships between T min vs. averaged surface temperature at the isolation site TS , maximal surface temperature at the isolation site T + S and minimal surface temperature at the isolation site T S . Latitude at the isolation site is expressed with the color-bar. The star on top of the vertical axis represents a statistical significant relationship (p-value < 0.05).

warm cold M A M I E L L A C E A E C O C C A C E A E B A T H Y M A M I E L L A L

  mum T opt temperature (26.3±1.01 to 29.3±1.2 • C, respectively). Strains isolated in regions where the average SST fluctuates between 16.0 and 18.0 • C presented a lower optimal growth rate (0.9±0.03d -1 ) at T opt = 22.6 ± 3.08 • C) and maintained positive growth from 4.2±5.6 • C to 28.7±4.63 • C. In strains isolated at sites with Submitted to the ISME Journal an average temperature between 10.1 and 13.6 • , µ opt still reached 0.87±0.08d -1 at T opt =23.8±0.62 • C and cells demonstrated an ability to grow over a very wide temperature range (from -0.7±7.46 • C to 29.4±1.55 • C ). Last, Arctic strains (RCC2306 and RCC2257) revealed both the narrowest growth temperature range (-7.0±0 • to 15.1±0 • C) and lowest growth rates (0.45±0.03 d -1 ) at 7.5±0 • C.

  and 21), faster adaptation reduces the average diversity erosion from 4.5 species lost per latitude degree (slow adaptation) to one species lost or even 2 species gained per latitude (fast adaptation, Figure5b and d). Thermal adaptation performed within 200-300 generations might be sufficient to mitigate the impacts of climate change on phytoplankton diversity. In contrast, an adaptation scale beyond 10 4 generations will not counteract the deep impacts of climate change on phytoplankton diversity.The adaptation time scale of the thermal tolerance of different phytoplankton taxa has been closely related to their respective thermal environments (measured with T opt or the Net Primary Production)[START_REF] Huertas | Warming will affect phytoplankton differently: evidence through a mechanistic approach[END_REF][START_REF] Padfield | Rapid evolution of metabolic traits explains thermal adaptation in phytoplankton[END_REF][START_REF] Bonnefond | Continuous selection pressure to improve temperature acclimation of Tisochrysis lutea[END_REF][START_REF] Schaum | Adaptation of phytoplankton to a decade of experimental warming linked to increased photosynthesis[END_REF]. Phytoplankton taxa that ought to efficiently adapt to temperature are encountered in highly variable thermal environments[START_REF] Huertas | Warming will affect phytoplankton differently: evidence through a mechanistic approach[END_REF], typically found at latitudes beyond 40 • , where we found positive change in future diversity. These regions are also the main areas of CO 2 mitigation and carbon export in the ocean[2,[START_REF] Lutz | Seasonal rhythms of net primary production and particulate organic carbon flux to depth describe the efficiency of biological pump in the global ocean[END_REF]. The deeper alteration of phytoplankton diversity in the tropics might prove less critical for the efficiency of the biological pump at the global scale. Future research should be addressed to understand the impact of microbial diversity on carbon export[START_REF] Guidi | Plankton networks driving carbon export in the oligotrophic ocean[END_REF].ConclusionThis study describes niche partitioning in the marine pico-phytoplankton Micromonas. We showed that this genus evolved into different thermotypes that discriminate according to their sensitivity to temperature. Our model predictions were validated by in situ data from the Tara Oceans scientific expedition and suggest that temperature is a robust descriptor of Micromonas distribution at Submitted to the ISME Journal mesoscale and on a yearly basis. The diversity within this genus is highly correlated to the diversity pattern of the whole phytoplankton community. It is crucial to dedicate specific efforts to monitor the evolution of this sentinel genus in order to keep a real-time high fidelity picture of the phytoplankton diversity across the oceans. It is likely that Micromonas genus comprises even more thermotypes. More refined laboratory assessments including more thermotypes, should they exist, would enhance the representation of the global phytoplankton distribution. In particular, new experiments with smaller temperature increments and including more points at low and high temperatures would provide with a much higher resolution in the predicted capabilities and better assessment of T min and T max .Although decisive, the ability of phytoplankton to adapt in a warming ocean is the yet uncertain parameter. Adaptation is directly or indirectly affected by a variety of factors such as local nutrient availability, predation, virus lysis, mixing regime, etc. All of them are affected by the local physical dynamics and will also be impacted by global warming. More research is thus required to understand the adaptation mechanisms of this sentinel organism, and especially the adaptive dynamics of the different thermotypes. Such an approach will progressively refine the picture of phytoplankton evolution in a changing ocean with the possibility to more rapidly detect tipping points.

  Evolutionary placements. RCC 18S amplicon sequences shorter than 1000 nt (n = 24) were placed onto the ML phylogeny using the Evolutionary Placement Algorithm (EPA) implemented in RAxML v8[START_REF] Berger | Aligning short reads to reference alignments and trees[END_REF]. Short RCC sequences were aligned with MAFFT v7 against the previously generated updated reference Mamiellophyceae 18S alignment (i.e., composed of reference Mamiellophyceae and long RCC amplicon sequences). The aligned short sequences were then placed Submitted to the ISME Journal onto the reference phylogeny using RAxML in EPA mode with the HKY85+I+G model.Thermal niche partitioning analysis. Thermal environment dataset. Using SST from the National Oceanic and Atmospheric Administration's (NOAA), we built a dataset gathering the environmental temperatures at the isolation site of the eleven experimental and 46 Micromonas collection strains referenced in the RCC. At each strain's isolation site, we retrieved the yearly average SST (T S ), minimum SST (T - S ), maximal SST (T + S ) and thermal amplitude (T + S -T - S ) corresponding to a 10-year average (2005 to 2014).

  Micromonas-likeV9-18S barcode sequences were retrieved based on the original taxonomic classification from the Tara Oceans consortium, which was conducted with the Protist Ribosomal Reference database[START_REF] Guillou | The Protist Ribosomal Reference database (PR2): a catalog of unicellular eukaryote small sub-unit rRNA sequences with curated taxonomy[END_REF] for the protist barcode subset. The resulting 1084 non-redundant barcodes classified as Micromonas-like, and which represented a total of 95755 occurrences across the V9-18S Tara Oceans sampling (334 samples from 47 stations), were then re-classified using a phylogenetic placement procedure. The non-redundant Micromonas-like V9-18S barcodes were aligned against a reference Mamiellophyceae alignment using the same methodology than for the short 18S amplicon sequences from the Micromonas RCC strains, as aforementioned. The V9-18S barcode sequences were then placed onto the Mamiellophyceae and RCC reference tree using RAxML EPA with the HKY85+I+G model. Based on the placement of the Tara Oceans barcode onto the Micromonas reference subtree, the corresponding taxonomic information (thermotype level) was assigned to the environmental barcode.

Figure 2 :

 2 Figure 2: Original thermal environments and growth response to temperature for Micromonas species. (a) Two-dimensional ordination space derived from a Non-Metric MultiDimensional Scaling (NMDS) procedure displaying the thermal dissimilarities site (T - S , T S , T +S and T + S -T - S ) between the original isolation sites of the eleven experimental and 46 Micromonas collection strains. The stress value (goodness-of-fit of the NMDS) is inferior to 0.05, indicating high dimensional relationships among samples. (b) Average growth response to temperature for each phylogenetic group computed from 100,000 possible response curves simulated within the ranges observed in each phylogenetic group. The black line represents the overall, optimal growth response envelope[START_REF] Grimaud | Modelling the temperature effect on phytoplankton: from acclimation to adaptation[END_REF] of Micromonas computed as µ opt vs. T opt , where µ opt and T opt are given by the average response of each thermotype. The grey shaded area is the standard deviation around

Figure 3 :

 3 Figure 3: Micromonas thermotypes relative abundance patterns as estimated from the 18S rRNA V9 region during the Tara Oceans cruise. (a) Map of the Tara Oceans transect (dashed black line)showing station for which 18S rRNA V9 region data were available from Vargas et al. (2015) [27]: Mediterranean Sea (Med S), Red Sea (Red S), Indian Ocean (Ind O), South Pacific Ocean (S Pac O), Southern Ocean (S O) and South Atlantic Ocean (S Atl O). (b) Two-dimensional ordination space derived from an NMDS analysis displaying Bray-Curtis distance between the Micromonas species assemblages of the Tara Oceans stations, fitted by significant environmental variable (p-value < 0.05). The stress value (goodness-of-fit of the NMDS) is 0.15, indicating fair dimensional relationships among samples. (c) Relative abundance of the 6 thermotypes per station, plotted according to yearly SST at station coordinates: data (circles) and polynomial regression (solid line) fitted with the 95% confidence interval (shaded area). Number of observations for the 6 thermotypes are represented in histograms, plotted according to yearly SST at station coordinates.

Figure 4 :

 4 Figure 4: Estimated and predicted interspecific diversity within the Micromonas genus in the global ocean. (a) Estimated and predicted interspecific diversity within the Micromonas genus along the Tara Oceans transect as estimated from the Micromonas OTUs read abundances (blue circles) and as predicted from our diversity model (red circles), fitted by a polynomial regression with a 95% confidence interval. (b) Thermotypes proportions (%) from Tara Oceans dataset for different oceanic regions: Mediterranean Sea (Med S), Red Sea (Red S), Indian Ocean (Ind O), South Pacific Ocean (S Pac O), Southern Ocean (S O) and South Atlantic Ocean (S Atl O). (c) Predicted Shannon diversity index (H) calculated with the equation 1 using annual averages SST (Copernicus Marine Environment Monitoring Service, 1993 to 2012 satellite data). (d) Comparison of the latitudinal average diversity for all phytoplankton (from Thomas et al. 2012. black line) with that estimated by our Micromonas model. Shaded area represents the standard deviation from the mean along latitudes.

Figure 5 : 1 ) 3 )CP 4 ) 7 ) 6 )

 513476 Figure 5: Micromonas diversity changes in a warming ocean for two evolution hypotheses: (a-b) Specialist-generalist with constant thermal niche and (c-d) Specialist-generalist with dynamical thermal niche. (a-b) Latitudinal averaged diversity erosion calculated as the difference between diversity in present period (2001 to 2010) and future (2091 to 2100). Black line represents the diversity erosion from Thomas et al. 2012, red and blue line are the diversity erosion for the fast adaptation scenario (Na = 100) and slow adaptation scenario (Na = 10 6 ) respectively. Filled area represent the standard deviation to the mean along latitude. (c-d) Averaged diversity erosion per latitude calculated for different adaptation kinetic (from Na = 1 to Na = 10 6 generations): model results (black circles) and polynomial regression (blue line) fitted. The Tipping point is calculated as the inflexion point for the derivative of blue curve. The 20% loss point is calculated as 20% evolution from the lowest erosion scenario (Na = 1).

Supplementary Figure 16 :

 16 Example of evolution of cardinal parameters with the dynamical model (eq. 2 in the main manuscript). Here we presented the generic cardinal parameter T C (T min , T opt , T max ) between time t (T t C ) and time t + 1 (T t+1 C ) for two evolution time scales: fast (Na = 100 generations) and slow (Na = 4000 generation).

  • C and 35 • C (41.52 ± 30.61 generations on average, Supplementary Table

  • C), M. polaris was the most stenotherm species. M. commoda cold and M. bravo cold showed very similar responses at colder temperatures but discriminated in regard to the optimum growth rate and maximum temperature. Their thermal niche of 25.42 ± 3.75 • C and 27.10 ± 0.91 • C, respectively, was representative of cold-temperate environments. Contrary to the cold species, and although they both live in warmer biotopes, the warm thermotype of species M. commoda and M. bravo showed very distinct thermal niches (34.00 ± 1.19 • C and 26.02 ± 5.11, respectively). Last, M. pusilla was found in both cold-and warm-temperate areas and showed an intermediate thermal response compared to the other Micromonas species, with a thermal niche of 28.85 ± 5.32 • C. With the most variable response to temperature, M. pusilla did not seem to speciate into different thermotypes; yet it clearly differentiated from other groups and would be the most eurytherm.

Acknowledgements

We thank members of the RCC for their help and work to maintain many phytoplankton strains in culture. We are also grateful to Mridul Thomas for providing his phytoplankton diversity data [3] and to Keith Paarporn for proofreading the text. Finally, we thank the reviewers and editor for their constructive comments on a previous version of this manuscript. This research was funded by the Inria Project Lab Algae in silico and by the ANR funding agency REVIREC (grant no.

12-BSV7-0006-01).

Conflit of interest statement

Authors declare no conflit of interest for this study. Micromonas experimental strains used in this study, plotted against yearly average SST for the year 2014 (from NOAA). (b) Growth rate vs. temperature curves for strains isolated in environments with different annual average temperature (T S ), fitted by the BR model [START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF]. Error bars are standard deviations (n 3).

Author contributions

Table 1: Linear relationship between cardinal parameters and environmental parameters (average temperature at the surface of isolation site, T S , and the latitude, Lat) for the eleven Micromonas experimental strains tested in this study.

1

Picoeukaryotes of the Micromonas genus: sentinels of a warming ocean Supplementary Figure 1: Summary of the scientific approach Supplementary Table 1: Information regarding the Micromonas strains used in the study. The thermal environment at the isolation sites is expressed in C: yearly averaged temperature ( TS ), minimal temperature (T S ), and maximal temperature (T + S ). The latitude of each isolation site (Lat) is expressed in degree and the growth temperature of cultures (T RCC ) is in C.

RCC # Species

Thermotype Supplementary Table 8: Cardinal parameters (in C), optimal growth rate (in day 1 ) and thermal niche width (in C) of the six thermotypes with their associated standard deviation.

Thermotypes