Paul Frémont
email: pfremont@genoscope.cns.fr

Marion Gehlen
email: marion.gehlen@lsce.ipsl.fr

Mathieu Vrac

Jade Leconte

Tom O Delmont

Patrick Wincker

Daniele Iudicone

Olivier Jaillon
email: ojaillon@genoscope.cns.fr

Restructuring of genomic provinces of surface ocean plankton under climate change

The impact of climate change on diversity, functioning and biogeography of marine plankton is a major unresolved issue. Here, niche theory is applied to plankton metagenomes of 6 size fractions, from viruses to meso-zooplankton, sampled during the Tara Oceans expedition. Niches are used to derive plankton size-dependent structuring of the oceans south of 60°N in climato-genomic provinces characterized by signature genomes. By 2090, assuming the RCP8.5 high warming scenario, provinces would be reorganized over half of the considered oceans and quasi-systematically displaced poleward. Particularly, tropical provinces would expand at the expense of temperate ones.

Compositional shifts among planktonic grazers and nitrogen-fixing bacteria suggest impacts on the nitrogen and carbon cycles. Sea surface temperature is identified as the main driver of the changes (~51%) followed by phosphate (11%) and salinity (10%). These results demonstrate the potential of integration of genomics with physico-chemical data for higher scale modeling and understanding of ocean ecosystem dynamics. Planktonic communities are composed of complex and heterogeneous assemblages of small animals, single-celled eukaryotes (protists), bacteria, archaea and viruses -that drift with currents. They contribute to the regulation of the Earth system through primary production via photosynthesis 1 , carbon export to the deep oceans 2,3 and form the base of the food webs that sustain the whole trophic chain in the oceans and beyond 4 .

The composition of communities varies over time at a given site with daily 5 to seasonal fluctuations 6 following environmental variability 7 . Overlying these relatively short scale .

spatio-temporal variations, a more macroscale partitioning of the ocean has been revealed by different combinations of biological and physico-chemical data [8][9][10] , and recently documented at the resolution of community genomics 11 . The basin scale biogeographical structure has been proposed to result from a combination of multiple bio-physico-chemical processes named the seascape 7 . These processes include both abiotic and biotic interactions 12 , neutral genetic drift 13 , natural selection 14,[START_REF] Delmont | Single-amino acid variants reveal evolutionary processes that shape the biogeography of a global SAR11 subclade[END_REF] , temperature variations, nutrient supply but also advection and mixing along currents 11,13 .

Today, knowledge of global scale plankton biogeography at the DNA level is in its infancy. We lack understanding and theoretical explanations for the emergence and maintenance of biogeographical patterns at genomic resolution. Omics data (i.e. the DNA/RNA sequences representative of the variety of coding and non-coding sequences of organisms) provide the appropriate resolution to track and record global biogeographical features 11 , modulation of the repertoire of expressed genes in a community in response to environmental conditions 2,16,17 and eco-evolutionary processes [13][14][START_REF] Delmont | Single-amino acid variants reveal evolutionary processes that shape the biogeography of a global SAR11 subclade[END_REF] . Moreover, metagenomic sequencing can be consistently analyzed across plankton organisms as recently demonstrated by global expeditions [18][19][20][21] . The strong links between plankton and environmental conditions suggest potentially major consequences of climate change on community composition and biogeography 22,23 .

Time series observations have highlighted recent changes in the planktonic ecosystem attributed to anthropogenic pressures, such as changes in community composition [START_REF] Michelangeli | Probabilistic downscaling approaches: Application to wind cumulative distribution functions[END_REF][START_REF]The probable error of a mean[END_REF] or poleward shifts of some species [START_REF] Watson | A database of global marine commercial, small-scale, illegal and unreported fisheries catch 1950-2014[END_REF] . These changes are expected to intensify with ongoing climate warming [START_REF]Maritime Boundaries Geodatabase: Maritime Boundaries and Exclusive Economic Zones (200NM)[END_REF] and could lead to major reorganization of plankton community composition 22 , with a potential decline in diversity [START_REF] Delmont | Functional repertoire convergence of distantly related eukaryotic plankton lineages revealed by genome-resolved metagenomics[END_REF][START_REF] Biecek | DALEX: explainers for complex predictive models[END_REF][START_REF] Barton | Anthropogenic climate change drives shift and shuffle in North Atlantic phytoplankton communities[END_REF] . Another major consequence of global reorganization of the seascape on biological systems would be a decrease of primary production at mid-latitudes and an increase at higher latitudes [START_REF]Maritime Boundaries Geodatabase: Maritime Boundaries and Exclusive Economic Zones (200NM)[END_REF] .

Here we report the global structure of plankton biogeography south of 60°N based on metagenomic data using niche models and its putative modifications under climate change. First, we define environmental niches 31 , i.e. the envelope of environmental parameters suitable for an organism or a population, at the scale of genomic provinces across 6 organism size fractions representing major plankton groups from nano-(viruses) to meso-zooplankton (small metazoans). Then, we spatially extrapolate their niches into climato-genomic provinces to derive the structure of plankton biogeography for each size fraction individually and for all combined. Next, considering the same niches, we assess the spatial reorganization of these provinces under climate change at the end of the century, with a focus on associated compositional shifts among copepods (planktonic grazers important for the carbon cycle) and nitrogen-fixing bacteria (important for both nitrogen and carbon cycles). Finally, we quantify the relative importance of the environmental drivers explaining projected changes.

Niche models and signature genomes from genomic provinces

We define and validate environmental niches using 4 machine learning techniques for 27 previously defined genomic provinces 11 ; they correspond to 529 metagenomes for 6 size fractions (ranging from 0 to 2000 µm) sampled at 95 sites from all oceans except the Arctic (Supplementary information 1, Supplementary Figs. 12). Predictor variables of the niches are sea surface temperature (SST), salinity, dissolved silica, nitrate, phosphate and iron, plus a seasonality index of nitrate.

The signal of ocean partitioning is likely due to abundant and compact genomes whose geographical distributions closely match provinces. Within a collection of 1778 bacterial, 110 archaeal and 713 eukaryotic environmental genomes 32,33 characterized from Tara Oceans samples without cultivation, we find a total of 324 signature genomes covering all but 4 provinces, and displaying taxonomic signal coherent with the size fractions (Fig. 1 for eukaryotes and Supplementary Fig. 3 for Bacteria and Archaea). Some of the signature genomes correspond to unexplored lineages with no cultured representatives, highlighting the knowledge gap for organisms that structure plankton biogeography and the strength of a rationale devoid of any a priori on reference genomes or species.

Structure of present day biogeography of plankton

To extrapolate the niches to a global ocean biogeography for each size fraction, we define the most probable provinces, named hereafter as dominant and assigned to a climatic annotation (Supplementary Table 1), on each 1°x1° resolution grid point using 2006-13 WOA13 climatology 34 (Supplementary Fig. 4 and Fig. 2).

