Major restructuring of marine plankton assemblages under global warming

Marine phytoplankton and zooplankton form the basis of the ocean’s food-web, yet the impacts of climate change on their biodiversity are poorly understood. Here, we use an ensemble of species distribution models for a total of 336 phytoplankton and 524 zooplankton species to determine their present and future habitat suitability patterns. For the end of this century, under a high emission scenario, we find an overall increase in plankton species richness driven by ocean warming, and a poleward shift of the species’ distributions at a median speed of 35 km/decade. Phytoplankton species richness is projected to increase by more than 16% over most regions except for the Arctic Ocean. In contrast, zooplankton richness is projected to slightly decline in the tropics, but to increase strongly in temperate to subpolar latitudes. In these latitudes, nearly 40% of the phytoplankton and zooplankton assemblages are replaced by poleward shifting species. This implies that climate change threatens the contribution of plankton communities to plankton-mediated ecosystem services such as biological carbon sequestration.

Month-to-month changes in composition were estimated through the Jaccard dissimilarity index and its two additive components (i.e. turn-over and nestedness). Jaccard's dissimilarity index is quantitative and scales between 0 and 1, with a value of 1 indicating a full re-shuffling in community composition. Lower Jaccard index values thus indicate the regions where month-to-month changes in species composition are weaker throughout the year. Contour lines indicate the 0.25 isopleth in mean annual Jaccard dissimilarity index. Ordinary linear regressions were used in (B)-(D)-(F) to test the strength of the linear relationships between mean annual species richness and month-to-month changes in species composition. The linear fits evidence the negative relationships observed between annual mean species richness and monthly changes in species composition. The plots indicate that the stronger seasonality in environmental conditions (and especially temperature) that prevails in temperate latitudes leads to a seasonal succession of species that prevents the establishment of a high species richness on the annual scale. Polar regions display both low mean annual species richness and low mean annual changes in month-to-month composition since their annual species richness is limited by cold temperatures. Consequently, these regions depart from the overall trend of decreasing annual richness with higher mean month-to-month changes in composition. Fitting an ordinary linear regression without accounting for polar regions reinforces the explanatory power and the significance of the linear model. N = 35170 grid cells for all maps and ordinary linear regressions.
Supplementary Note 1: Analyzing the skills of species distribution models (SDMs) for various sets of environmental predictors of decreasing complexity (i.e. increasingly lower number of predictors), and justification of final choice of predictors set for modelling the phyto-and zooplankton species.
Supplementary Note 1-1: Distribution of the adjusted R 2 values of the models used to explore the explanatory power of 14 sets of predictors for predicting phytoplankton species distributions (presence/background data). Distributions are shown for the three main phytoplankton groups. Boxplots (A) display the interspecific variations in R 2 based on Generalized Linear Models (GLM), while (B) display the interspecific variations in R 2 for Random Forest (RF) models. The lower, middle and upper boundaries of the boxplots correspond to the 25th, 50th and 75th percentiles. The lower and upper whiskers extend from the hinges to the lowest or largest value no further than 1.5*IQR (inter-quartile range) from the lower and upper hinges. The sample size (N species modelled) of each boxplot is N = 163 for Diatoms, N = 157 for Dinoflagellates, N = 24 for Haptophytes.
The sets of predictors are as follows:

A B
The sets of predictors are as follows:  To select the sets of predictors to be used in the SDMs, preliminary Generalized Linear Models (GLM) and Random Forest (RF) models were performed for each species with several predictors sets of decreasing complexity (from 10 to 5). The explanatory power (adjusted R 2 ) of the models and the ranking of predictors within each set were extracted and their changes in distribution between the various sets were examined. This was done for total phytoplankton and total zooplankton separately, and then between their respective constituting groups: Bacillariophyta (Diatoms), Dinoflagellata (Dinoflagellates) and Haptophyta (Coccolithophores) for phytoplankton, and Copepoda, Chaetognatha, Pteropoda, Malacostraca, Jellyfish, Chordata and Foraminifera for zooplankton. This way, we identified the most important predictors for modelling the species distributions and to evaluate if a decrease in the models' skill was linked to the removal of certain variables. Group patterns allowed to test whether different groups differed in their main environmental drivers. For phytoplankton species, 14 sets of variables were examined. The first nine aimed to test the impact of: (i) alternative choices between variables that were identified as collinear (Wind vs. MLPAR vs. MLD+PAR); (ii) progressively discarding variables that presented lower ranks (logEKE, Si * ); and (iii) choosing logNO3 over logSiOH4, two variables embodying global macronutrients availability and that presented relatively high correlation coefficient (⍴ = 0.59). The last five sets of predictors (10)(11)(12)(13)(14) allowed to test the impact of alternatively removing those variables that presented relatively high ranks in the previous sets: SST, dSST, N * , logSiOH4, logChl, PAR. Similarly, 15 sets of variables were tested for zooplankton. The first ten aimed to test the impact of: (i) choosing Wind over MLPAR or over MLD+PAR; (ii) selecting PAR over MLD; (iii) discarding Si * , N * , logEKE; and (iv) choosing logNO3 over logSiOH4 (⍴=0.64). The last five sets of predictors (11)(12)(13)(14)(15) aimed to test the impact of alternatively discarding the top 5 predictors: SST, dSST, dO2, logSiOH4, and logChl.
All the resulting R 2 and predictors ranks are shown in the Supplementary Note 1-1 above. Pairwise non parametric variance analyses (Kruskal-Wallis tests) were carried out to test for significant increase/decrease in the models' R 2 between pairs of predictors set. We now describe those results to clarify our final choices of predictors sets mentioned in the main Methods section. First, we describe the variations in the models' R 2 for the phytoplankton species (and their three main groups) across the predictors sets. Second, we summarize the overall predictors ranking for all phytoplankton and their three groups. The results relative to the zooplankton species will be described right after.

