Major restructuring of marine plankton assemblages under global warming

The species diversity of marine plankton governs some of the most important marine ecosystem services. In the sun-lit layers of the oceans, photoautotrophic phytoplankton are responsible for about 50% of Earth’s annual net primary production and are grazed by heterotrophic zooplankton, sustaining global fisheries production. Together, these trophic levels drive the biological carbon pump, a key determinant of the ocean-atmosphere CO2 balance. Diversity of phyto- and zooplankton modulates this pump and influences fisheries recruitment. Most studies indicate plankton diversity is largely controlled by climate with temperature as the main driver; warm temperatures promote diversity by enhancing speciation and metabolic rates. However, ocean warming forces poleward range shifts and is expected to trigger species extirpations and community restructuring, with potentially deleterious consequences for food-web functioning and biogeochemical cycles. Yet, responses of plankton diversity and associated ecosystem functions to future warming remain poorly understood across clades and trophic levels. Earth system models with embedded ecosystem models lack trait and species resolution, and historical observations were too sparse to enable empirical global projections, leading past studies to rely on virtual taxa or limited observations. Consequently, the extent to which global phyto- and zooplankton diversity may be affected by future climate change remains unclear. We address these limitations by modeling the monthly and mean annual diversity patterns from distributions of 860 plankton species via an ensemble of species distribution models (SDMs) trained on nearly one million occurrences and matched with environmental climatologies. Assuming niche conservatism, we project habitat suitability into late-century conditions under RCP8.5 using five CMIP5 ESMs and derive species richness and turnover, assessing uncertainties across 80 ensemble members. Our ensemble projects marked increases in phytoplankton richness over most basins (except the Arctic) and strong zooplankton richness gains in temperate to subpolar latitudes with slight tropical declines, alongside substantial species turnover, implying threats to plankton-mediated ecosystem services.

Study scope and data: Compiled 934,696 species occurrence records for 860 open-ocean plankton species (336 phytoplankton; 524 zooplankton) spanning 13 phyla, 71 orders, and 324 genera. Occurrences aggregated monthly on a 1°×1° grid; regions with seafloor depth <200 m excluded. Phytoplankton data from GBIF, OBIS, MAREDAT, Villar et al.; zooplankton dataset newly compiled from OBIS, GBIF, plus group-specific sources (e.g., MAREDAT pteropods, Cornils et al. for copepods). Taxonomy harmonized (AlgaeBase, WoRMS; Razouls et al. for copepods); meroplankton excluded; duplicates removed; depth-limited to ≤500 m to represent euphotic/mixed layer communities. Environmental predictors: Prepared monthly climatologies (1°×1°) for candidate predictors: SST, annual SST range (dSST), dO2 at 175 m (zooplankton), macronutrients (NO3, PO4, Si(OH)4; log-transformed), chlorophyll a (logChl), photosynthetically active radiation (PAR), mixed-layer depth (MLD), MLPAR, wind stress, eddy kinetic energy (logEKE), and nutrient stoichiometric indices (N*=[NO3]-16[PO4], Si*=[Si]-[NO3]). Discarded SSS and pCO2 due to sampling/basin biases. Addressed collinearity via Spearman rank correlations, retaining parsimonious sets. SDMs and background: Presence-only SDMs with pseudo-absence/background using the target-group approach to match sampling bias; total target background adopted for consistency across taxa. Four SDM algorithms: GLM (binomial with linear/quadratic terms, stepwise), GAM (binomial, smoothers with k=5, no interactions), Random Forest (750 trees, min node size 10, mtry ≈ predictors/3), Artificial Neural Networks (tuned hidden units/decay via cross-validation). Predictor pools: four final sets for phytoplankton and four for zooplankton selected after hierarchical screening and variable ranking exercises. Only species with >75 presences modeled; model skill evaluated via repeated 80/20 splits (10 repeats) using TSS and AUC; retained species with mean TSS>0.30 (resulting in 860 species). Contemporary projections: For each species, SDM, and predictor set, projected monthly habitat suitability indices (HSI) and averaged across evaluation runs. Stacked species HSIs to derive monthly and mean annual species richness (SR) for total plankton, phytoplankton, and zooplankton. Future projections: Used five CMIP5/MAREMIP ESMs (CESM1, GFDL-ESM2M, IPSL-CM5A-LR, CNRM-CM5, MIROC5) forced by RCP8.5. Computed monthly climatologies for baseline (2012–2031) and end-of-century (2081–2100), derived anomalies, and added them to observation-based climatologies (delta method). Projected all SDMs into future conditions for each ESM and predictor set. Generated 80 ensemble members (4 SDMs × 4 predictor sets × 5 ESMs) and computed mean annual SR and ensemble statistics (mean, interquartile range). Turnover and associations: Computed beta diversity using Jaccard dissimilarity, decomposed into true species turnover (ST) and nestedness. Converted HSI to presence-absence using SDM-specific probability thresholds determined to best match probability-based SR (GLM/GAM/ANN: 0.25–0.40; RF: 0.10–0.25). Estimated contemporary vs. future ST for each ensemble member and averaged. Identified changes in species associations via a likelihood ratio-based text analysis of co-occurrence (LLR), retaining top positive associations (≥75th percentile; negative LLR treated as non-significant), and quantified gains/losses. Uncertainty and novel climates: Quantified ensemble variance and sensitivity; SDM sensitivity ranking GLM < GAM < ANN < RF; ESM sensitivity GFDL-ESM2M < CNRM-CM5≈CESM1 < IPSL-CM5A-LR < MIROC5. Identified novel environmental conditions using MESS diagnostics at species level monthly, summarized annually; SST increases in tropics primary driver of non-analog conditions. Link to ecosystem services and traits: Clustered global ocean into six regions using PCA of five impact metrics (%ΔSR for phyto and zoo; ST for phyto, zoo; total turnover), applying k-medoids on Euclidean distances (PC1–PC4 explain >95%). Defined a severity index from weighted PC scores and compared against proxies of ecosystem services: megafauna SR, small pelagic catches, NPP, export flux (FPOC), export efficiency (FPOC/NPP), and plankton size index. Assessed size-structure changes for diatoms (cell volume, S/V, C content) and copepods (body length) by HSI-weighted metrics for contemporary and future conditions.