In agreement with previous observations 11 , provinces of large size fractions (>20 µm) are wider and partially decoupled from those of smaller size fractions, probably due to differential responses to oceanic circulation and environmental variations, different life cycle constraints, lifestyles 7,11 and trophic network positions 35 . Biogeographies of small metazoans that enrich the largest size fractions (180-2000 and 20-180 µm) are broadly aligned with latitudinal bands (tropico-equatorial, temperate and (sub)-polar) dominated by a single province (Fig. 2ab). A more complex oceanic structuring emerges for the smaller size fractions (<20 µm) (Fig. 2c-f) with several provinces per large geographical region. For size fraction 0.8-5 µm enriched in small protists (Fig. 2d), distinct provinces are identified for tropical oligotrophic gyres and for the nutrient-rich equatorial upwelling region. A complex pattern of provinces, mostly latitudinal, is found for the bacteria (Fig. 2e, 0.22-3 µm) and the virus enriched size classes (Fig. 2f, 0-0.2 µm) although less clearly linked to large-scale oceanographic regions. A single province extending from temperate to polar regions emerges for size fraction 5-20 µm enriched in protists (Fig. 2c), for which fewer samples were available (Supplementary Fig. 2b-c), which probably biases this result. A consensus map combining all size fractions, built using the PHATE algorithm 36 , summarizes the main characteristics of the biogeographies described above (Supplementary information 2 and Supplementary Fig. 5).

Finally, we compare genomic biogeographies and existing ocean partitionings [8][9][10] .

Though each of them is unique (Supplementary Fig. 8910), common borders highlight a global latitudinal partitioning independent of the type of data (Supplementary information 3).

Future changes in plankton biogeography structure

We assess the impacts of climate change on plankton biogeography at the end of the century following the Representative Concentration Pathway 8.5 (RCP8.5) 37 greenhouse gas concentration trajectory. To consistently compare projections of present and future biogeographies, we use bias-adjusted mean of 6 Earth System Model (ESM) climatologies (Supplementary Table 2, Supplementary Fig. 11). The highest warming (7.2°C) is located off the east coast of Canada in the North Atlantic while complex patterns of salinity and nutrient variations are projected in all oceans (Supplementary Fig. 12). Following this trajectory, future temperature at most sampling sites will be higher than the mean and maximum contemporary temperature within their current province (Supplementary Fig. 13).

Our projections indicate multiple large-scale changes in biogeographical structure including expansions, shrinkages and shifts in plankton organism size-dependent provinces (Fig. 2a-d, Supplementary Fig. 141516). A change in the dominant province in at least one size fraction would occur over 60.1% of the ocean surface, ranging from 12% (20-180 µm) to 31% (0.8-5 µm) (Fig. 2, Table 1).

Centroids of provinces with dominance areas larger than 10 6 km 2 within a basin would be moved at least 200 km away for 77% of them, 96% of which move poleward (Supplementary Figs. 15 and16). While a few longitudinal shifts larger than 1000 km are projected, the distribution of latitudinal shifts is largely concentrated around the mean (290 km) with no shifts superior to 1000 km (Supplementary Fig. 16b). These important longitudinal shifts corroborate existing projections 22,38,39 and differ from trivial poleward shifts due to temperature increase, reflecting more complex spatial rearrangements of the other environmental drivers (Supplementary Fig. 12). The average displacement speed of the provinces' centroids is 76 ± 79 km.dec -1 (latitudinally mean of 34 ± 82 km.dec -1 , longitudinally 59 ± 82 km.dec -1 ).

Projected shifts in phytoplankton enriched provinces corroborate previously published shifts of North Pacific phytoplankton biomes: provinces C4 and C9 are projected to shift respectively at speeds of 118 km.dec -1 and 195 km.dec -1 comparable to 100 km.dec -1 and 200 km.dec -1 for the subtropical and tropical biomes of Polovina et al. 39 . For all size fractions, climate change would lead to a poleward expansion of tropical and equatorial provinces at the expense of temperate provinces (Supplementary information 4, Supplementary Table 2, Supplementary Figs. 14 and17). This is illustrated by the temperate province F5 of size fraction 180-2000 (Supplementary Fig. 16), which is projected to shrink in the five major basins. In the North Atlantic, its centroid would move approximately 800 km to the northeast (Supplementary Fig. 16c). Similar trends are found comparing present day and end of the century consensus maps (Supplementary information 2, Supplementary Fig. 17).

We calculate a dissimilarity index (equation (3)) at each grid point between probabilities of future and present dominant provinces for all size fractions combined (Fig. 4a).

Large dissimilarities are obtained over northern (25° to 60°) and symmetrically southern (-25 to -60°) temperate regions (mean of 0.29 and 0.24 respectively) mostly reflecting the poleward retraction of temperate provinces (red arrows, Fig. 4a). In austral and equatorial regions, despite important environmental changes (Supplementary Fig. 12) and previously projected changes in diversity [START_REF] Delmont | Functional repertoire convergence of distantly related eukaryotic plankton lineages revealed by genome-resolved metagenomics[END_REF][START_REF] Biecek | DALEX: explainers for complex predictive models[END_REF][START_REF] Barton | Anthropogenic climate change drives shift and shuffle in North Atlantic phytoplankton communities[END_REF] and biomass 40 , the contemporary provinces remain the most probable at the end of the century (mean dissimilarities of 0.18 and 0.02 respectively).

To further study the future decoupling between provinces of different size fractions, we analyze the assemblages of dominant provinces of each size fraction. By using two differently stringent criteria, from 45.3 to 57.1% of ocean surface, mainly located in temperate regions, would be inhabited in 2090-99 by assemblages that exist elsewhere in 2006-15 (Fig. 4b versus Fig. 4c). Contemporary assemblages would disappear on 3.5 to 3.8% of the surface, and, conversely, novel assemblages, not encountered today, would cover 2.9 to 3.0% of the surface. While these changes are limited to a relatively modest area, they include important economic zones (Fig. 4b, Supplementary Fig. 19).

Over 41.8% to 51.8% of the surface of the main fisheries and 41.2% to 54.2% of Exclusive Economic Zones, future assemblages would differ from those present today (Supplementary Fig. 19).

Future changes in the distribution of grazers and nitrogen-fixing bacteria

In order to elucidate the potential biogeochemical impact of biogeographical restructuring, we focus on compositional changes among copepods and nitrogen-fixing bacteria (a.k.a, diazotrophs), two groups considered important for the carbon and nitrogen cycles 41,42 that are well represented among the Tara Oceans environmental genomes. Focusing on marine areas where dominant provinces are projected to be replaced, we compare the present and future distribution of environmental genomes corresponding to 198 copepods and 27 diazotrophs 32,33 .

Copepods are cosmopolitan small crustaceans. These abundant grazers contribute to the biological pump and their body size is considered to be a key trait for carbon export 41 . They feed on smaller plankton 41 and are prey to higher trophic levels. We Therefore, large copepods are preferentially from clade A in our data.