Phytoplankton
The first three sets allowed to test between sets including Wind (Set 1), or MLPAR (Set 2) or MLD+PAR (Set 3). No differences in R 2 were detected for Diatoms and Dinoflagellates, but R 2 were higher for Haptophyta when using Set 3 so the latter was kept. Sets 4 and 5 tested the differences between PAR and MLD. Set 4 exhibited significantly higher R 2 than Set 5 for Diatoms only. So PAR was preferred over MLD. Sets 3 and 4 tested the effect of removing Si*. No changes in R 2 were observed when removing Si* so it was discarded for the next tests because we favoured more parsimonious sets. Sets 5 and 6 tested whether using logNO3 (Set 5) over logSiOH4 (Set 6) lead to higher models of explanatory power. We found no significant changes in R 2 between the two sets, so either one of the two predictors will be used in the alternative final predictors sets. Comparing Sets 4, 5 and 6 together also showed that using one of the two or both at the same time does not change R 2 distributions. Sets 6 and 7 tested the effect of including logEKE. No significant differences in R 2 were found so logEKE was discarded for phytoplankton. This was also supported by testing changes in R 2 between Sets 8 (with logEKE) and 9 (without logEKE). Comparing Set 7 to Set 9, and Set 6 to Set 8 allowed to test the effect of removing N*. A significant drop in R 2 was found for Dinoflagellates (R 2 did not vary for the two other groups) so N* was kept. Comparing Set 7 to the last five sets (10)(11)(12)(13)(14) allowed to test the effect of removing those five variables that seemed essential from the previous tests (and the variable importance rankings): SST (Set 10), dSST (Set 11), logNO3/logSiOH4 (Set 12), logChl (Set 13) and PAR (Set 14). According to GLM, removing either SST, logChl and logNO3/logSiOH4 decreased model skill significantly for Diatoms. Similarly, the models built for Dinoflagellates exhibited significant drops in R 2 when discarding SST and logChl. For Haptophyta, removing PAR and SST lead to models of lower explanatory power. According to RF models, removing dSST decreased model skill for all groups.
Examining the distribution of the predictors' relative importance in the GLM and RF allowed us to establish the following ranking of variables for phytoplankton (rankings within phytoplankton groups are not explicitly detailed here but are shown in Supplementary Note 1-3,4): • GLM: SST >> N* > logChl > logSiOH4 > Wind > dSST/logNO3 > logEKE/MLPAR/ Si* > MLD > PAR • RF: SST > N* >> logChl > dSST > logSiOH4/Wind > PAR > logEKE > logNO3 > Si* > MLPAR/MLD Consequently, in combination with the results described above, four final sets of predictors chosen for the SDMs were: Zooplankton Like for phytoplankton, the first three sets allowed to test between sets including Wind (Set 1), MLPAR (Set 2) and MLD+PAR (Set 3). R 2 slightly increased from Set 1 to Set 3 for all groups, but the sharpest increase was observed for Foraminifera according to GLM and RF. Comparing Sets 4 (with PAR instead of MLD) and 5 (MLD instead of PAR) enabled us to identify PAR as the main driver of these variations in model skill, as removing PAR decreased R 2 . Therefore, PAR was kept as a potential predictor as it seemed key to model Foraminifera species. Comparing Set 3 to Sets 4 and 5 allowed to test the effect of removing Si* as the latter is absent in Sets 4 and 5 but not 3. A slight but significant decrease in the R 2 of GLM was observed for Copepoda, Jellyfish and Foraminifera. Therefore we kept it in the list of potential predictors and waited to further examine its relative importance in the models.
Comparing Set 4 to Set 6 allowed to test the effect of removing PAR. No significant changes in R 2 were found so PAR was discarded to favour predictor parsimony.
Similarly, MLD was also removed since comparing Set 5 to Set 6 (no MLD) showed that it did not improve the models R 2 .
Comparing Set 6 to Set 7 allowed to test the effect of choosing logNO3 (Set 7) over logSiOH4 (Set 6). The GLM built for the Chaetognatha showed a slight decrease in R 2 when choosing logNO3. But this result was not confirmed by the RF that showed a slight increase in R 2 for Chordata when picking nitrates over silicates. Like for phytoplankton, both macronutrients variables were kept and could be alternatively chosen in the final predictors sets. Comparing Set 6 to Set 8, and Set 9 to Set 10, allowed to test the impact of removing logEKE. GLMs and RF showed a decrease in R 2 for Copepoda and Chaetognatha, so logEKE was kept. Comparing Set 6 to Set 9, and Set 8 to Set 10, allowed to test the impact of removing N*. Neither GLM nor RF showed a change in R 2 so N* did not appear as a key predictor for most of the species. However, we considered keeping it because it could be used as a predictor to embody gradients in nitrates availability when logSiOH4 is used instead of logNO3, without decreasing model skill.
Comparing Set 9 to the last five sets (11)(12)(13)(14)(15) allowed to test the effect of removing those five variables that seemed essential from the previous tests (and the variable importance rankings): SST (Set 11), dSST (Set 12), dO2 (Set 13), logChl (Set 14) and logSiOH4/logNO3 (Set 15). Removing SST significantly decreased the R 2 of GLM for Copepoda, Chaetognatha and Malacostraca. Surprisingly, no changes in the R 2 of the RF models were found. On the contrary, both RF and GLM found significant decreases in models R 2 for Copepoda, Jellyfish, Chaetognatha, Pteropoda, Malacostraca and Foraminifera when removing dSST, thus suggesting this variable is quite relevant for modelling the species' distribution. Removing dO2 significantly decreased the models' R 2 for some groups, but those differed according to the type of model: Malacostraca and Chordata with GLM, Jellyfish, Malacostraca and Foraminifera with RF. Omitting logChl from the predictors weakened the GLM explanatory power for Copepoda and Malacostraca (none with the RF). And finally, removing either logSiOH4 or logNO3 significantly decreased the R 2 of GLM for: Copepoda, Chaetognatha, Jellyfish, Malacostraca and Foraminifera. Again, RF were better able to accommodate this loss than GLMs as no significant decrease in their R 2 could be observed.
Examining the distribution of the predictors' relative importance in the GLM and RF allowed us to rank the variables for total zooplankton (rankings within zooplankton groups are not explicitly detailed here but they are shown in Supplementary Note 1-7,8). GLM and RF provided slightly different rankings, which help explain why the models' R 2 sometimes decreased according to one type of model but not the other. Therefore, we chose to clarify the rankings for both types of algorithms: Please note how the distributions of ranks were much less even between the predictors when using RF instead of GLM, especially for the variables that usually ranked lower than SST. This showed why the R 2 of the RF models were much less affected when removing a predictor other than SST, dSST and dO2.
Consequently, in combination with the results described above, four final sets of predictors chosen for the SDMs were: The ensemble projections of climate change impacts on plankton species richness and composition were summarized using a Principal Component Analysis (PCA; see Methods section) whose principal components were used to compute a (B) quantitative index of the relative severity of climate change impacts on plankton diversity whose distribution across the newly-defined regions was used to rank them in decreasing order. The spatial distribution and the regional variations of (M)-(X) six proxy variables of marine ecosystem services were examined per region: (M)-(N) Oceanic Biodiversity (normalized species richness of the marine megafauna that mainly inhabit the open ocean [1]), (O)-(P) logged mean annual reported and unreported catch rates of small (< 30cm) pelagics fisheries [2], (Q)-(R) mean annual Net Primary Production (NPP; mg C m -2 d -1 ) [3], (S)-(T) mean annual flux of particulate organic carbon (FPOC; mg C m -2 d -1 ) below the euphotic zone [3], (U)-(V) mean annual FPOC/NPP ratio (e ratio) [34], and (W)-(X) mean annual plankton size index (i.e. inverse of the slope of the particles size distribution) [4]. Nonparametric variance analyses (Kruskal-Wallis tests) were performed to test the significance level of the regional variations in plankton diversity changes and proxy variables of marine ecosystem services (0.01 significance threshold). Normality and homoscedasticity tests were performed prior to the Kruskal-Wallis tests. Pairwise post hoc tests of multiple comparisons of mean rank sums were performed to identify which pairs of regions displayed those significant variations. Bonferroni corrections were applied for p-values adjustment.