• Contemporary diversity gradients: Strong latitudinal gradient with total plankton SR decreasing from equator to poles. Phytoplankton SR peaks near equator/tropical upwellings; zooplankton SR peaks in subtropics (20°–30°) with a tropical dip. Patterns across 10 functional groups (PFGs) match prior observations.
• End-of-century changes (2081–2100, RCP8.5): Global median increase in total plankton mean annual SR of 5% (IQR 1.9–9.6%; p<1e-15), driven by strong gains at 40°–55° latitudes (+22%; 17–62%; p<1e-15), offset by smaller tropical gains (+4%; 0–32%; p<1e-8).
• Phytoplankton SR: Global +16% (11–22%; p<1e-15); tropics +21% (12–100%; p<1e-15); temperate/subpolar +13% (10–38%; p<1e-15); Arctic north of 70°N modest decrease −11% (−17 to +24%; p<1e-9).
• Zooplankton SR: Near-zero global median change (+0.4%), masking regional contrasts: +24% (19–69%; p<1e-15) at 40°–55°; slight tropical decline −4% (−6 to +7%; p<1e-6); Southern Ocean −3% (−9 to +5%; p<1e-12).
• Temperature as main driver: Predictor dominance analyses identify SST as primary control of contemporary SR and future changes. Log(SR) vs. thermal energy slopes consistent with MTE: phytoplankton ~0.33 (eV)−1 above 22 °C; zooplankton ~0.66 (eV)−1 between 11–20 °C; non-linearities (e.g., zooplankton inflection ~25 °C) explain tropical zooplankton SR declines and high-latitude reductions.
• Range shifts: 79% of all species shift poleward; more common in zooplankton (87%) than phytoplankton (67%). Median poleward shift velocity 35 ± 22 km/decade (phyto 34 ± 28; zoo 36 ± 20), consistent with isotherm shifts.
• Species turnover: Global median true turnover (ST) 18% (±10) for total plankton; higher for zooplankton than phytoplankton (Kruskal–Wallis Chi2=2577; p<1e-3). Latitudinal gradient: tropics 16% (±7.9); >60° latitudes 45% (±16). Arctic communities exhibit ST >45% broadly, indicating borealization.
• Species associations: ~40% of species associations reshuffled. Phytoplankton: 28% gained, 10% lost; zooplankton: 10% gained, 26% lost. Median 27% of future phyto–zoo associations are potentially novel.
• Functional reshuffling and size: Divergent SR responses across PFGs (e.g., tropical declines for haptophytes, copepods, euphausiids, jellyfish, chordates). Projected replacement of larger high-latitude diatoms/copepods by smaller warm-water species, implying reduced export efficiency.
• Uncertainty: Largest from SDM choice; ESM choice secondary; predictor pool minor. Non-analog climates emerge mainly in tropical oceans due to SST increases; sensitive ensemble members predict large phytoplankton SR increases where monthly SST exceeds contemporary ranges.
• Ecosystem services overlap: Regions with most severe diversity and composition changes (temperate and Arctic) currently have high export efficiency and support small pelagic fisheries and megafauna, indicating elevated risk to ecosystem services.