The relative abundances of unclassified clade B in size class 20-180 μm (Fig. 3e) and of mesozooplankton clade A in size class 180-2000 μm, (Supplementary Fig. 20a) are projected to increase in regions where the temperate province is replaced by the tropico-equatorial province. In areas where the polar province is replaced by the temperate province, a greater relative abundance in microzooplankton clade A (1% to 10%) (Fig. 3e) is projected. No significant compositional differences in the provinces of size fraction 5-20 μm are found (Supplementary Fig. 20b). These results highlight potential significant compositional shifts by the end of the century in the different clades and sizes of main grazers.

Diazotrophy, the biotic fixation of atmospheric nitrogen, is an important process for both nitrogen and carbon cycles. It supports biological productivity in the nitrogen-limited tropical oceans 42 . Marine diazotrophs include cyanobacteria (e.g. Trichodesmium) described as the principal nitrogen fixers 43 and various heterotrophic bacterial diazotrophs (HBDs) that lack cultured representatives or in situ imaging 33,44 . Forty-eight of the bacterial environmental genomes are diazotrophs (8 cyanobacteria and 40 heterotrophic bacterial diazotrophs (HBDs)), encapsulating 92% of the metagenomic signal for known nifH genes (a reference marker for nitrogen fixation 45 ) at the surface of the oceans. Twenty-seven are found in at least 5 samples and have been used for the analysis.

In size fraction 0.8-5 μm, we project significantly higher relative abundances in cyanobacteria in the tropico-equatorial regions of the Pacific ocean (C9 to C11 and C4, Fig. 3f), which might point towards an increase in nitrogen fixation in this region as previously suggested by other models 42 . Supporting this result, we find similar significant compositional changes towards an increase in cyanobacteria in size fraction 5-20 μm (Supplementary Fig. 21c). We also project significant compositional changes for some clades of HBDs (e.g. increase in gammaproteobacteria: C8 to C3, Fig. 3f) though no global trend can be identified here. In the other size classes (20-180, 180-2000 and 0.22-3 μm) compositional changes are not significant (Supplementary Fig. 22). To summarize, although we cannot estimate nitrogen fixation rates using genomic data alone, genomic measurements of nitrogen-fixing cyanobacteria are in agreement with an increase in nitrogen fixation in the tropics, echoing results from other models 42 .

Drivers of plankton biogeography reorganization

We quantify the relative importance of environmental predictors (sea surface temperature, salinity, dissolved silica, phosphate, nitrate, iron and seasonality of nitrate) into niche definition and in driving future changes of the structure of plankton biogeography (equation ( 4)). Among these environmental properties, temperature is the first influential parameter (for 19 niches out of 27) but only at 22.6% on average (Supplementary Fig. 22a).

The relative impact of each environmental parameter is calculated 22 for each site presenting a significant dissimilarity between 2006-15 and 2090-99 (Fig. 5a). Overall, SST would be responsible for the reorganization of the provinces at 50% followed by Phosphate (11%) and Salinity (10.3%) (Supplementary Fig. 23). Over the majority of the ocean, SST is the primary driver of the reorganization (Fig. 5a). In some regions, salinity (e.g. eastern North Atlantic) and Phosphate (e.g. equatorial region) dominate (Fig. 5a).

When excluding the effect of SST, salinity and phosphate become the primary drivers of the reorganization of the provinces (Fig. 5b). The impact of SST varies across size classes with a significantly higher contribution in large size classes (>20 µm) compared to the small ones (mean of ~73% versus ~49%, t-test p<0.05; Fig. 5c). Though the contribution of combined nutrients to niche definition is similar for small and large size classes (mean of ~56% versus ~61%, Supplementary Fig. 22, Supplementary Table 3), their future projected variations have a higher relative impact on the reorganization of biogeographies of small organisms (mean of ~39% versus ~20%, t-test p<0.05, Supplementary Fig. 22, Supplementary Table 3). For instance, in the tropical zone, the shrinkage of the equatorial province C9 (size fraction 0.8-5 µm, Fig. 2b,d, Supplementary Fig. 24e) is driven at 24% by reduction of dissolved phosphate concentrations and at 25% by SST increase. In contrast, SST drives at 56% the shrinkage of the temperate province F5 (size fraction 180-2000 µm, Supplementary Fig. 24d). Finally, non-poleward shifts are found only within small size fractions (<20 µm) (Supplementary Figs. 14, 15) highlighting differential responses to nutrients and SST changes between large and small organisms, the latter being enriched in phytoplankton that directly rely on nutrient supplies.

Discussion

We propose a novel partitioning of the ocean in plankton size dependent climatogenomic provinces, complementing previous efforts based on other bio-physicochemical data [9][10][11] . Though initially built at genomic scale, our biogeographies paradoxically reveal basin scale provinces that are larger than BGCP 10

and Fay and

McKingley biomes 9 . These provinces are probably relatively stable across seasons suggesting limited effects of seasonality on the position of frontiers of BGCPs provinces 10 . We propose that this apparent paradox emerges from the combination of the scale, nature and resolution of sampling. First, two proximal samples from the Tara Oceans expedition are separated by ~300 km on average sampled over three years. This relatively large spatio-temporal scale overlies shorter scale compositional variations previously observed 5,6 . Second, our estimates of plankton community dissimilarities are highly resolutive as they are computed at genomic scale with billions of small DNA fragments 11,18 thus smoothing out the more discrete species level signal. Together, from these combinations of processes and patterns occurring at multiple scales emerge basin scale provinces associated with coherent environmental niches and signature genomes.

These climato-genomic provinces are structured in broad latitudinal bands with smaller organisms (<20 μm) displaying more complex patterns and partially decoupled from larger organisms. This decoupling is the result of distinct statistical links between provinces based on organism size fractions and environmental parameters and could reflect their respective trophic modes 35 .

Complex changes of the parameters defining the niches are projected under climate change leading to the reorganization of size-dependent provinces. Assuming a constant relationship to environmental drivers that define the climato-genomic provinces, climate change is projected to restructure them over approximately 50% of surface oceans south of 60°N by the end of the century (Fig. 4). The largest reorganization is detected in subtropical and temperate regions in agreement with other studies [START_REF] Biecek | DALEX: explainers for complex predictive models[END_REF]39 and is accompanied by appearance and disappearance of size-fractionated provinces' assemblages. Out of contemporary range and novel environmental conditions are projected for tropico-equatorial and austral regions. While some studies extrapolate important diversity and biomass changes in these zones [START_REF] Delmont | Functional repertoire convergence of distantly related eukaryotic plankton lineages revealed by genome-resolved metagenomics[END_REF][START_REF] Biecek | DALEX: explainers for complex predictive models[END_REF][START_REF] Barton | Anthropogenic climate change drives shift and shuffle in North Atlantic phytoplankton communities[END_REF]40 , here we project shifts of their boundaries and maintenance of their climatic label. The present approach does not account for putative changes in community composition or the emergence of novel niches over these regions for which novel environmental selection pressure is expected.