Supplementary Note 2 -Spatial overlap analysis between climate change impacts on global surface plankton diversity and the current provision of marine ecosystem services
For the ensemble projections of plankton diversity changes, all five variables displayed significant variations across the six regions (all p < 0.001) and post hoc tests found all pairwise comparisons to be significant except for (F) the mean percentage difference in annual phytoplankton species richness between regions 3 and 6. For the proxy variables of marine ecosystem services, all six variables displayed significant variations across the six regions (all p < 0.001) and post hoc tests found all pairwise comparisons to be significant except for the variations in: (N) Oceanic Biodiversity between regions 3 and 4 and between regions 1 and 5; (P) annual small pelagics catch rates between regions 3 and 4; (T) mean annual FPOC between regions 1 and 4; and (V) e ratio between regions 1 and 2. Because the majority of low catch rates for the Southern Ocean are due to a lack of reported data, region 5 was discarded when comparing mean annual fisheries catch rates.

Supplementary Note 3: List of the characters used to identify the datasets in GBIF and OBIS issued from sediment cores that provide sediment or fossil species records.
Zooplankton groups that are characterized by a calcareous or a siliceous shell (e.g. Foraminifera, Thecosomata, and Radiolaria) are often found as fossils in sediment cores where they are used as paleo-indicators of past environmental conditions. As a result, many of these groups' biological observations are archived in online biodiversity repositories come from deep sediment cores. This is particularly relevant for Foraminifera with the advent of global sediment core databases such as CLIMAP, MARGO or ForCenS [5]. In the context of our study, that aims to model extant species diversity from observations collected in the surface ocean, it is critical to remove such observations from the OBIS and GBIF datasets we used.
To do so, we first examined the names of the datasets we downloaded from OBIS ('bibliographicCitation' column) and GBIF (using the 'datasets' function of the rgbif R package). Then, we excluded those that contained at least of the following keywords: Finally, the names of the remaining datasets were inspected again, and those that were not discarded but still presented equivocal names (e. As shown by the differences between maps (A), (B) and (C), including SSS or Longitude as a predictor leads to peaks in richness being concentrated in the Atlantic Ocean, whereas removing it leads to more even richness values across the basins without significantly decreasing SDMs skill. We interpret these differences in diversity patterns as the result artifactual effect caused by the concentration of zooplankton occurrences in the Atlantic Ocean and an ensuing bias towards higher SSS values in the SDMs. We do not neglect here the potential role of SSS in structuring plankton communities in the surface ocean. But, in light of the biases in environmental space evidenced here that are leading to obvious differences in the resulting species richness patterns, we chose to discard SSS and pCO2 from the list of potential environmental predictors.

Supplementary Note 5: ODMAP protocol
Major restructuring of marine plankton assemblages under global warming.

Focal Taxon
Focal Taxon: Marine plankton (main phytoplankton and zooplankton groups).

Location
Location: Global surface ocean

Assumptions
Model assumptions: our SDMs rely on the following main assumptions: (i) niche conservatism through time; (ii) species distributions are not strongly limited by dispersal at a macroecological scale, an assumption valid for plankton considering the very strong connectivity of ocean basins through surface current on decadal scales, which enables plankton species to display very large spatial ranges; (ii) at the scale of the study, species spatial distributions are primarily shaped by the combinations of environmental factors (not biotic interactions) that define the conditions allowing a species to develop.

Algorithms
Modelling techniques: glm, gam, randomForest, ann Model complexity: To cover range of models types (regression, classification trees, neural networks) and model complexity.
Model averaging: Models projections (habitat suitability indices) were averaged without weights.

Workflow
Model workflow: See online Methods section upon publication of the study.

Software
Software: biomod2 Code availability: All codes available on the GitHib page of the first author: https://github.com/benfabio Data availability: All occurrence data obtained from publicly available datasets.

Data Biodiversity data
Taxon names: 860 Plankton species (336 phytoplankton species and 524 zooplankton species) -see Supplementary Data 2 for the list of species names and their taxonomic classification.
Cleaning: See online extended Methods upon publication of the study (+Supplementary Table  S7).
Absence data: No absence data available.
Validation data: Randomly selected 20% withheld from model fitting.
Test data: Not available.

Model Variable pre-selection
Variable pre-selection: See online extended Methods upon publication of the study.

Model estimates
Coefficients: See Model settings above.
Variable importance: See online extended Methods upon publication of the study.

Model selection -model averaging -ensembles
Model averaging: Arithmetic mean of the models projections of species habitat suitability indices (no weighting by SDM evaluation metric).

Analysis and Correction of non-independence
Spatial autocorrelation: Thinning of species occurrences over 100km.
Nested data: Not applicable.

Threshold selection
Threshold selection: Non applicable. Diversity estimates were based on the species' continuous average habitat suitability indices.

Performance statistics
Performance on training data: TSS, AUC, Kappa

Performance on validation data: TSS, AUC, Kappa
Performance on test data: Not available.

Plausibility check
Response shapes: Examination of SDM-specific mean univariate response curve for every species modelled.
Expert judgement: Examination of SDM-specific mean annual habitat suitability maps for every species modelled.

Prediction Prediction output
Prediction unit: Monthly and mean annual species richness estimated through sum of species habitat suitability indices.
Post-processing: Not applicable.

Uncertainty quantification
Parameter uncertainty: Not applicable.
Scenario uncertainty: Uncertainties in ensemble projections were examined through the standard deviation computed from the ensemble model members projections (n = 16 for contemporary conditions; n = 80 for end of century conditions since 5 Earth System Models were used).
Novel environments: SDMs projections uncertainties (see section above) and non analog climate conditions (i.e. identified and quantified through species-level MESS maps) were all illustrated on the main Figure  On (A) all PFGs score with PC1 > 0 as they all display a global latitudinal diversity gradient in mean annual SR (i.e. increase in SR from the poles towards the equator). Differences between phyto-and zooplankton PFGs appear on PC2 as phytoplankton score with PC2 < 0 while the zooplankton score with PC2 > 0. This discrepancy in PC space stems from the spatial separation of the SR maxima: the SR of phytoplankton PFGs peaks near the equator and in tropical upwellings whereas the SR of zooplankton PFGs (except Chaetognaths) peaks in the subtropics.
On (B) differences between and within phyto-and zooplankton PFGs are clearer: zooplankton PFGs (and Haptophytes) score with PC1 > 0 whereas phytoplankton groups (and Chaetognaths) score with PC1 < 0 and are projected on PC2 > 0 more strongly. PFGs scoring with PC1 > 0 are those presenting decreasing SR in the tropics and increasing SR at higher latitudes (~40-50°N). PFGs scoring with PC1 < 0 and PC2 > 0 are those presenting increased SR in the tropics and temperature latitudes and slighter decrease in SR near the poles. Dinoflagellates do not score with PC1 because they display decreasing SR in tropical upwellings country to Diatoms and

A B
Chaetognaths. Foraminifera is the only group slightly scoring with PC2 < 0 because it displays stronger decrease in SR than other PFGs.
The PFGs constituting our projections of phyto-and zooplankton SR here show distinguishable mean annual SR spatial patterns in the contemporary and future ocean. This supports our assumption that phytoplankton and zooplankton can be separated into PFGs whose diversity responds to environmental variations in different ways. These LDGs were used to further validate our ensemble projections of mean annual phyto-and zooplankton diversity by comparing them to previously published plankton groups LDG that were based on observations. Please note that direct global comparisons to previous studies are made difficult because: (i) previous diversity estimates are often based on single sampling cruises and thus do not span the same spatio-temporal scales as our estimates (e.g. single point measurements vs. 1°x1° monthly estimates); (ii) previous cruises estimate plankton alpha diversity through approaches that are quite different from ours (e.g. high throughput sequencing of environmental DNA; sum of observed taxa accounting for their relative abundances); (iii) the global LDG of some of the PFG have simply not been documented from observations (e.g., Dinoflagellates, Jellyfish, Chordates) which emphasizes the importance of our study. Therefore, it is important to keep in mind that we are here looking for similarities in overall latitudinal patterns. Below, we cite the studies that were used to validate our PFG diversity patterns and briefly highlight the elements that validate our mean annual estimates while discussing discrepancies. . The authors also found Diatom richness to decrease from 30°S to 45°S but then to increase again from 45°S to 60°S. This provides further support to our projections of mean annual Diatom richness which increase from 40°S to > 60°S.

Supplementary
• Dinoflagellates: no previous study; closest estimate is the LDG of "photosynthetic/mixotrophic (P) Protists" from Ibarbalz et al. [8] - Fig. 2A (same as for Diatoms). The authors also found a peak in species richness near the equator and decrease towards poles.
• Haptophyta (mainly Coccolithophores): O'Brien et al. [11] - Fig. 3C & Fig. 6A. The authors modelled global Coccolithophores species richness from in situ observations and monthly climatologies through neural networks. The authors found mean annual Coccolithophores richness to peak in the tropics and decrease towards high latitudes, with notable decreases in richness near tropical upwelling regions, which is very similar to our mean annual estimates.
• Copepods:  Rombouts et al. [12] - Fig. 7. The authors also found Copepod species richness to increase from the poles to the equator with peaks around 20°-30° latitude. The authors found Copepod richness to peak in the subtropical Atlantic Ocean likely because of the overrepresentation of this basin in their dataset (their Fig. 1) and the inclusion of sea surface salinity as one of the main predictors of Copepod diversity (but see Supplementary Note 4 of our study). The authors projected zooplankton diversity to peak in the subtropics as a result of the selection of species based on their relative thermal tolerances. Their modelling approach was validated against the observations of Rombouts et al. [12]. The authors also found zooplankton diversity to decrease >23°C.  Fig. 1L. The authors also found Euphausiids species richness to decrease from the equator to the poles with peaks ~30° latitude. Overall, and despite the much lower spatial resolution of their estimates (880km), both our studies found very similar hotspots in Euphausiid diversity and a decrease in diversity at the equator.
• Chordates (Salps, Doliolids and Larvaceans): no previous study for this PFG. The authors also found Chaetognatha species richness (estimated from species counts in the Indo-Pacific) to decrease from tropics to high latitudes, with peaks in the western equatorial Pacific and west of Australia. The authors also found zooplankton richness to decrease > 23°C. However, their observed patterns in Chaetognatha richness differs from ours in the Indian Ocean. This may be linked to the differences of our approaches but also by the high seasonality in Chaetognatha diversity they documented for this basin (Fig. 7). The authors also found Foraminifera species richness to decrease from the equator to the pole with peaks in the subtropics, whereas we found Foraminifera richness to peak near the equator (i.e. in the tropical western Pacific) though with substantial longitudinal variations. However, the authors also found zooplankton richness to decrease with SST beyond 23°C.
 Tittensor et al. [1] - Fig. 1M. The authors heavily relied on the data of Rutherford et al. [17] so they found similar Foraminifera richness patterns.
 Yasuhara et al. [18] - Fig. 2. The authors also found Foraminifera species richness to decrease from the tropical Atlantic Ocean to the Arctic Ocean. Plus, the authors also found zooplankton diversity to decrease with temperature above 23°C.

Supplementary Note 6-3:
Comparison of the present ensemble projections of contemporary and future changes in mean annual phytoplankton and copepod species richness (SR) to the contemporary and future changes in mean annual species diversity estimates (Shannon index, H') of Ibarbalz et al. [8].
We were able to directly compare their model projections of copepods and photosynthetic-mixotrophic protists (labelled as "Protists (P)" in Ibarbalz et al. [8]) species diversity to our results. These two groups provide the most comparable projections to our present phytoplankton ( Figure 1) and copepod (Supplementary Note 6-2 above) diversity (i.e. mean annual species richness) estimates. In short, the following pairs of fields were compared: -Contemporary mean annual phytoplankton species richness (i.e. from ensembles of habitat suitability; Fig. 1d Before a direct comparison could be made, we normalized both ours and their estimates of contemporary phytoplankton and copepod diversity by their respective maximum values. For each pair of variables, bivariate plots were drawn with our own model estimates on the y axes and Spearman rank correlation coefficients (rho) were computed. This way we evaluate the similarity of the spatial patterns in diversity between the results of Ibarbalz et al. [8] and ours. Then, using the bivariate plots illustrating the amplitude of the future differences in diversity, we identified the regions were our projections agree or disagree on the sign of the response of protist/copepod diversity to future climate changes. Plus, by drawing the 1:1 line of these two plots, we also identified the regions where our changes in diversity are predicted to be larger or weaker than those of Ibarbalz et al. [8], relative to the respective contemporary conditions. All the plots and the results from the correlation analyses are summarized below.

Supplementary Note 6-4:
Comparison between our present estimates (always on the y axis) to those of Ibarbalz et al. [8] (always on the x axis) for (a) contemporary mean annual phytoplankton/ Protists (P) species diversity (species richness SR vs. Shannon diversity index H') and (c) mean difference (future decadecontemporary decade) in mean annual phytoplankton/ Protists (P) species diversity. (b) same as (a) and (d) same as (c) but for copepod species diversity instead of phytoplankton/ Protists (P). Each point corresponds to a 1°x1° grid cell and was colored as a function of its latitude. The dashed line represents the 1:1 line, and the dotted lines show where changes in species diversity are equal to zero.
Results from the Spearman' rank correlation tests associated to the data plotted above: a. Rho = 0.852 ; p-value < 2.  [8] as well as the relative amplitude of the projected differences in species diversity when both estimates agree on the direction of change. Cells in red and green thus correspond to those regions where our projections disagree with those of Ibarbalz et al. [8].
In short, we find significant correlations coefficients for four pairs of variables which indicates that our results are overall in line with those of Ibarbalz et al. [8]. Yet, the correlations vary in strength. Mean annual Protists (P)/phytoplankton and copepod species diversity display a similar patterns between the two studies (all rho > 0.83), indicating that both studies find very similar latitudinal diversity gradients of species diversity for Protists (P)/phytoplankton and copepods for the contemporary ocean. Although the correlation coefficient is relatively weak (rho = 0.248), both studies find Protists (P)/phytoplankton diversity to increase in the future, but they strongly disagree on the response of diversity in high latitudes of the northern hemisphere (increase vs. decrease in our case). The predicted future changes in global copepod species diversity are more similar between the two studies compared to phytoplankton, although Ibarbalz et al. [8] found more regions where copepod diversity is likely to increase in the future.
We would like to underline that there are several major methodological differences, which makes it difficult to pinpoint the reasons behind the differences shown above. Indeed, a key issue besides the broader spatio-temporal coverage of our data, is the way plankton diversity is measured and modelled. Ibarbalz et al. [8] derived Shannon diversity indices (H') from numbers of reads of operational taxonomic units (OTUs) obtained from high throughput 16S and 18S RNA sequencing. Therefore, it is hard to evaluate the similarity of the phytoplankton (or the "Protists (P)") community they sample compared to ours. They directly estimate species diversity at 189 sampling stations and then model it as a response variable through GAMs. Then, the authors use the latter GAMs to project species diversity in space and time as a function of sea surface temperature (SST) and chlorophyll a only. Meanwhile, we estimate species diversity as a property that emerges from stacking several hundreds of different species for which we can reliably model the global habitat suitability patterns. Furthermore, we use a broader range of model types and complexity whereas Ibarbalz et al. [8] rely on GAMs only. Therefore, our respective approaches differ in a multitude of ways: (i) our biological observations span much broader spatial and temporal scales (several decades and all ocean basins vs. two cruises); (ii) species diversity is directly measured and modelled in their case contrary to our approach; (iii) we use a broader range of statistical models that better covers the commonly used range of algorithms and thus actually accounts for this major source of uncertainty in diversity forecasts; (iv) they rely on two environmental predictors (SST and chlorophyll a) to model the diversity of all the plankton groups whereas we made sure to use four different sets of predictors that span several niche dimensions adapted to each trophic level; (v) the baseline and end-of-century periods defined to compute the future environmental fields on which the statistical models are projected on differ too ( The modelled relationship between mean annual zooplankton species richness and temperature (as shown in Figure 2) obtained from the rarefied dataset of zooplankton occurrences. The modelled species richness estimates are lower than those obtained from the unrarefied dataset, as expected from the necessarily lower number of zooplankton species modelled (~50%). Yet, the diversity-thermal energy relationship observed is very similar to the one obtained from the full dataset. As temperature is the main driver of the modelled diversity patterns, this strongly suggests that our main results are not an artefact of sampling biases and that the emergent decrease of zooplankton diversity towards the highest temperatures is not linked to the lower sampling effort observed near the equator. The solid curve illustrates the 3 rd degree polynomial fit that best explains the variation in log(Zooplankton SR) as a function of mean annual available thermal energy. The colored isopleths illustrate the density of ocean grid cells, based on two dimensional kernel density estimates, and were used to highlight the parts of the gradients driving the observed non-linear relationships.
Following the MTE framework, we examined the relationship between the logged annual species richness (SR) and available thermal energy to better understand why phyto-and zooplankton SR respond differently to climate change in the tropical ocean. The SR of both groups displays non-monotonic relationships with thermal energy that depart from the linear slope predicted by the MTE (~|0.32| for phytoplankton; ~|0.65| for zooplankton). For both groups, different modes of SR variations are modelled across the thermal energy gradient as illustrated by the polynomial fits. First, a mode of decreasing SR is found from the coldest temperatures (-1.8°C) to approximately 10.5°C, with slope values steeper than -0.13 for phytoplankton (R 2 = 0.11; p < 0.001), and -0.31 for zooplankton (R 2 = 0.60; p < 0.001) where mean annual sea surface temperature (SST) is < 5°C. This decrease in SR is not predicted by the Metabolic Theory of Ecology (MTE) and can be ascribed to the strong environmental seasonality prevailing in such temperate regimes because it prevents the establishment of a high phyto-and zooplankton SR on the annual scale, as supported by relatively high rates of month-to-month species turn-over (Supplementary Figure 6). Second, a mode of steep SR increase was found to start ~11°C above which the fitted linear slopes are closer to the predictions of the MTE. For phytoplankton, the linear slopes fitted from 11°C to increasingly higher thermal energies range between 0.33-0.40 until ~23°C (R 2 = 0.16-0.58; all p < 0.001), after which a third mode of weaker SR increase takes over, as evidenced by slopes < 0.29. The linear slope best matching the MTE predictions is the one fitted for oceans regions where mean annual SST is > 22°C (slope = 0.33, R 2 = 0.164; p < 0.001). For zooplankton, the linear slopes fitted from 11°C to 13-19°C exceed the slope predicted by the MTE (all slopes range between 0.68-1.00; R 2 = 0.31-0.77; all p < 0.001). The best fit to the MTE was found for the regions where mean annual SST is comprised between 11°C-22°C (slope = 0.66; R 2 = 0.78; p < 0.001). Beyond a maximum SST of 22°C, the linear slopes fitted from 11°C become lower than MTE predictions as zooplankton SR clearly starts decreasing towards highest temperatures. This pattern of decreasing diversity at high temperatures has been reported for other animal clades. This third mode of decreasing zooplankton SR does not match the corresponding mode found for phytoplankton where SR keeps increasing. Such different arrangements of SR to the warmest range of the thermal gradient explain why SR peaks are not co-located for phyto-and zooplankton and why their responses differ in the tropics (Figure 1). In regions where climate change will create conditions of mean annual SST >25°C, zooplankton SR is predicted to decrease abruptly because the novel temperature regime might exceed the thermal tolerances of many species, which lowers habitat suitability. Meanwhile, phytoplankton SR increases nearly linearly when SST is >11°C, and is thus promoted by future warming and very high temperatures remain suitable for most phytoplankton taxa.

Supplementary Note 9
Supplementary Note 9-1: Heatmaps showing the Spearman's rank correlation coefficients computed between the (A) mean annual baseline species richness and (B) the % difference in mean annual species richness (future-baseline) of phyto-and zooplankton as well as their ten main Plankton Functional Groups (PFGs). Total sample size is N = 35023.

A B
Supplementary Note 9-3: Heatmap showing the Spearman's rank correlation coefficients (rho) computed between the proxy variables of ecosystem services (ES), a few key environmental covariates and the mean annual baseline species richness (SR) of phytoand zooplankton, and their ten main Plankton Functional Groups (PFGs). Total sample size is N = 35023.
While we cannot infer mechanistic links based on our modelling approach and from the data at hand, we use the present correlation analysis to assess the strength of associations between the PFGs SR and the proxy variables of ES, on an annual scale. This way, we can test relevant hypotheses that link plankton diversity to ecosystem functioning. Because we rely on non parametric rank correlation coefficients that can detect weak signals as significant, we carefully highlighted with black outlines the strongest correlation coefficients (rho > |0.4|) that we deem interesting for our study as Supplementary Note 9-2: Boxplots showing the distribution of the % difference in mean annual species richness (future-baseline) between the six regions of climate impacts severity for the ten main Plankton Functional Groups (PFGs) studied. The lower, middle and upper boundaries of all the boxplots correspond to the 25th, 50th and 75th percentiles. The lower and upper whiskers extend from the hinges to the lowest or largest value no further than 1.5*IQR (inter-quartile range) from the lower and upper hinges. N = 35023 grid cells for the total sample size (Nregion1 = 2344; Nregion2 = 5954; Nregion3 = 6748; Nregion4 = 8290; Nregion5 = 7637; Nregion6 = 4050). they support potentially robust Biodiversity-ES relationships. All correlations displaying a rho > |0.1| tested as significant (p-value < 0.01).
Please note that the variable used here represent the mean annual conditions of the contemporary offshore global ocean. Oceanic biodiversity quantifies the species richness of higher trophic levels (tunas, sharks, mammals, squids) [1]. FPOCex (mg Carbon m -2 day -1 ) quantifies the amount of particulate organic matter exported beyond the euphoric zone [3]. NPP (mg Carbon m -2 day -1 ) quantifies the rate at which biomass is stored in the phytoplankton and made available to grazers [3]. The e ratio is the FPOCex to NPP ratio and thus quantifies the efficiency of the biological carbon pump in exporting carbon out of the euphotic zone [3]. PSI is an index of surface plankton particles size derived the from the slope of the plankton particles spectrum [4]. Fish catches (tons km -2 year -1 ) estimates the reported and unreported catch rates of small (< 30cm) pelagic fishes over the 1990-2019 period [2]. NO3 and SiO2 (µM) quantify surface macronutrients concentrations (nitrates and silicates respectively). SST (°C) corresponds to sea surface temperature. Chlorophyll (mg C m -3 ) was used as a proxy of surface phytoplankton biomass. Chlorophyll, NO3, SiO2, Fish catches, NPP and FPOCex were log transformed.
Our analysis shows that the SR of PFGs covaries positively with the SR of higher trophic levels, meaning they might respond to similar drivers (i.e. temperature, productivity) in the same way. This also brings further support to the validity of our projections as our models do seem to capture the main processes that drive marine biodiversity. We also find that the efficiency of the biological carbon pump decreases with the SR of all PFG, which supports the view that species-rich and functionally diverse communities better retain the biomass produced within surface layers while species-poor communities tend to leak higher rates of biomass out of the surface ocean. This can also be mediated through rearrangement of community size structure as species rich communities are composed of more numerous smaller species (but see Supplementary Notes 10 & 11). This is supported by the fact that plankton size (i.e. PSI) shows similar correlation patterns to the e-ratio. Haptophyta (i.e. Coccolithophores) SR decreases with plankton size but also silicates. Plus, we evidenced that Haptophyta SR strongly decreases where Diatoms SR peaks (see Supplementary Note 6). This could stem from our models partly capturing competitive effects leading to the dominance of Diatoms over Coccolithophores under conditions of silicates replenishment, though this is hard to evaluate at the scale of our study. The ensuing decrease/increase in Coccolithophores/Diatoms habitat suitability could help explain why plankton side distribution shifts towards larger particles (microphytoplankton). Such competitive effects could also help explain why Diatoms SR and Haptophyta SR display such contrasting responses to future climate change (Supplementary Note 9-2 above). Zooplankton SR decreases with plankton particles size, suggesting that higher densities of large phytoplankton cells might not favour the diversity of grazing PFGs (i.e. Copepods and Malacostraca) as only a fraction of their taxa may be able to graze upon larger cells due to size-based limitations in prey capture. Similarly, we find zooplankton SR to decrease with phytoplankton biomass (i.e. Chlorophyll). This implies that higher food availability may lead to a decrease in SR through competitive exclusion, instead of promoting niche partitioning within the grazing zooplankton. Ultimately, we find phytoplankton SR to decrease with higher nutrients availability this can be interpreted by two mutually non-exclusive hypotheses: (i) species rich communities draw down nutrients concentrations more efficiently, meaning phytoplankton diversity optimizes nutrients use efficiency in the ecosystem which would translate into less biomass/energy leaving the surface ecosystem (corroborated by negative correlations between phytoplankton SR and e-ratio); (ii) less seasonally varying and warmer and nutrients-poor ecosystems (e.g. tropical gyres) sustain species-rich communities as competition under very low nutrients availability leads to niche partitioning that enables the co-existence of many taxa, potentially characterized by diverse functional traits. Under either hypotheses, our results support the existence of important links between key ecosystem functions (i.e. nutrients use efficiency and carbon export) and plankton species richness, and warrants further studies examining at the relationships between ES, species richness, community traits expression and functional diversity. the data at hand, we use the present correlation analysis to assess the strength of associations between Copepod community size structure and the proxy variables of ES. This way, we test relevant hypotheses that link zooplankton size structure and species diversity to ecosystem functioning. Because we rely on non parametric rank correlation coefficients that can detect weak signals as significant, we carefully highlighted with black outlines the strongest correlation coefficients (rho > |0.4|) that we judge interesting for our study as they support potentially robust zooplankton size-ES relationships. All correlations displaying a rho > |0.1| tested as significant (p-value < 0.01).
Please note that the variable used here represent the mean annual conditions of the contemporary offshore global ocean. Oceanic biodiversity quantifies the species richness of higher trophic levels (tunas, sharks, mammals, squids) [1]. FPOCex (mg Carbon m -2 day -1 ) quantifies the amount of particulate organic matter exported beyond the euphoric zone [3]. NPP (mg Carbon m -2 day -1 ) quantifies the rate at which biomass is stored in the phytoplankton and made available to grazers [3]. The e ratio is the FPOCex to NPP ratio and thus quantifies the efficiency of the biological carbon pump in exporting carbon out of the euphotic zone [3]. PSI is an index of surface plankton particles size derived the from the slope of the plankton particles spectrum [5]. SST (°C) corresponds to sea surface temperature. Chlorophyll (mg C m -3 ) was used as a proxy of surface phytoplankton biomass. Chlorophyll, NPP and FPOCex were log transformed.
Overall, we find that Copepod SR decreases with Copepod median body size and body size diversity, suggesting that communities characterized by a larger proportion of largebodied species still comprise smaller bodied ones. Meanwhile, tropical communities are characterized by numerous smaller-bodied species and very few large-bodied species. Interestingly, Copepod SR and size structure also show contrasting correlation patterns with phytoplankton biomass (Chlorophyll), plankton particles size (PSI) and the efficiency of the biological carbon pump (e ratio). Indeed, higher phytoplankton biomass promotes lower Copepod species diversity by favoring larger-bodied species which could outcompete smaller ones as they are able to graze on larger phytoplankton cells (e.g. corroborated by the positive correlation between Copepod median size and PSI). This supports the view that higher primary production leads to competitive exclusion (and not niche partitioning) within the zooplankton community, which in turn leads to a decrease in species diversity. However, we cannot exclude that SST also strongly influences the emerging pattern of Copepod SR through the selection of species with broader thermal tolerances [13]. Yet, both hypotheses are not mutually exclusive. Our results also suggest that Copepod size structure is a major factor in regulating the efficiency of the biological carbon pump but not the amount of carbon exported outside of the euphotic zone (i.e. no correlation with FPOCex). This can be explained by the fact that larger copepod species present many functional traits that tend to scale with body size (e.g. production of larger fecal pellets, deeper of vertical migrations, higher biomass etc. [23][24][25]) and that promote carbon export. Overall, this gives us further confidence that our species diversity and community composition estimates are relevant for studying ecosystem functioning and the provision of marine ecosystem services.
Supplementary Note 11-3: Patterns of % difference (future-baseline) in (A) mean annual Copepod species richness (SR), (B) annual median Copepod body size, and (C) annual median Copepod body size diversity.
As in Supplementary Note 11-1, the future monthly projections of Copepod species habitat suitability (HSI) from each model member (n = 80) were combined with the species-level estimates of average adult female body size measurements to estimate monthly HSI-weighted median Copepod body size for the future surface global ocean. By looking at the % difference between the future and the baseline estimates of Copepod community size structure, we investigate how anthropogenic climate change might reshuffle the size structure of the zooplankton community. More specifically, we here test whether future warming will lead to a decrease in median body size through the replacement of large-bodied Copepod species by smaller ones, a process expected under global warming, particularly towards higher latitudes [24][25][26].
In line with Figures 1 and 3, our ensemble projections show that anthropogenic climate change will trigger the poleward migration of tropical Copepod species which results in large increase in SR in high latitudes combined with large turn-over rates (i.e. species replacement, see Methods) as the warm-water species replace the cold-water ones (already discussed in the main text). Here, we further show that such changes in SR and community composition will also affect the size structure of the Copepod communities. In agreement with our expectations, we find that global warming will decrease the median body size of the Copepod community in high latitudes (except the Southern Ocean), which is driven by the replacement of large-bodied species by smaller ones. Median Copepod body size variance (i.e. body size diversity) follows a similar pattern: future changes in Copepod SR and composition will decrease body size diversity Together with Supplementary Note 11-2, our results confirm our expectations that future changes in plankton diversity will alter ecosystem functioning and the efficiency of the biological carbon pump, particularly at high latitudes (regions 1 and 2 in Figure 4). Our result show that future increases in Copepod SR, and the ensuing turn-over in community composition, will decrease the median body size of the surface zooplankton community which might decrease the efficiency of the biological carbon pump (i.e. the fraction of carbon produced in the surface that is exported below the euphotic zone).