The study demonstrates that warming-driven changes in ocean temperature will reorganize global plankton biodiversity, with substantial poleward shifts and increased species richness in temperate to subpolar regions but declines or weak gains in the tropics (especially for zooplankton). Temperature emerges as the dominant driver, with non-linear diversity–energy relationships explaining spatial heterogeneity in responses. High true turnover rates (up to ~40–45% at high latitudes) indicate major species replacement, beyond what richness changes alone reveal, implying profound restructuring of community composition and potential alterations of trophic interactions. Approximately 40% of species associations are reshuffled, and over a quarter of phyto–zoo associations become novel, suggesting significant rewiring of plankton interactomes. These biodiversity changes overlap with regions that provide key ecosystem services (carbon export, fisheries, megafauna biodiversity), implying risks to the biological carbon pump and fisheries yields. The projected shift towards smaller diatom and copepod sizes at high latitudes likely reduces export efficiency and may diminish carbon sequestration. While phytoplankton richness generally increases, zooplankton declines in the tropics raise concerns about food-web support in warm regions. The findings address the research question by quantifying, across trophic levels and functional groups, how climate change is likely to alter plankton diversity patterns and community structure globally, highlighting temperature’s central role and the potential ecosystem service consequences.

Using an ensemble of SDMs for 860 plankton species and RCP8.5 climate projections, the study provides global, trophically explicit projections of plankton biodiversity change. Key contributions include: (i) quantifying differential responses of phytoplankton and zooplankton richness, (ii) documenting widespread poleward range shifts at ~35 km/decade, (iii) revealing high species turnover and widespread reshuffling of species associations, and (iv) linking biodiversity changes to potential impacts on carbon export and fisheries through functional and size-structure shifts. The results suggest that climate change will substantially restructure plankton communities, particularly in temperate and polar regions that are critical to ecosystem services. Future work should improve monitoring of plankton diversity (including rare taxa), integrate abundance and trait data, better resolve biotic interactions, and develop models that couple biodiversity dynamics to ecosystem functioning under multiple stressors. Enhanced treatment of novel climatic conditions and uncertainties across SDMs/ESMs will further refine projections.

• Taxonomic and sampling biases: Focus on species with sufficient occurrence records likely overrepresents dominant and larger taxa, underrepresenting rare and small-sized taxa (e.g., Cercozoa, Radiolaria), potentially underestimating true diversity.
• Presence-only data: Lack of comparable abundance/presence–absence data limits inference on population dynamics and interaction strengths; habitat suitability not converted to strict occupancy.
• Niche conservatism: Projections assume stable niches over time, not accounting for adaptation or rapid evolution.
• Predictor limitations and biases: Discarded pCO2 and SSS due to sampling/basin biases; collinearity addressed but residual biases may persist. Oxygen considered at depth for zooplankton; other factors (e.g., micronutrients, pH/acidification) not explicitly included.
• Model uncertainty: Largest uncertainties from SDM choice; ESM choice secondary; emergence of non-analog climates (especially in tropics) can reduce forecast reliability.
• Spatial scope: Open ocean only (depth >200 m); coastal/shelf systems excluded.
• Biotic interactions: Associations inferred from co-occurrence; direction and mechanisms not resolved; food-web impacts not explicitly quantified.
• Temporal aggregation: Use of climatologies may smooth seasonal/episodic dynamics; baseline observations span decades.
• Conversion to binary for turnover: Thresholding HSIs introduces additional uncertainty in turnover estimates.