Despite these limitations, genomic data allow us to project compositional shifts due to the reorganization of the provinces at a high phylogenetic and functional resolution.

Large (size class 180-2000 µm) and small (size class 20-180 µm) copepods are expected to be in greater relative abundance in subtropical and subpolar regions, respectively. Though the Tara Oceans environmental genomes collection might not reflect the full diversity in copepods, this projected restructuration of the communities in small and large copepods could modify carbon export as grazers' size is a key trait for this process 41 . Secondly, it might lead to novel prey-predator interactions e.g. in regions where new assemblages of communities are projected.

Important compositional shifts among most abundant marine nitrogen-fixing bacteria are also found. The relative share of nitrogen-fixing cyanobacteria is projected to increase in tropico-equatorial regions congruent with modeling studies 42 . Though their contribution to nitrogen fixation is not characterized, heterotrophic bacterial diazotrophs are abundant in genomic samples 33 and significant compositional shifts in certain clades are reported here. Ultimately, further studies associating omics and physico-chemical data should be informative for biogeochemical modeling of nitrate fixation rates and carbon export fluxes.

Overall, our projections for the end of the century do not take into account possible future changes of major bio-physico-chemical factors such as the dynamics of community mixing, trophic interactions through transport 46 , the dynamics of the genomes [13][14][START_REF] Delmont | Single-amino acid variants reveal evolutionary processes that shape the biogeography of a global SAR11 subclade[END_REF] (adaptation or acclimation) and biomass variations 40 . New sampling in current and future expeditions 47 , as well as ongoing technological improvements in biophysico-chemical characterization of seawater samples 32,47,48 , will soon refine functional 16,49 , environmental (micronutrients 50 ) and phylogenetic 32 characterization of plankton ecosystems for various biological entities (genotypes, species or communities) and spatio-temporal scales 47 . Ultimately, integrating this varied information will allow a better understanding of the conditions of emergence of ecological niches in the seascape and their response to a changing ocean. is the number of samples where the MAG is present and match a sample of the province. and are respectively the number of samples where the MAG is not present in a sample of the province and inversely. A MAG is considered to be signature of a province if the Jaccard index is superior to 0.5 with this province and inferior to 0.1 for all other provinces of the given size class (Fig. 1 and Supplementary Fig. 3).

World Ocean Atlas data

Physicochemical parameters proposed to have an impact on plankton genomic Nutrients, such as NO 3 and PO 4 , display a strong collinearity when averaged over the global ocean (correlation of 0.95 in WOA13) which could complicate disentangling their respective contributions to niche definition. However, observations and experimental data allow identification of limiting nutrients at regional scale characterized by specific plankton communities [START_REF] Moore | Processes and patterns of oceanic nutrient limitation[END_REF] . The projection of niches into future climate would yield spurious results when the present-day collinearity is not maintained [START_REF] Brun | The predictive skill of species distribution models for plankton in a changing climate[END_REF] but there is up to now no evidence for large scale changes in global nutrient stoichiometry [START_REF] Redfield | On the Proportions of Organic Derivatives in Sea Water and Their Relation to the Composition of Plankton[END_REF] .

Earth System Models and bias correction

