OverviewOur study investigates the patterns and drivers of global marine plankton diversity by simultaneously modeling the spatial distribution of 860 phyto- and zooplankton species, based on the widest and most recent compilations of in situ observations available. These observations were associated with various sets of relevant predictors to train a range of statistical species distribution models (SDMs) on a monthly resolution. The SDMs were used to estimate contemporary and future levels of global surface species richness (SR) for total plankton, phytoplankton and zooplankton. We explore how, and why, global phyto- and zooplankton SR and community composition are affected by future climate change under the RCP8.5 scenario of greenhouse gas (GHG) emissions. We also summarize regional patterns of climate change impacts on plankton diversity by clustering the global ocean and examine how hotspots of climate change impacts might overlap with the current provision of marine ecosystem services. All data manipulation and analyses were performed under the R programming language39. The R packages used are mentioned below in their corresponding section.Plankton species observationsFirst, to model global, open ocean plankton diversity from species-level field observations, comparable datasets of phytoplankton and zooplankton occurrences (i.e., presences) had to be compiled. We refer to as open ocean all those regions where the seafloor depth exceeds 200 m. We made use of the large dataset of phytoplankton occurrences recently compiled by Righetti et al.11. For zooplankton, a new dataset was compiled following the same methodology. Both occurrence datasets were based on publicly available data from online biodiversity repositories, as well as some additional published datasets. The R packages mainly used for implementing the datasets are those constituting the tidyverse package.Phytoplankton occurrences
For the phytoplankton occurrences used here, Righetti et al.40 compiled data from various sources: the Global Biodiversity Information Facility (GBIF; https://www.gbif.org), the Ocean Biogeographic Information System (OBIS; https://www.obis.org), the data from Villar et al.41, and the MAREDAT initiative42. Righetti et al.40 gathered >106 presences from nearly 1300 species sampled through various methodologies within the monthly climatological mixed-layer depth, at an average depth of 5.41 ± 6.95 m (mean ± sd), between 1800 and 2015. The species names were corrected and harmonized following the reference list of Algaebase (http://www.algaebase.org/) and were further validated by expert opinion. The final species list spanned most of the extant phytoplankton taxa composing the biodiversity of the euphotic zone of the global ocean. Fossil records, sedimentary records, and occurrences associated with senseless metadata were removed. This dataset has been mined to effectively obtain phytoplankton SR estimates for the global open ocean that were: (i) robust to sampling spatial-temporal biases11, and (ii) validated against independent data11.
Zooplankton occurrences
A new dataset of global zooplankton species occurrences was compiled in a comparable fashion to that put together for phytoplankton. Prior to retrieving the occurrence data online, we first identified the phyla (Order/Class/Family) that comprise the bulk of extant oceanic zooplankton communities: Copelata (i.e., appendicularians), Ctenophora, Cubozoa (i.e., box jellyfish), Euphausiidae (i.e., krill), Foraminifera, Gymnosomata (i.e., sea angels, pteropods), Hydrozoa (i.e. jellyfish), Hyperiidea (i.e., amphipods), Myodocopina (i.e., ostracods), Mysidae (i.e., small pelagic shrimps resembling krill), Neocopepoda, Podonidae and Penilia avirostris (i.e., cladocerans), Sagittoidea (i.e., chaetognaths), Scyphozoa (i.e., jellyfish), Thaliacea (i.e., salps, doliolids and pyrosomes), Thecosomata (i.e., sea snails, pteropods), and four families of pelagic Polychaeta (i.e., worms) that are often found in the zooplankton and whose species are known to display holoplanktonic lifecycles (Tomopteridae, Alciopidae, Lopadorrhynchidae, Typhloscolecidae). The presence data associated with species belonging to these groups were retrieved from OBIS and GBIF between the 12/04/2018 and the 18/04/2018 using online queries via the R packages RPostgreSQL, robis and rgbif. Since the Neocopepoda infra-class comprise several thousands of benthic and parasitic taxa43, a preliminary selection of the non-parasitic planktonic species had to be carried out prior to the downloading using the species list of Razouls et al.43 as a reference. The spatial distributions of the groups cited above were first inspected using GBIF’s and OBIS’s online mapping tools to evaluate the potential number of overlapping observations between the two databases. As a result of their relatively low contributions to total observations/diversity, and very high overlap between databases, the occurrences of Cladocera and Polychaeta were retrieved from OBIS only (which usually harbors more occurrences). On top of the data collected from OBIS and GBIF, the copepod occurrences from Cornils et al.44 and the pteropod occurrences from the MAREDAT initiative45 were added to the dataset. This initial collection of zooplankton observations gathered 4,899,151 occurrences worldwide.
Then, similar criteria as Righetti et al.11,40 were applied to progressively remove those presences that would be discordant with estimates of contemporary open ocean zooplankton diversity. The number of observations and species discarded after each main step and for each initial dataset are reported in Supplementary Data 1. We discarded records that: (i) presented at least one missing spatial coordinate, (ii) were associated with an incomplete sampling date (d/m/y), (iii) were associated with a year of collection older than 1800, (iv) were not associated with any sampling depth, (v) were not identified down to the species level, and (vi) were issued from drilling holes or sediment core data. For step (vi), a list of keywords (Supplementary Note 3) was used to identify the names of those original datasets that contained either fossil or sedimentary records. These first steps resulted in the removal of 1,766,783 occurrences (~36%). Like for phytoplankton, the remaining occurrences were associated with surface salinity values from the World Ocean Atlas (WOA) 201346 and bathymetry levels from the National Oceanic and Atmospheric Administration (NOAA) using the marmap R package. Occurrences associated with salinity levels 500 m. The average depth was used when maximal depth was not provided in the metadata. Therefore, the maximal depth of a zooplankton species occurrence allowed in our dataset is 500 m. This way, we tried to account for the zooplankton community that frequently performs diel vertical migration across the euphotic zone or the mixed layer, and that often co-occurs with species inhabiting surface layers. This removed 109,582 (~6%) occurrences. Next, for each phylum, OBIS and GBIF datasets were merged and the list of species names were extracted. Every species name was then carefully examined and compared to the taxonomic reference list of the World Register of Marine Species (WoRMS; http://www.marinespecies.org) for all taxa but copepods, for which we used the taxonomic reference of Razouls et al.43. This way, we rigorously harmonized and corrected the species names across all datasets. In addition, we used the notes and attributes of WoRMS to identify whether species were holoplanktonic or meroplanktonic (i.e., those species that have at least one benthic phase in their life cycle). Jellyfish species usually display a fixed polyp phase during their life cycle, therefore we used the dataset of Gibbons et al.47 to remove the species that were not holoplanktonic. Overall, these steps discarded 37,234 occurrences (~2%). One of two duplicate occurrences were removed from the dataset if they displayed the same species name, sampling depth, sampling date, and if they occurred within the same 0.25° x 0.25° cell grid. This last step removed 900,446 occurrences (~54%), highlighting the high overlap between the two main data sources. The remaining 766,033 presences were binned into the monthly 1° x 1° grid cell of the WOA to match the spatial resolution of the environmental predictors. The average maximum (±std) sampling depth was 73 ± 109 m and the average sampling year was 1985 ± 21. Observation densities were spatially biased towards the North Atlantic Ocean and the Southern Ocean (Supplementary Figure 1). The data reflected the historical seasonal sampling bias towards spring and summer. In the northern hemisphere, observations were equally distributed from March to October but constituted 78% of the data. In the southern hemisphere, 75% of occurrences were sampled between November and March. The final dataset gathered occurrences for 2034 different species (576 genera, 161 families) spanning all the major zooplankton phyla and several size classes. The only notable missing taxa are those belonging to the Cercozoa and Radiozoa because they present little to no species-level observations in online biodiversity repositories as they have been historically overlooked by traditional sampling techniques48.
Contemporary environmental conditionsA comprehensive set of environmental variables that are known to affect the physiology and constrain the distribution of plankton was prepared to define the candidate predictors for the SDMs. The R packages mainly used were raster and ncdf4. First, twelve primary variables that are relevant for modeling the distribution of both phytoplankton and zooplankton taxa were identified49,50,51,52. These were then aggregated into twelve monthly climatologies at a 1° x 1° resolution (i.e., the spatial cell grid of the WOA). The first six primary variables were sea surface temperature (SST, °C), sea surface salinity (SSS), nitrates (NO3-), phosphates (PO43-) and silicic acid (Si(OH)4) surface concentrations (µM), as well as dissolved oxygen concentration (dO2, ml l−1). Oxygen limitations and oxygen minimum zones (OMZ) are key factors controlling the horizontal and vertical distribution of zooplankton53,54. However, the effects of oxygen are often confounded with those of temperature because surface oxygen scales linearly with SST on a global scale. Therefore, dO2 at 175 m depth was used instead of surface dO2. For all six variables, the twelve monthly climatologies of the WOA13v2 (https://www.nodc.noaa.gov/OC5/woa13/woa13data.html) were used. In addition, satellite observations stemming from the Sea-viewing Wide Field-of-view Sensor (SeaWiFS; https://oceancolor.gsfc.nasa.gov/) over the 1997 to 2007 time period were used to derive monthly climatologies of Photosynthetically Active Radiation (PAR; µmol m−2 s−1) and chlorophyll (Chl; mg m−3), the latter serving as a proxy for surface phytoplankton biomass. Monthly climatologies of mixed-layer depth (MLD, m) based on the temperature criterion of55 from the Argo floats data (http://mixedlayer.ucsd.edu/) were also considered. Climatologies of surface wind stress (m s−1) were obtained from the Cross-Calibrated Multi-Platform56, using data from 1987 to 2011 (https://podaac.jpl.nasa.gov/). Climatologies of surface carbon dioxide partial pressure (pCO2; atm) were obtained from the Surface Ocean CO2 Atlas (SOCATv2; https://www.socat.info/) and made available by Landschützer et al.57. Lastly, a variable depicting sub-mesoscale dynamics and the strength of sea currents was derived from the daily satellite altimetry observations over the 1993–2012 period (https://cds.climate.copernicus.eu/#!/home): mean Eddy Kinetic Energy (EKE, m2 s−2). EKE was computed from the northward and eastward components of surface geostrophic seawater velocity (assuming the sea level as geoid), following the method of Qiu & Chen58. Such variable enabled us to account for the potentially important role of sub-mesoscale activity in structuring plankton biodiversity59.Then, nine secondary predictors were derived from some of the predictors described above. PAR over the MLD (MLPAR, µmol m−2 s−1) was calculated following Brun et al.49 An estimate of annual range of SST (dSST) was added by computing the difference between the warmest and the coldest temperature across the 12 months. The excess of nitrate to phosphate (N*, µM) relative to the Redfield ratio was computed as [NO3-] −16[PO43-]. Changes in N* represent varying conditions of denitrification and remineralization from N2-fixing organisms7. The excess of silicates to nitrates (Si*=[Si(OH)4]-[NO3-], µM) was also computed to represent regions where silicates are in excess compared to what diatoms would need to use up the nitrates7. Si* > 0 are indicative of conditions where diatoms can grow healthy. Since the distribution of macronutrients concentrations, chlorophyll concentration, and EKE values were all skewed towards lower values, we considered their logarithmic values (logNO3, logPO4, logSiOH4, logEKE, and logChl), based on either natural log or base 10, as additional predictors because they were much closer to a normal distribution.Species distribution modelingSDMs refer to a wide range of statistical algorithms that link an observed biological response variable (i.e., presence-only, presence/absence, abundance) to contextual environmental variables in the form of a response curve60. The latter is used to explore how a species’ environmental niche is realized in space and time24. In short, SDMs mainly rely on the following assumptions: (i) 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 scales61,62, which enables plankton species to display very large spatial ranges;11,47,60 (ii) species distributions are primarily shaped by the combinations of environmental factors that define the conditions allowing a species to develop. The latter assumption has been supported on macroecological scales, where the imprint of biological interactions (and dispersal) has been found to be relatively small63,64. Neither comparable abundance data nor presence/absence data were available from our datasets. In addition, presence-only data are less sensitive to discrepancies in species detection across various plankton sampling techniques. Therefore, based on species presences, we developed an exhaustive SDM framework to estimate plankton diversity patterns from an ensemble modeling approach65 that addresses the underlying main sources of uncertainties66,67.We follow the methodology developed in ref. 11, but simplify this approach to accommodate the limited predictor availability in the future model projections (see section “Choice of environmental predictors”), and the large number of diverse species we model in this work (sections “Background data (pseudo-absences)”, “Choice of environmental predictors”, and “SDMs evaluation and projections of monthly plankton species community composition”). We further derive SR based on habitat suitability rather than observed presence–absence data (section “Ensemble projection of global plankton species richness”). All methodological choices led to a minimization of computational cost and model complexity, while preserving all crucial patterns reported in ref. 11. Each methodological choice was carefully evaluated against other options, see sections below.
Background data (pseudo-absences)
Since we aimed at training correlative SDMs to model species distributions from presence data and environmental predictors, background data had to be simulated to indicate those conditions where a species is likely not to occur (i.e., pseudo-absences68). The generation of background data is a critical step in niche modeling experiments, and though no single optimal method has been identified by the niche modeling community, this step must address the important spatial and temporal sampling biases inherent to field-based observations. To do so, we made use of the target-group approach of Philipps et al.69, which has been shown to efficiently model phytoplankton distributions11. This method was found appropriate for our study because it generates background data according to the density distribution of the presence data, and therefore it: (i) does not induce additional bias to the initial biases in the presence data; and (ii) does not misclassify regions lacking observations (e.g., South Atlantic and Subtropical Pacific, Supplementary Fig. 1) as regions of absences.
For phytoplankton, we followed the background selection procedure described in Righetti et al.11. The authors used either the total pool of occurrences as a target-group, or defined three target groups based on the taxonomic groups contributing most to species diversity and observations (Bacillariophyta, Dinoflagellata, and Haptophyta). Background data of each species were randomly drawn based on the monthly resolved 1° x 1° occurrences of both their corresponding target groups, after applying an environmental stratification based on the SST and MLD gradients. This way, a species’ background is located at the sites where its lack of presence is most likely to reflect an actual absence. For each species, ten times more background data than presences were generated following the guidelines of Barbet-Massin et al.68. The amount of background data sampled from a specific SSTxMLD stratum was proportional to the number of monthly 1° x 1° cells provided by the target-group in this very stratum, thereby reflecting original sampling efforts. Both the total target-group background data (drawn from all sampling sites together) and the group-specific target background data were considered for our study, but both led to comparable estimates of phytoplankton species diversity (but see Fig. S3B of ref. 11).
The same method was applied for zooplankton species. First, we defined target groups based on their sampling distribution and broad taxonomic classification: Arthropoda (mainly copepods, but also krill and amphipods, that are sampled through similar techniques), Pteropoda, Chaetognatha, Cnidaria, Ctenophora, Chordata, Foraminifera, and Annelida. Unfortunately, the last three target groups displayed too few occurrences for drawing ten times more background data than presences. Consequently, their background data were drawn from the total pool of occurrences. Ctenophora also showed very few observations so they were merged with the Cnidaria as they are often considered together as jellyfish and collected in similar ways.
Total and target-group background data were drawn for all zooplankton species presenting more than 100 presences to run preliminary SDMs based on a preliminary set of predictors (Supplementary Fig. 2). These SDMs were then used to project preliminary diversity patterns for the four months that represent each season in the northern hemisphere (April, July, October, and January). These projections were examined for every group, and the predictive skills of the SDMs were evaluated using a repeated ten times split-sample test (see below). These tests showed that the total target-group background and the group target-group background converged towards SDMs of comparable skills and similar diversity patterns, except for the Chaetognatha, for which the target-group background leads to models of much poorer predictive skills (Supplementary Fig. 2). Furthermore, both background choices led to very similar latitudinal diversity gradients for phytoplankton (Fig. S3B in ref. 11). Therefore, to generate diversity patterns that are robust and consistent across the two trophic levels, the total target background data were used as standard background. The phytoplankton diversity pattern obtained with the total target background approach was only slightly lower in the Indo-Pacific and at very high latitudes11. Once they were generated, all background data were matched with monthly values for the 21 environmental variables described above.
Choice of the SDMs algorithms
The choice of the statistical method is a main source of uncertainty when projecting biodiversity scenarios through niche modeling66,67. Therefore, an ensemble forecasting strategy was adopted based on four types of SDMs that cover the range of algorithms types and model complexity that are commonly used:67,69 Generalized Linear Models (GLM), Generalized Additive Models (GAM), Random Forest (RF), and Artificial Neural Networks (ANN). The level of complexity of those models was constrained to avoid model overfitting70, a common pitfall when dealing with noisy and spatially biased data. SDMs including numerous predictors and parameterization features are more likely to fit spurious relationships and to be less transferable70,71. Consequently, the number of predictors was limited relative to the number of presences (see below) and the SDMs were tuned to fit relatively simple response curves. The GLM followed a binomial logit link, including linear and quadratic terms, and a stepwise bi-directional predictor selection procedure. The GAM also followed a binomial logit link. Smoothing terms with five dimensions, estimated by penalized regression splines without penalization to zero for single variables, were applied. Interaction levels between environmental predictors were set to zero for both GLM and GAM. The RF included 750 trees, and terminal node size was fixed at 10 to avoid having single occurrences as end members of some trees. The number of variables randomly sampled as candidates at each tree split (mtry parameter) was equal to the number of predictors used divided by three. The numbers of units in the hidden layers of the ANN, as well as the decay parameter, were optimized through five different cross-validations and a maximum of 200 iterations. Background data were weighted inverse-proportional to that of presence data (total weight = 1).
Choice of environmental predictors
To select for parsimonious and ecologically relevant sets of environmental predictors, a three-stage hierarchical selection framework was developed: (i) the distribution of the predictors’ values fitted to the presences were compared to their realized distribution between the main ocean basins to check whether one predictor could bias SDMs outputs towards a particular basin; (ii) pair-wise rank correlations between variables were examined, and one of two collinear variables was discarded where necessary; and (iii) models were trained to evaluate the explanatory power of several predictors sets of increasing parsimony, and rank the predictors within those sets at species-level. This selection procedure was carried out by separating phytoplankton from zooplankton since: (i) the two groups show different sampling distributions, and (ii) their niche dimensions might differ because of differences in their lifecycles (few days for phytoplankton, months to years for zooplankton) and biological requirements (photoautotrophy vs. heterotrophy and respiration). For those tests, only well-observed species with >100 occurrences were selected (nphytoplankton = 328; nzooplankton = 372). Ultimately, to account for the uncertainty in predictors choice, several final sets of predictors were defined based on the steps of the selection framework, and ensemble forecasting was adopted again (i.e., diversity estimates will be averaged across the sets of predictors).
Removal of variables impacted by sampling imbalances across ocean basins: Imbalance of sampling effort in geographical space can lead to sampling imbalance in environmental space if portions of an environmental gradient are strongly connected to an ocean basin that has been surveyed more extensively than others. To avoid such issues, the distributions of the annual values of the predictors were examined between the main basins (Arctic, Southern, Pacific, Indian and Atlantic Oceans). The most spatially imbalanced predictors were SSS and pCO2: the former is on average higher in the Atlantic Ocean, while the latter exhibits many of its most extreme values in the Peruvian upwelling system (Supplementary Note 3). The Peruvian upwelling is a hotspot of phytoplankton observations with clearly skewed observations (and the number of species sampled) towards pCO2 values > 400 atm (Supplementary Note 3). Plus, the pCO2 data do not cover the Arctic Ocean, the Mediterranean Sea, and the Red Sea. Consequently, pCO2 was discarded from the list of predictors to avoid strong sampling bias effects on SDM projections. A majority of the zooplankton data are concentrated in the Atlantic Ocean (Supplementary Fig. 1). As a result, the distribution of SSS values fitted to zooplankton occurrences is skewed towards SSS values >35 (Supplementary Note 3). As SSS is commonly used as a predictor for modeling the distribution of zooplankton47,48,49, we wanted to further examine its potential to act as a basin indicator rather than a predictor meant to represent an actual environmental control on species distribution. To do so, we performed ensemble SDMs projections for the zooplankton species, based on three variables sets: (i) without SSS, (ii) with SSS, and (iii) with Longitude (0°–360°) instead of SSS. Variables sets (ii) and (iii) led to very similar global zooplankton SR patterns with hotspots in the Atlantic Ocean. On the contrary, (i) led to more balanced zooplankton SR between basins without significantly lowering SDMs skills (Supplementary Note 3). We interpreted this as a bias in environmental space towards the conditions prevailing in the Atlantic Ocean, therefore we chose to discard SSS from the list of predictors.
Removal of collinear variables: Strong correlations among predictors can mislead the ranking of variable importance in SDMs72, so it has become common practice to exclude one of two variables that are highly collinear. Pair-wise Spearman’s rank correlation coefficients (⍴) were computed based on the predictors’ values fitted to the presences. When two variables exhibited a |⍴| > 0.70, the one displaying the distribution closest to a normal distribution was kept. From phytoplankton occurrences, we identified two clusters of strongly correlated variables: one comprising MLD, PAR, MLPAR (by construction), and wind stress (but PAR and MLD were only correlated at ⍴ = −0.66); and the other one comprising [NO3-], [PO43-] and their logged values. Similar clusters were found from zooplankton data, except that PAR was slightly less correlated to Wind stress (⍴ = −0.66) and MLD (⍴ = −0.58), and that [NO3-], [PO43-], plus their logged versions, showed stronger correlations with SST (⍴ = −0.80). As [NO3-] is a key factor for structuring planktonic systems7, and because we aimed to keep the variables sets as consistent as possible between species, logNO3 was kept as a candidate predictor. The variables retained for phytoplankton were: SST, dSST, logEKE, Si*, N*, logSiOH4, logNO3, logChl, wind stress, PAR, MLPAR, and MLD. The last four variables were kept to explore the outcome from alternative choices in the variables sets (but see below). The variables retained for zooplankton were the same but with the addition of dO2.
Examination of the explanatory power of predictors sets and ranking of predictors: To further evaluate which subset of these variable subsets are key to model species distributions, GLM and RF were performed for each species for several sets of decreasing complexity (from ten to five predictors), and the adjusted R2 of the models, as well as the ranking of predictors within each set, was extracted (Supplementary Note 1). GLM and RF were used here because they are part of the SDMs that will be used for projections afterwards and because they represent maximally different inherent model complexities among the SDM types used73. For GLM, predictor importance was determined according to their absolute t-statistic using the caret R package. For RF, predictor importance was based on the Gini index, which measures the mean decrease in node impurity by summing over the number of splits (across all trees) that includes a variable, proportionally to the number of samples it splits. The ranger R package was used for assessing variable importance with RF models. To keep the variables ranks comparable across predictors sets, rank values were normalized to their maximum. For each model type, the distributions of the models’ R2 and the distribution of the predictors’ ranks were examined for phytoplankton and zooplankton separately. The same was done between the main groups constituting the phytoplankton (Bacillariophyta, Dinoflagellata and Haptophyta) and the zooplankton (Copepoda, Chaetognatha, Pteropoda, Malacostraca, Jellyfish, Chordata and Foraminifera). This allowed us to identify the most important predictors for modeling the species distributions and to evaluate if a decrease in the models’ skill was linked to the removal of certain variables. Group patterns allowed us to test whether different groups differed in their main environmental drivers.
For phytoplankton species, 14 sets of variables were examined (Supplementary Note 1). The first nine aimed to test: (i) the impact of alternative choices between variables that were identified as collinear (wind vs. MLPAR vs. MLD + PAR); (ii) the impact of progressively discarding variables that initially presented lower ranks (logEKE, Si*), and (iii) the impact of choosing logNO3 over logSiOH4, two variables representing global macronutrients availability and that present relatively high correlation coefficient (⍴ = 0.59). The last five sets of predictors (10–14) aimed to test the impact of alternatively removing those variables that presented relatively high ranks in the previous sets: SST, dSST, N*, logSiOH4, logChl, PAR. In a similar fashion, 15 sets of variables were tested for zooplankton (Supplementary Note 1). The first ten aimed to test: (i) the impact of choosing wind stress over MLPAR or over MLD + PAR; (ii) the impact of selecting PAR over MLD; (iii) the impact of discarding Si*, N*, logEKE; and (iv) the impact of choosing logNO3 over logSiOH4 (⍴ = 0.64). The last five sets of predictors (11–15) aimed to test the impact of alternatively discarding the top five predictors: SST, dSST, dO2, logSiOH4, and logChl.
GLM and RF converged towards similar median variable rankings and evidenced high inter-species variability (Supplementary Note 1). For total phytoplankton, GLM identified the following median ranking across all species: SST > N* >logChl > logSiOH4 and dSST > logNO3 > logEKE > Si* > the PAR/MLD/MLPAR/wind stress cluster. RF ranked predictors in the following median order: SST and N* >logChl > dSST > logSiOH4 > logEKE and logNO3 > Si* > the PAR/MLD/MLPAR/wind stress cluster. Yet, both GLM and RF also identified PAR as a major predictor for Haptophyta, which does not appear in the rankings for total phytoplankton because Haptophyta represented ~9% of species composition only. Since adding PAR does not alter the models’ R2 for the Bacillariophyta and Dinoflagellata, it was retained for the final predictors sets. For total zooplankton, GLM ranked predictors in the following median order across all species: SST > dSST and logSiOH4 > logChl and logEKE > dO2 and logNO3 > N* >Si* and MLPAR > wind stress > PAR and MLD. RF identified the following median ranks: SST > dSST and dO2 > logNO3 > logSiOH4 > logChl and logEKE > N* and Si* > PAR > wind stress > MLD and MLPAR. These median rankings reflected those of the Copepoda since they represented >70% of all zooplankton species. Again, rankings displayed high variance, reflecting high inter-species variability. Overall, based on all the results shown above, eight different final predictors sets were kept for modeling the distribution of phytoplankton (n = 4) and zooplankton (n = 4) consistently. In contrast to ref. 11, predictor ensembles were defined across all species rather than for each species. This was due to multiple reasons: (i) predictor availability for future model projections was limited and did not allow for species-specific variable choices, (ii) computational constraints with regard to the total number of ensemble members that could be projected, (iii) the five sets already contain those predictors that explain a majority of the variability in most models, (iv) recent findings from Righetti et al. (in prep.) that the uncertainty due to predictor choice is low for models with optimized background selection.
Phytoplankton:
1.
SST, dSST, logChl, N*, PAR, and logNO3
2.
SST, dSST, logChl, N*, PAR, and logSiOH4
3.
SST, dSST, logChl, N*, PAR, logNO3 and Si*
4.
SST, dSST, logChl, PAR, and logNO3
Zooplankton:
1.
SST, dSST, dO2, logChl, and logNO3
2.
SST, dSST, dO2, logChl, and logSiOH4
3.
SST, dSST, dO2, logChl, logSiOH4, and N*
4.
SST, dSST, dO2, logChl, logNO3, and Si*
SDMs evaluation and projections of monthly plankton species community composition
Only species with more than 75 presences were considered for modeling plankton species distributions (nphytoplankton = 348; nzooplankton = 541) because we aimed to achieve a relatively high presence-to-predictors ratio (~15, which is the ratio achieved for a species with 75 presences and five to six predictors) to be more conservative than Righetti et al.11 (i.e., minimum 24 presences) since we aimed to project the SDMs in future conditions based on a pool of species for which we have high confidence. This is in line with Guisan et al.60 who suggest to maintain at least a ratio of ten. For each species, each SDM, and each set of predictors, presences and background data were randomly split into a training set (80%) and a testing set (20%) and these evaluation tests were repeated ten times. Therefore, 160 (four SDM types x four predictor sets x ten separate evaluation runs) models were trained per species, resulting in a total of 142,240 SDMs. Model skill was evaluated based on two widely used metrics: the True Skills Statistic (TSS74) and the Area Under the Curve (AUC60). TSS values range between −1 and 1, with null values indicating that models perform no better than at random. AUC ranges between 0 and 1, with values 0.30 were retained for the final ensemble projections (Supplementary Fig. 3).
In total, 860 species were considered as successfully modeled (nphytoplankton = 336; nzooplankton = 524; Supplementary Data 2). For those, each of the 160 SDMs was projected onto the twelve monthly climatologies of its corresponding predictors set and the projections were averaged over the ten cross-evaluation runs. This way, we obtained global maps of monthly mean presence probability for each of the 16 SDM x predictor set combinations. These maps are to be interpreted as habitat suitability patterns that highlight the regions of the global ocean where the environmental conditions are most favorable for a species to develop. Habitat suitability maps were not converted to binary presence–absence maps as probabilistic outputs provide more gradual responses that should better reflect the very dynamic occupancy patterns of plankton and that are better suited than threshold approaches for our purposes75,76,77. For each grid cell of the global ocean, the probabilistic estimates of species habitat suitability were stacked to obtain monthly estimates of species composition. All SDMs were trained, evaluated, and projected using the biomod2 R package.
Ensemble projection of global plankton species richness
For every SDM x predictor set combination and every month, we summed the species habitat suitability to estimate monthly SR. Then, annual average SR was estimated for each cell. The annual average was preferred over the annual integral because of the high latitudes that presented a lot of missing values in winter because of the lower coverage of satellite products. Since this diversity estimate is the sum of habitat suitability indices, it is to be interpreted as the amount of SR that the monthly/annual average environmental conditions should be able to sustain (i.e., potential SR). This way we obtained 16 estimates of annual SR for total plankton (all 860 species together), phytoplankton and zooplankton. Ensemble projections of annual SR were then obtained for these three categories by averaging the annual SR estimates.
As the biological data used to train the SDMs span several decades (mostly between the 1970s and 2000s), our diversity estimates are integrative of changes in SR and species composition (i.e., changes in beta diversity) that occurred during these decades. The phytoplankton species modeled are mainly members of the Bacillariophyceae (45.8%), and the Dinoflagellata (45.8%), which usually rank among the large marine microalgae. Therefore, the phytoplankton SR estimates shown here should be mainly representative of the microphytoplankton (20–200 µm) rather than smaller size fractions. Nearly half (51.9%) of the zooplankton species modeled are Copepoda, making it the most represented groups in the zooplankton SR patterns followed by: Malacostraca (crustacean macrozooplankton such as Euphausiids and Amphipods; 13.9%), Jellyfish (13.1%), Foraminifera (5%), Chaetognatha (5%), Pteropods (4%), Chordata (3%). The last 4% of species modeled are a mix of Annelids and Branchiopods. The full list of the species modeled as well as their taxonomic classification is given in the Supplementary Data 2.
We underline that the 860 species modeled are those for which we have enough records for training reliable SDMs. Therefore, these species are likely to be those that are the most frequently detected by conventional sampling and identification techniques, either because: (i) they are the ones dominating total plankton abundance, which makes their collection more likely; or (ii) they are larger species (in terms of cell volume or body length), which would facilitate their sampling and identification under the binocular or recent imaging systems. We acknowledge that our approach does not allow us to account for rare taxa and thus under samples the true diversity of the marine plankton. Nonetheless, we argue that our approach does allow us to estimate global plankton diversity patterns as the species dominating plankton abundances are those carrying biogeographical information30, meaning their distribution and abundance patterns can be correlated to environmental gradients. Meanwhile, the patterns of rare and non-dominant species, which constitute the majority local SR, exhibit no biogeographical signature30. This has been supported in the previous study of Righetti et al.11 where the authors showed that global SR patterns were robust to the progressive exclusion of taxa with relatively few records.
We also acknowledge our estimates of species distribution might be biased by imbalances in species detection and sampling effort between sampling cruises, as those rely on a wide range of collection and identification methodologies. We argue that such biases are particularly significant when relying on abundance data, and that we mitigate them by: (i) converting all observations to presences and aggregating them onto a 1° x 1° grid; (ii) modeling SR as an emergent property overlapping the distribution of single species with equal weighting rather than modeling SR directly, in which case diversity estimates would be highly sensitive to sampling effort imbalances; (iii) by designing the SDMs in a way that accounts for spatial and temporal sampling biases in geographical and environmental space; and (iv) by tuning down the complexity of the SDMs (i.e., reduced number of features and predictors) in order to avoid model overfitting70.
Future environmental conditions in the global surface oceanThe future monthly fields of the selected environmental predictors were obtained from the projections for the 2012–2100 period of five ESM simulations for the IPCC’s RCP8.5 scenario from the MARine Ecosystem Model Intercomparison Project (MAREMIP, http://pft.ees.hokudai.ac.jp/maremip/index.shtml78) and/or the Coupled Model Intercomparison Project 5 (CMIP579). The model ensemble contained the following five ESMs (with their embedded ocean and ecosystem models indicated indicated in brackets after the semicolons): Community Earth System Model version 1 (CESM1, POP-BEC), Geophysical Fluid Dynamics Laboratory Earth System Model with Modular Ocean Model version 4 (GFDL-ESM2M; MOM-TOPAZ), Institut Pierre Simon Laplace Climate Model version 5A-LR (IPSL-CM5A-LR; NEMO-PISCES), Centre National de Recherches Météorologiques Climate Model version 5 (CNRM-CM5; NEMO-PISCES) and the Model for Interdisciplinary Research on Climate version 5 (MIROC5; MRI.COM-MEM). All ESMs were fully-coupled except for MIROC5 for which the ocean model was forced by the atmospheric component. All of the projections used here were benchmarked, quality-controlled and described in the previous multi-model comparison studies of Laufkötter et al.26,80. Considering the scope of the present study, we refer to these authors’ previous extensive descriptions for the full detail of the ESMs used here. Taken together, the present five ESM ensemble gathers models of various sensitivity to future climate forcing, and thus provides a wide range of alternative environmental conditions projected for the future surface ocean. With the present ESM ensemble, we account for the variability in the choice of the climate model, which is known to be a significant source of uncertainty in biodiversity projections; this source being consistently lower than those associated with SDM choice, though28,66,67.The monthly projections of the five selected ESMs were interpolated on the 1° x 1° cell grid of the WOA (i.e., the one used to train our SDMs) over the 2012–2100 period for all the nine chosen environmental predictors. To obtain future monthly climatologies that span a comparable amount of temporal variability as the in situ climatologies used to train the SDMs (~20 years), a baseline and an end-of-century time periods were first defined (2012–2031 and 2081–2100, respectively) for every ESM projection run. The 12 monthly climatologies were derived based on the models’ monthly projections and monthly anomalies were computed by subtracting the baseline values to the end-of-century ones. For dSST (i.e. annual range of SST), the annual maximum of SST was derived from the monthly climatologies and the difference between the baseline and the end-of-century dSST provided the delta value. These anomalies can be either positive or negative and they represent the difference in the predictors’ condition due to future climate change under the RCP8.5 GHG concentration scenario25. To obtain the final conditions prevailing in the surface ocean for the end-of-century period, the delta values were simply added to the in situ climatologies representing the conditions in the contemporary ocean. The SDMs of the 860 plankton species successfully modeled were then projected onto these future monthly climatologies for each of the ESM. This way, we estimate the monthly probability-based species composition in the future global ocean for each of the 80 combinations of SDMs (n = 4), ESMs (n = 5), and predictor set (n = 4). Overall, our ensemble forecast approach65 generates an unprecedented set of 825,600 species-level estimates of global future habitat suitability patterns. Finally, mean annual SR and community composition were calculated for total plankton, phytoplankton and zooplankton for each of the 80 possible combinations of projections, as described in section “Ensemble projection of global plankton species richness”.Analyses
Ensemble projections of changes in species richness, community composition turnover, and changes in species associations between the contemporary and the future ocean
For each of the 80 projection combinations described above, the mean annual SR estimates for the contemporary ocean were subtracted to their corresponding mean annual SR estimates for the future ocean to compute the percentage difference in mean annual SR (%∆SR) for total plankton, phyto- and zooplankton. The %∆SR represents the emergent change in SR caused by future climate change(s) through changes in species-level habitat suitability patterns. While changes in SR indicate climate change impacts on plankton alpha diversity, these do not inform us on the potential impacts on beta diversity (i.e., changes in community composition81). A community that experiences the replacement of all its constituting species by an equivalent number of newcomers will display a 100% rate of community turnover but no changes in SR. To investigate the amplitude of global plankton species turnover triggered by climate change, we examined future total turnover in annual species composition using Jaccard’s dissimilarity index while decomposing its two additive components: nestedness (i.e., changes in SR) and true turnover (ST), which indicates the % of species that will be replaced in a community27 using the betapart R package.
To do so, the mean annual species habitat suitability patterns used to estimate the ensemble changes in SR had to be converted to presence–absence maps as the Jaccard’s dissimilarity index requires binary inputs80. A range of thresholds (0.10 to 0.80, by steps of 0.01) was first explored for each SDM type (GLM, GAM, ANN, and RF) to infer threshold-based annual SR patterns. Then, we quantified the similarity of the threshold-based annual SR vs. the probability-based annual SR using Spearman’s rank correlation coefficient (⍴) and ordinary linear regressions (R2) to identify the range of thresholds that best match the probability-based estimates. The 0.25–0.40 range provided the most similar global SR patterns for GLMs, GAMs and ANNs (all ⍴ > 0.95, and all R2 > 0.90). The 0.10–0.25 range was chosen for RF models. These ranges largely overlap with the species mean probability thresholds that maximize the TSS/AUC evaluation metrics, which are commonly used to convert habitat suitability into presence–absence maps. However, the maximizing-TSS approach tends to underestimate the natural gradual response of organisms to environmental variations, which is particularly problematic when dealing with SR patterns of widely-dispersed and climate-sensitive ectotherms such as the plankton75,76,77. Therefore, we chose to rely on a range of thresholds instead as it enables us to account for a wider range of possible realizations of community composition and better reflect the dynamic occupancy patterns inherent to planktonic taxa. Consequently, we derived ST estimates in annual species composition for each of the SDM-dependent threshold mentioned above and every of the 80 annual projections combinations, for total plankton, phyto- and zooplankton separately. Again, the ensemble projection in annual ST was derived by averaging those projections.
We further examined how climate change could impact not only community composition but also those species associations within the community that represent potential biotic interactions, which support ecosystem functioning3,20,21. Based on the mean annual species composition estimates used to compute ST rates, a text analysis algorithm82,83 was used to identify pairs of species, which co-occur more frequently than expected given their individual occurrence. In short, the text analysis algorithm assigns an association score to each possible pair of two plankton species in all grid cells based on a likelihood ratio (LLR)82. The latter compares the probability of two species co-occurring together to the probability of one species occurring without their partner (i.e., two alternative probabilities), or when both are projected as absent, based on a combination of Shannon’s entropy indices (H’). LLR values are >0 and they scale with the significance level of the projected species association, whether it is a positive (co-occurrence) or a negative (one-sided occurrence or co-absence) association. To disentangle between those two cases, when a species pair displayed an observed co-occurrence frequency lower than the product of the one-sided occurrence frequencies normalized to sample size, its LLR value was multiplied by −1. This way, we can identify those species pairs whose co-occurrence probability is lower than the products of the two one-sided occurrence probabilities (LLR More