Outputs from 6 Earth System Models (ESM) (Supplementary Gradient Boosting Machine (gbm) [START_REF] Ridgeway | gbm: Generalized Boosted Regression Models[END_REF] , Random Forest (rf) [START_REF] Breiman | Breiman and Cutler's random forests for classification and regression[END_REF] , fully connected Neural Networks (nn) [START_REF] Venables | Modern Applied Statistics with S Fourth edition by[END_REF] and Generalized Additive Models (gam) [START_REF] Wood | Stable and efficient multiple smoothing parameter estimation for generalized additive models[END_REF] . Hyper parameters of each technique (except gam) are optimized. These are (1) for gbm, the interaction depth (1, 3 and 5), learning rate (0.01, 0.001) and the minimum number of observations in a tree node (1 to 10); (2) for rf, the number of trees (100 to 900 with step 200 and 1000 to 9000 with step 2000) and the number of parameters used for each tree (1 to 7); (3) for nn, the number of layers of the network (1 to 10) and the decay (1.10 -4 to 9.10 -4 and 1.10 -5 to 9.10 -5 ). For gam the number of splines is set to 3, respectively 2 only when not enough points are available (for fraction 0-0. ("presence" or not of the province in the sample) and the predictors consisting of the environmental variables for each sample. A fixed probability threshold of 0.5 for presence/absence detection is used to calculate the AUC, i.e. samples with P, the probability of presence calculated by the model, superior to 0.5 are assigned to 1 and inversely samples with P inferior to 0.5 are assigned to 0. Fixing the probability threshold allows optimization of all models according to this threshold so that within a size fraction the dominant province has a reasonably high probability of presence (at least in regions with similar environmental parameters to the training dataset) and for the four types of models we use (gbm, nn, rf and gam). The best combination of hyper parameters is the one for which the mean AUC over the cross-validation is the highest.

A model is considered valid if at least 3 out of the 4 techniques have a mean AUC superior to 0.65, which is the case for 27 out of the 38 provinces (Supplementary Fig. 2a). A climatic annotation is given to the 27 validated niches (Supplementary Table 2).

Final models are trained on the full dataset and only the techniques that have a mean AUC higher than 0.65 are considered to make the projections. The vast majority (23) of the 27 validated niches is validated by all four models and 4 by only 3 models. Relative influences of each parameter in defining environmental niches are calculated using the feature_importance function from the DALEX R package [START_REF] Biecek | DALEX: explainers for complex predictive models[END_REF] for all four statistical methods (Supplementary Fig. 22a). To evaluate the consistency and coherence of environmental niche models, we first make global projections on the 2006-13 WOA2013 climatology.

Projections are consistent with sampling regions for provinces encompassing vast oceanic areas. For example, the genomic province sampled in temperate Atlantic regions of size fraction 180-2000 µm is projected to be present in the north and south temperate Atlantic but also other temperate regions (Supplementary Fig. 4). For model training and projections, physicochemical variables are scaled to have a mean of 0 and a variance of 1. For this scaling, the mean and standard deviation of each WOA13 variable (+ PISCES-v2 Fe) co-localized with Tara Oceans stations with a value available are used. This standardization procedure allows for better performance of nn models.

Finally, as statistical models often disagree on projection sets whereas they give similar predictions on the training set (Supplementary Fig. 6,7), we use the ensemble model approach for global-scale projections of provinces [START_REF] Jones | Multi-model ensemble projections of climate change effects on global marine biodiversity[END_REF] i.e. the mean projections of the validated machine learning techniques.

Combined size class provinces and ocean partitioning comparisons

To combine all size classes' provinces, we use the PHATE algorithm 36,64 from the R package phateR. This algorithm allows visualization of high dimensional data in the requested number of dimensions while best preserving the global data structure [START_REF] Vallejos | Exploring a world of a thousand dimensions[END_REF] . We choose to train PHATE separately on WOA13 projections and present day and end of century projections including presence probabilities of non dominant provinces. We use 3 dimensions and set hyper parameter k-nearest neighbors (knn) and decay respectively to 1000 for WOA13 and 2000 for model data as in this case there are twice as many points. The hyper parameter knn reflects the degree to which the mapping of PHATE from high to low dimensionality should respect the global features of the data.

We argue that 1000 and 2000 are good choices as it will be sufficient to have a highly connected graph, conserve global structure, allow visualization of structures of the size of the provinces (mean number of points in a province: 4867) and have a reasonable computational time. Decay is set to 20 in both cases. Then we cluster the resulting distance matrix using the k-medoïds algorithm [START_REF]Clustering by means of Medoids[END_REF] and the silhouette average width criteria [START_REF] Rousseeuw | Silhouettes: A graphical aid to the interpretation and validation of cluster analysis[END_REF] is used as an indicator of good fit. The silhouette criterion is maximal for 2, 3 and 4 clusters and 2 peaks are found at 7 and 14 clusters (the peak at 7 is slightly less high than the one at 14, data not shown). We choose to present the 4 and 7 cluster geographical patterns as they seem more relevant with respect to the resolutions of our environmental datasets (WOA13 and climate models). We compare the three polar clusters of the 7 cluster geographical patterns with Antarctic Circumpolar Currents fronts [START_REF] Orsi | On the meridional extent and fronts of the Antarctic Circumpolar Current[END_REF] by overlying them on the map (black lines Supplementary Fig. 5b).

To visualize the global biogeography structure, the resulting 3 vectors of PHATE are plotted using an RGB color code. Each coordinate of each vector is respectively assigned to a given degree of color component between 0 and 255 (8 bits red, green or blue) using the following formula (Supplementary Figs. 5,17):

(

) 2
is the ith component of the PHATE axes. Respectively, components 1, 2 and 3 are assigned to red, green and blue.

To compare the six size fraction provinces, the combined size class with existing biogeochemical partitions of the oceans 9,10 and with each other, we use the adjusted rand index 68 (Supplementary Fig. 8910) and overlay their masks above our partitions. In this case, presence probabilities of dominant provinces are not used anymore. Instead, each ocean grid point is assigned to the dominant provinces or the phate clusters.

Centroids and migration shifts

The centroid of each province is defined as the average latitude and longitude for which the probability of presence is superior to 0.5 and weighted by both the probability of presence at each grid point and the grid cell area. It is calculated for both present day conditions and end of the century conditions. The migration shift is calculated as the distance between the present day and the end of the century centroids considering the earth as a perfect sphere of radius 6371 km. For consistency (i.e. avoid long distance aberrant shifts), it is only calculated for provinces with an area of dominance larger than 10 6 km 2 in the given basin.

Bray-Curtis dissimilarity index

Climate change impact on global projections is calculated at each grid point as the Bray-Curtis dissimilarity index [START_REF] Somerfield | Identification of the Bray-Curtis similarity index: Comment on Yoshioka[END_REF][START_REF] Bloom | Similarity Indices in Community Studies: Potential Pitfalls[END_REF] defined as follows:

(

Where ( and ) are respectively the probability of presence of the province n in present day and at the end of the century. Only the probabilities of dominant provinces are non-null and all others are set to zero. The mask of main fisheries [START_REF] Watson | A database of global marine commercial, small-scale, illegal and unreported fisheries catch 1950-2014[END_REF] (chosen as the first 4 deciles) and Exclusive Economical Zones 72 is overlaid on the Bray-Curtis map.

Change in province assemblages

A province assemblage is defined as the assemblage of dominant provinces of each size fraction at a given grid point of the considered ocean. We consider two criteria of change in province assemblage between present day and end of the century conditions.

The first one, more straightforward and less stringent, considers that a province assemblage occurs when a change of dominant province is found in at least one size fraction. In a more stringent way, a change of assemblage is considered significant for (previous section). This threshold corresponds to an idealized case where each dominant province has a probability of one and a change of dominant province is found in only one size fraction. For example, the dominant province assemblage goes from vector (F5,E6,D3,C8,B7,A7) (with the size fractions in decreasing order) corresponding to all temperate provinces to vector (F8,E6,D3,C8,B7,A7). This example corresponds to the replacement of the temperate province of size fraction 180-2000 µm (F5) by the tropico-equatorial province (F8). This criterion allows us to discard assemblage changes for which the changes in probability of presence of dominant provinces are very low. With this criterion, only a small oceanic area is found to have no changes of assemblage (Fig. 4c light blue zones).

Composition of provinces in bacterial diazotrophs and marine copepods

We characterize the composition in Marine Hexanauplia (copepods) and marine diazotrophs of provinces by considering the mean relative abundances of groups of MAGs 32,33 characterized taxonomically (for both copepods and diazotrophs) and by size

Driver analysis

To assess the relative importance of each driver in province changes, the methodology from Barton et al. [START_REF] Michelangeli | Probabilistic downscaling approaches: Application to wind cumulative distribution functions[END_REF] is adopted. For a set N of n provinces (individual provinces or all provinces together), the probability of presence of each province is recalculated for present day conditions except for driver d (from the set of drivers D) for which the end of the century condition is used ( ). The set of driver D can be either all drivers (Fig. 5a,c) or all drivers except SST (Fig. 5b). The relative importance of driver d at a given grid point for the set of N of provinces is computed as follows:

(4) is computed at grid cells where and calculated with either the set of all drivers (Fig. 5a,c) or all drivers except SST (Fig. 5b). When RI(d) is calculated for individual provinces (Fig. 5c and Supplementary Fig. 24d,e), it is computed only at grid cells where and the concerned province is either dominant in present day and/or end of century conditions. For each plankton size class, indexes of presence enrichment (equation ( 1)) for 713 genomes of eukaryotic plankton 38 in corresponding provinces are clustered and represented in a color scale.

Signature genomes (see Methods) are found for almost all provinces, their number and taxonomies are summarized (detailed list in Supplementary Table 6).

  characterize the composition of provinces using 198 environmental genomes detected in at least 5 samples and annotated as Marine Hexanauplia of clade A and B. They are further divided into five subgroups based on their differential abundances across the large size classes (>5 μm): 8 mesozooplankton clade A, 19 unclassified clade A, 71 microzooplankton clade A, 10 unclassified clade B and 90 microzooplankton clade B.

  provinces11 are used to define environmental niches: sea surface temperature (SST), salinity (Sal); dissolved silica (Si), nitrate (NO 3 ), phosphate (PO 4 ), iron (Fe), and a seasonality index of nitrate (SI NO 3 ). With the exception of Fe and SI NO 3 , these parameters are extracted from the gridded World Ocean Atlas 2013 (WOA13)34 .Climatological Fe fields are provided by the biogeochemical model PISCES-v2[START_REF] Aumont | PISCES-v2: An ocean biogeochemical model for carbon and ecosystem studies[END_REF] . The seasonality index of nitrate is defined as the range of nitrate concentration in one grid cell divided by the maximum range encountered in WOA13 at the Tara Oceans sampling stations. All parameters are co-located with the corresponding stations and extracted at the month corresponding to the Tara Oceans sampling. To compensate for missing physicochemical samples in the Tara Oceans in situ data set, climatological data (WOA) are preferred. The correlation between in situ samples and corresponding values extracted from WOA are high (r 2 : SST: 0.96, Sal: 0.83, Si: 0.97, NO 3 : 0.83, PO 4 : 0.89). In the absence of corresponding WOA data, a search is done within 2° around the sampling location and values found within this square are averaged.

(

  for copepods) in each province for all size fractions except the viral one (for diazotrophs) and for size fractions >5 μm (for Marine Hexanauplia). Marine Hexanauplia are annotated taxonomically (either as belonging to clade A (109 MAGs) or clade B (105 MAGs)) from which 198 are found in at least 5 samples that we use (98 clade A and 100 clade B). To attribute a preferential size class to these MAGs, mean relative abundances over all sites of each of them are compared across size fraction 180-2000, 20-180 and 5-20 μm using Welch ANOVA 73 (p-value<0.05). When the Welch ANOVA test is significant the MAG is either annotated as mesozooplankton (when most abundant in size class 180-2000 μm) or microzooplankton (when most abundant in size class 20-180 or 5-20 μm). When the Welch ANOVA test is not significant, the MAG is annotated as unclassified. Then for each group of MAG (defined by the preferential size class plus the clade), the mean of the sum over the MAGs from this group is calculated in each province to characterize the province. The same procedure is applied for 27 out of 48 bacterial diazotrophs (present in at least 5 samples) from the prokaryotic MAG collection 33 distinguishing groups at the phylum level (Gammaproteobacteria n=8, Cyanobacteria n=8, Deltaproteobacteria n=2, Alphaproteobacteria n=4, Planctomycetes n=3, Verrucomicrobiota n=2). Finally, significant differential composition between provinces are annotated using Holm corrected 74 pairwise Mann-Whitney U test 75 (p<0.05 for significance) comparing the abundance distributions of each group of either diazotrophs or copepods between each pair of provinces from a same size class.

Fig. 1 |

 1 Fig. 1 | Eukaryotic signature genomes of provinces of eukaryote enriched size classes.

Fig. 2 |Fig. 3 |

 23 Fig. 2 | Global biogeographies of size-structured plankton provinces projected on WOA2013 dataset. (a) Metazoans enriched (180-2000 µm) (b) Small metazoans enriched (20-180 µm) (c) protist enriched (5-20 µm) (d) protist enriched (0.8-5 µm) (e) Bacteria enriched (0.22-3 µm) (f) Viruses enriched (0-0.2 µm). (a-f) Dotted areas represent uncertainty areas where the delta of presence probability of the dominant province and an other (from the same size fraction) is inferior to 0.5. Simple biogeographies are observed in large size fractions (>20 µm) with a partitioning in three major oceanic areas: tropico-equatorial, temperate and polar. More complex geographic patterns and patchiness are observed in smaller size fractions with the distinction of pacific equatorial provinces and provinces associated with oligotrophic tropical gyres

Fig. 4 |

 4 Fig. 4 | (a) Bray-Curtis dissimilarity index map comparing present day with end of the century projections of dominant provinces. Maps of trans-kingdom assemblage reorganization of dominant provinces (b) and with a criterion of significance (c). (a) Bray-Curtis dissimilarity index (equation (3)) is calculated by integrating all the dominant provinces presence probabilities over the six size fraction. Most important changes appear in subtropical, temperate and subpolar regions. These changes are due to the displacement of tropical and temperate provinces towards the pole but also the geographical decoupling between large and small size plankton. The mean change in niche dissimilarity index is 0.25. (b) An assemblage is the combined projected presence of the dominant province of each size class. Assemblage reorganization (present day versus end of the century) is either mapped on all considered oceans or with a criterion on the Bray-Curtis dissimilarity index (BC>1/6, see Materials and Methods) (c). Depending on the criterion from 60.1% (b, dark blue) to 48.7% (c, dark blue) of the oceanic area is projected to change of assemblage. New assemblages are expected to appear in 2090 (purple+blue) whereas some 2006 specific assemblages are projected to disappear (red+blue). New assemblages as well as lost assemblages are mostly found in temperate, subtropical and tropical regions where most of the rearrangements are projected.

Fig. 5 | 61 1106Supplementary Fig. 21 |

 56121 Fig. 5 | Map of most impacting drivers on dominant province changes (a), most impacting driver without considering temperature change (b) and relative importance of the drivers in the different size fractions (c). (a) Temperature appears as the top impacting driver on the majority of the projected ocean with a significant change of province (Fig. 4). (b) Salinity and dissolved phosphate are found to be the second and third drivers of province reorganization notably at tropical and subpolar latitudes. Note the importance of nitrate at temperate southern latitudes. (c) Temperature is found to be the most important driver for all size classes but has a more important impact in large size classes (>20 μm). Nutrients have on average a relatively more important impact in small size classes in driving province reorganization.

  Mediterranean Sea. Ecography (Cop.). (2018). doi:10.1111/ecog.02434 24. Beaugrand, G. et al. Prediction of unprecedented biological shifts in the global ocean. Nat. Clim. Chang. 9, 237-243 (2019). 25. McMahon, K. W., McCarthy, M. D., Sherwood, O. A., Larsen, T. & Guilderson, T. P. Millennial-scale plankton regime shifts in the subtropical North Pacific Ocean. Visualizing structure and transitions in high-dimensional biological data. Nat. Biotechnol. (2019). doi:10.1038/s41587-019-0336-3 37. van Vuuren, D. P. et al. The representative concentration pathways: An overview. Polovina, J. J., Dunne, J. P., Woodworth, P. A. & Howell, E. A. Projected expansion of the subtropical biome and contraction of the temperate and equatorial upwelling biomes in the North Pacific under global warming. ICES J. :10.1093/gigascience/giw018 49. Busseni, G. et al. Meta-Omics Reveals Genetic Flexibility of Diatom Nitrogen Transporters in Response to Environmental Changes. Mol. Biol. Evol. (2019).

	(1)
	Clim. Change (2011). doi:10.1007/s10584-011-0148-z
	38. Marinov, I. et al. North-South asymmetry in the modeled phytoplankton doi:10.1093/molbev/msz157
	Science (80-. ). (2015). doi:10.1126/science.aaa9942 community response to climate change over the 21st century. Global 50. Grand, M. M. et al. Developing autonomous observing systems for micronutrient
	26. Pinsky, M. L., Worm, B., Fogarty, M. J., Sarmiento, J. L. & Levin, S. A. Marine Biogeochem. Cycles (2013). doi:10.1002/2013GB004599 trace metals. Frontiers in Marine Science (2019). doi:10.3389/fmars.2019.00035
	taxa doi:10.1126/science.1239352 track local climate 27. Bopp, L. et al. Multiple stressors of ocean ecosystems in the 21st century: velocities. Science (80-. ). (2013). 39. Science (80-. ). 368, eaay9514 (2020).
	bioRxiv 2020.10.15.341214 (2020). doi:10.1101/2020.10.15.341214 45. Gaby, J. C. & Buckley, D. H. A comprehensive evaluation of PCR primers to
	33. Delmont, T. O. et al. Heterotrophic bacterial diazotrophs are more abundant than amplify the nifH gene of nitrogenase. PLoS One 7, (2012).
	their cyanobacterial counterparts in metagenomes covering most of the sunlit 46. Iudicone, D. Some may like it hot. Nature Geoscience (2020).
	ocean. bioRxiv 2021.03.24.436778 (2021). doi:10.1101/2021.03.24.436778 doi:10.1038/s41561-020-0535-z
	34. Boyer, T. P. et al. WORLD OCEAN DATABASE 2013, NOAA Atlas NESDIS 72. 47. Gorsky, G. et al. Expanding Tara Oceans Protocols for Underway, Ecosystemic
	Sydney Levitus, Ed.; Alexey Mishonoc, Tech. Ed. (2013). doi:10.7289/V5NZ85MT Sampling of the Ocean-Atmosphere Interface During Tara Pacific Expedition

Projections with CMIP5 models. Biogeosciences (2013). doi:10.5194/bg-10-6225-2013 28. Thomas, M. K., Kremer, C. T., Klausmeier, C. A. & Litchman, E. A global pattern of thermal adaptation in marine phytoplankton. Science (80-. ). (2012). doi:10.1126/science.1224836 29. Busseni, G. et al. Large scale patterns of marine diatom richness: Drivers and trends in a changing ocean. Glob. Ecol. Biogeogr. (2020). doi:10.1111/geb.13161 30. Ibarbalz, F. M. et al. Global Trends in Marine Plankton Diversity across Kingdoms of Life. Cell (2019). doi:10.1016/j.cell.2019.10.008 31. Hutchinson, G. E. Concludig remarks. Cold Spring Harb. Symp. Quant. Biol. (1957). 32. Delmont, T. O. et al. Functional repertoire convergence of distantly related eukaryotic plankton lineages revealed by genome-resolved metagenomics. 35. Sunagawa, S. et al. Tara Oceans: towards global ocean ecosystems biology. Nature Reviews Microbiology (2020). doi:10.1038/s41579-020-0364-5 36. Moon, K. R. et al. Mar. Sci. (2011). doi:10.1093/icesjms/fsq198 40. Flombaum, P., Wang, W. L., Primeau, F. W. & Martiny, A. C. Global picophytoplankton niche partitioning predicts overall positive response to ocean warming. Nat. Geosci. (2020). doi:10.1038/s41561-019-0524-2 41. Richardson, A. J. In hot water: zooplankton and climate change. ICES J. Mar. Sci. 65, 279-295 (2008). 42. Wrightson, L. & Tagliabue, A. Quantifying the Impact of Climate Change on Marine Diazotrophy: Insights From Earth System Models . Frontiers in Marine Science 7, 635 (2020). 43. Luo, Y.-W. et al. Database of diazotrophs in global ocean: abundance, biomass and nitrogen fixation rates. Earth Syst. Sci. Data 4, 47-73 (2012). 44. Zehr, J. P. & Capone, D. G. Changing perspectives in marine nitrogen fixation. (2016-2018). Front. Mar. Sci. (2019). doi:10.3389/fmars.2019.00750 48. Istace, B. et al. de novo assembly and population genomic survey of natural yeast isolates with the Oxford Nanopore MinION sequencer. Gigascience (2017). doi

Table 2 )

 2 are used to project environmental niches under greenhouse gas concentration trajectory RCP8.5 37 .

	Environmental drivers are extracted south of 60° north for present day (2006-2015) and
	end of century (2090-2099) conditions for each model and the multi-model mean is
	computed. A bias correction method, the Cumulative Distribution Function transform,
	CDFt 56 , is applied to adjust the distributions of SST, Sal, Si, NO 3 and PO 4 of the multi-
	model mean to the WOA database. CDFt is based on a quantile mapping (QM)
	approach to reduce the bias between modeled and observed data, while accounting for
	climate change. Therefore, CDFt does not rely on the stationarity hypothesis and

present and future distributions can be different. CDFt is applied on the global fields of the mean model simulations. By construction, CDFt preserves the ranks of the simulations to be corrected. Thus, the spatial structures of the model fields are preserved.

Environmental niche models: training, validation and projections

Provinces with similar metagenomic content are retrieved from Richter et al.

11

. From a total of 48 initial provinces, 10 provinces are removed either because they are represented by too few samples (7 out of 10) or they are found in environments not resolved by ESMs (e.g. lagoons of Pacific Ocean islands, 3 out of 10). This narrows down the number of samples from 644 to 595 metagenomes. Four machine learning methods are applied to compute environmental niches for each of the 38 provinces:

Table 1 | Global statistics of covered areas and province changes and transitions

 1

. From 12% to 31% of the total covered area is estimated to be replaced by a different province at the end of the century compared to present day depending on the size fraction. In total, considering all size fractions this represents 60% of the total covered area with at least one predicted change of dominant province across the six size fractions.

Table 1 | Earth System models used to compute the mean model. 1160 1161

 1 64index of Nitrate (SI NO 3 ) have the lowest median relative importance in (a). Numbers above violins are median (and mean) relative influences. 67 size fraction 0.8-5. Considering only latitudes between the two tropics, changes are mainly driven by decreases in PO 4 (24%) in addition to SST (27%) (overall 34% STT and 20% PO 4 ).. CC-BY-NC-ND 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is madeThe copyright holder for this preprint (which this version posted July 22, 2021. ; https://doi.org/10.1101/2020.10.20.347237

	bioRxiv preprint	doi:

Table 2 | Genomic provinces climatic annotations and oceanic surfaces covered in present day (2006-15) and at the end of the century (2090-990).

 2

		Fraction (µm)	SST	Sal	Si	NO3	PO4	Fe	SI NO3
		180-2000	27,7	16,8	10,1	10,3	11,4	14	9,7
		20-180	29,7	13,3	8,9	10,7	16,3	13,2	7,9
		5-20	21,5	14,2	8,4	12,7	13,7	12,6	17
	Niche definition	0,8-5	22,6	17,3	9,2	12,3	13,4	14,5	10,8
		0,22-3	23,9	13,7	10,5	13	12,5	14,2	12,2
		0-0,2	27,1	16,4	9,3	9,5	17,9	12,1	7,6
		all	25	15,5	9,5	11,6	14,2	13,6	10,8
		180-2000	75,2	8,4	2,5	2,6	4,1	5,8	1,4
		20-180	71,7	4,2	7,3	3,9	8,7	3,3	1,0
		5-20	57,9	6,1	4,4	6,6	9,0	6,5	1,0
	Climate change	0,8-5	54,9	13,4	7,2	5,5	9,7	5,9	3,6
		0,22-3	42,6	9,9	13,5	6,4	9,1	15,2	3,3
		0-0,2	42,3	18,9	6,6	6,4	12,1	10,8	2,9
		all	51,5	10,3	8,7	6	11	8,7	3,8
	Supplementary

Table 3 | Summary table of the relative importance of each environmental driver in niche definition and in driving geographical reorganization in response to climate change.

 3 Note that in both cases (niche definition and climate change), the row 'all' is not the mean over the size fractions. In the case of niche definition, this is due to a different number of niches in each size class. In the case of climate change, relative influence is either calculated for single provinces at a given grid point then recalculated for individual size class or calculated for all provinces together (row 'all').

System Models outputs, Tilla Roy for preparation of the data, LAGE (Laboratoire d'Analyses Génomiques des Eucaryotes, CEA) members for stimulating discussions on this project, Mahendra Mariadassou, Sakina Dorothée Ayata and Bruno Hay Mele for discussions on statistics and climate envelope models, Laurent Bopp for initial discussions on this project and on climate models and Noan Le Bescot (TernogDesign) for help with the Figures. We thank all members of the Tara Oceans consortium for maintaining a creative environment and for their constructive criticism. Tara Oceans would not exist without the Tara Ocean Foundation and the continuous support of 23 institutes (https://oceans.taraexpeditions.org/). This article is contribution number XX of Tara Oceans.

Acknowledgments PF was supported by a CFR doctoral fellowship and the NEOGEN impulsion grant from the Direction de la recherche fondamentale (DRF) of the CEA. This study received funding from the European Union's Horizon 2020 Blue Growth research and innovation program under grant agreement No 862923 (project AtlantECO). We thank the LSCE (Laboratoire des Sciences du Climat et de l'Environnement, CEA) for providing Earth

Competing interests

The authors declare no competing interests.

Author contributions

PF, OJ and MG conceived the study. MV wrote the bias correction algorithm. PF computed the results, compiled and analyzed the data. PF wrote the initial draft of the paper. JL, OJ and MG conducted a preliminary study. PF, OJ, MG, MV, TD, DI, and PW discussed the results and contributed to write the paper.

Online content

All data and codes used are available at http://www.genoscope.cns.fr/tara/SourceCodes/NCLIM-20102618_codes.tar.gz

Materials and methods

Genomic provinces of plankton

Environmental niches are computed for trans-kingdom plankton genomic provinces from Richter et al. 11 . They consist of the clustering of metagenomic dissimilarity matrices (based on the amount of DNA k-mers shared between pairs of samples) from 6 available size fractions with sufficient metagenomic data from the Tara Oceans dataset.

The six size fractions (0-0.2, 0.22-3, 0.

Genome signature of the provinces

We analyzed the distribution of 713 eukaryotic and 1888 prokaryotic genomes 32,33 within the genomic provinces. These genomes are Metagenome-Assembled Genomes (MAGs) obtained from Tara Oceans metagenomes. For each size class, we select MAGs that are present (according to a criteria defined in Delmont et al. 32 ) in at least 5 samples. We computed an index of presence enrichment of MAGs within provinces as the Jaccard index 51 , defined as follows:

Supplementary information 1 | Niche models of genomic provinces

To compute and test the validity of realized environmental niches, we train four machine learning techniques to probabilistically associate genomic provinces with environmental data: sea surface temperature, salinity, three macronutrients (dissolved silica, nitrate and phosphate), one micronutrient (dissolved iron) plus a seasonality index of nitrate.

The four machine learning techniques are Gradient Boosting Machine (gbm) 1 , Random Forest (rf) 2 , fully connected Neural Networks (nn) 3 and Generalized Additive Models (gam) 4 . Following a cross-validation framework, a valid environmental niche is obtained for 27 out of 38 initial provinces (71%) comforting their definition and covering 529 samples out of 595 (89%, Supplementary Fig. 2). Rejected provinces contain relatively few stations (mean of 6 ± 2.6 versus 19 ± 15.3 for valid provinces, p-value<10 -3

Wilcoxon test 5 ). For spatial and temporal extrapolations of the provinces presented below, we use the ensemble model approach 6 that considers mean predictions of machine learning techniques.

Supplementary information 2 | Integrated biogeographies using the PHATE algorithm

We use the PHATE dimension reduction algorithm 7 to combine all provinces for all size classes into a single consensus biogeography revealing 4 or 7 robust clusters (Supplementary Fig. 5, methods). The 4 cluster consensus biogeography is mainly latitudinally organized distinguishing polar, subpolar, temperate and tropico-equatorial regions. The 7 cluster consensus biogeography distinguishes the equatorial pacific upwelling biome and three subpolar biomes that most likely reflect the chemico-physical structuring of the Southern Ocean and known polar fronts 8 (red lines Fig. 2h). However, learning data are scarcer south of 60°S so these extrapolations need to be taken with caution.

PHATE is also used for modeled projections into two comparable consensus maps, one for present day and one for the end of the century (Supplementary Fig. 17). Some particularly visible patterns of geographical reorganization are common to several or even all size fractions and visible when comparing the two consensus maps (Supplementary Fig. 17 compared to Fig. 3a-d and Supplementary Fig. 14). For example, the tropico-equatorial and tropical provinces expand in all size fractions and the provinces including the pacific equatorial upwelling shrink for size fractions smaller than 20 µm.

Supplementary information 3 | Comparison of the biogeographies with existing partitions of the oceans

Previous ocean partitioning either in biomes [9][10][11] or biogeochemical provinces (BGCPs) 9,10 are based on physico-biogeochemical characteristics including SST [9][10][11] , chlorophyll a [9][10][11] , salinity [9][10][11] , bathymetry 9,10 , mixed layer depth 11 or ice fraction 11 .

Considering three of these partitions as examples we notice differences with our partitions (Supplementary Fig. 8- and with a seasonal variability of the frontiers 10 (Supplementary Fig. 8). Some of these transitions are very large and match entire BGCPs, for example in subtropical North Atlantic and subpolar areas where high annual variations are well known 10 .

Supplementary information 4 | Expansion and shrinkage of provinces in response to climate change

To quantify patterns of expansion or shrinkage of the provinces, we calculate the surface covered by the dominant provinces weighted by probabilities of presence (Supplementary Fig. 18, Supplementary Table 2). In this way, dominant provinces are defined on 100% of the surface ocean (327 million km 2 ) but their presence probabilities correspond to the equivalent of 45 to 74% (due to sampling variability and niche overlaps) of the surface ocean depending on the plankton size fraction (Table 1