Marine Ecosystem Changepoints Spread Under Ocean Warming in an Earth System Model

Sudden shifts in marine plankton communities in response to environmental changes are of special concern because of their low predictability and high potential impacts on ocean ecosystems. We explored how anthropogenic climate change influences the spatial extent and frequency of changepoints in plankton populations by comparing the behavior of a plankton community in a coupled Earth system model under pre‐industrial, historical 20th century, and projected 21st century forcing. The ocean areas where surface ocean temperature, nutrient concentrations, and different plankton types exhibited changepoints expanded over time. In contrast, regional hotspots where changepoints occur frequently largely disappeared. Heterotrophy and larger organism sizes were associated with more changepoints. In the pre‐industrial and 20th century, plankton changepoints were associated with shifts in physical fronts, and more often with changepoints for iron and silicate than for nitrate and phosphate. In the 21st century, climate change disrupts these interannual‐variability‐driven changepoint patterns. Together, our results suggest anthropogenic climate change may drive less frequent but more widespread changepoints simultaneously affecting several components of pelagic food webs.

Rapid changes in marine plankton populations may occur when changes in environmental conditions generate a sudden forcing on the community (Ardyna et al., 2014;Muller-Karger et al., 2019), or when gradual changes in environmental conditions provoke a nonlinear response (Scheffer et al., 2001;Stock et al., 2014a). With or without environmental forcing, rapid changes can also emerge from the internal dynamics of plankton communities (Barton et al., 2020;Di Lorenzo & Ohman, 2013;Huisman & Weissing, 1999). Because of the challenges of collecting long-term and broad-scale measurements of plankton populations (Benway et al., 2019), much remains unknown about how frequent abrupt changes are in the lower trophic levels of pelagic marine ecosystems, how abrupt changes differ with trophic status or organism size (Barton et al., 2020), how abrupt changes in plankton populations correspond to those in physical or chemical environmental conditions, or how their frequency and distribution will change (or already are changing) with climate change (Beaulieu et al., 2016).
Here, we address these questions in the context of a plankton community model integrated into a global Earth system model (ESM). The ESM includes a comprehensive ecosystem model, called Carbon, Ocean Biogeochemistry and Lower Trophics (COBALT), that captures regional and seasonal variations in integrated ecosystem properties (e.g., chlorophyll and primary production) as well as the emergent biogeographies of phytoplankton and zooplankton across contrasting body sizes, functional groupings, and predator-prey interactions (Stock & Dunne, 2010;Stock et al., 2014b). We investigate the occurrence of abrupt transitions in marine ecosystems through the analysis of changepoints, a general term defined as a time point where a change in a statistical property of a time series can be identified, but here more specifically meaning interannual, decadal, or multidecadal changes in trends or mean values (Reeves et al., 2007). We focused on centennial pre-industrial, historical, and projected climate change simulations for surface ocean temperature, nitrate concentrations, and phytoplankton and zooplankton of different sizes. Our objective is to map where changepoints occur and how frequent they are, with respect to oceanographic features, trophic levels, and climatic forcing. Our focus is specifically on plankton changepoints, and while it is unwieldy to establish the mechanisms underlying all individual plankton changepoints, we identify some common environmental drivers.

Numerical Model
Our analyses use simulations conducted with ESM2M-COBALT (Stock et al., 2014a). ESM2M-COBALT was derived from GFDL's ESM2M ESM (Dunne et al., 2012 by replacing the ocean biogeochemical component with the Carbon, Ocean Biogeochemistry and Lower Trophics (COBALT, Stock et al., 2014b) model, while preserving other Earth system components. The ocean component is GFDL's Modular Ocean Model version 4.1 (MOM4p1; Griffies, 2009) with a horizontal resolution of 1° and 50 vertical layers. The comprehensive ocean biogeochemistry and ecosystem COBALT model (Stock & Dunne, 2010;Stock et al., 2014aStock et al., , 2014b includes 33 prognostic tracers, including three phytoplankton groups, three zooplankton groups, and tracers representing the coupled elemental cycles of carbon, nitrogen, phosphorus, silicate, and iron, as well as alkalinity and lithogenic material. As described in Stock et al. (2014a), ESM2M-COBALT simulations follow the protocols of Phase 5 of the Coupled Model Intercomparison Project (CMIP5; Flato et al., 2014) enlisted herein, and include (a) a 1,500 yr spin-up simulation with 1,860 radiative forcings and potential vegetation, the last 100 yr of which are used as a pre-industrial (hereafter PI) control simulation, (b) a historical simulation from 1860 to 2004 featuring observed greenhouse gas concentrations, solar insolation, volcanic eruptions, ozone, and land use changes, from which we take the years 1901-2000 for the historical 20th century simulation, and (c) a future projection from (b) out to 2100 under the RCP8.5 scenario, the highest emission scenario among the set of Representative Concentration Pathway scenarios (Riahi et al., 2011), to which we append the years 2001-2004 from the historical simulation and hereafter refer to as the 21st century simulation . The ecological and biogeochemical properties from COBALT simulations compare well with global observations over the past few decades (Stock & Dunne, 2010;Stock et al., 2014aStock et al., , 2014bStock et al., , 2017. Plankton traits and ecological interactions in the COBALT model are tied to body size (Stock et al., 2014a) and parametrized by allometric trait relationships gleaned from large compilations of laboratory measurements. For example, smaller model phytoplankton have higher nutrient affinity compared to the larger phytoplankton (K. F. Edwards et al., 2012), but are grazed upon by smaller, more efficient zooplankton (Hansen et al., 1997). Similar allometric trait formulations and trait-tradeoffs underpin other size-based plankton community models (Baird & Suthers, 2007;Banas, 2011;Taniguchi et al., 2014;Ward et al., 2012), and allow the COBALT model to plausibly simulate biogeographical and phenological patterns for a range of phytoplankton and zooplankton sizes (Stock et al., 2014b). Phytoplankton growth in the model is an increasing exponential function of temperature (Eppley, 1972), such that growth is higher at warmer temperatures and growth acclimates instantaneously to changes in temperature for each phytoplankton type. A similar Q 10 was applied to zooplankton growth and grazing, although remineralization and particle aggregation and export were independent of temperature (Stock et al., 2014b). Plankton therefore do not have discrete temperature niches ). In this study, we use an identical model formulation to Stock et al. (2014b). This type of size-structured plankton community model is ideal for studying how changepoints vary across body size and trophic level because it not only simulates these state variables but also because it encodes allometric gradients in physiological rates, generation times, and interaction strengths that are likely to influence the occurrence of model changepoints. For the purposes of our analyses, we exclude nitrogen-fixing phytoplankton (diazotrophs); while their presence in COBALT is essential for maintaining ocean surface nitrogen inventories, their traits and ecological interactions are not as well constrained as for other model phytoplankton (but see Monteiro et al., 2011) and their contribution to total productivity is modest. In our analyses, we use depth-integrated (over the top 100 m) plankton biomass model output but temperature and nutrient data from the surface layer. We acknowledge that some changes could occur in depth distributions of plankton over the centuries of simulation that we do not resolve, but prefer instead to focus on horizontal and temporal patterns. The surface data for environmental properties are sufficient as we are interested in broad temporal and spatial gradients in these properties, rather than vertical profiles. Surface nutrients are used in conjunction with 100 m integrated biomass because surface nutrient values are most reflective of nutrient limitation in the euphotic zone. The plankton community model in COBALT represents differences in body size and associated physiological traits and ecological interactions, as well as patterns of biogeography and phenology for model plankton types that are similar to observations (Stock et al., 2014b).

Changepoint Analysis
We then identify changepoints for the time series at each model grid box and in each simulation for each model variable of interest. In essence, changepoint methods are designed to identify points in a time series where a statistical property of that time series changes, such as its mean or its trend; such methods have been widely adopted by the statistics community for their robustness and ability to handle for example, changes in time series with trends (Killick et al., 2020;Reeves et al., 2007). As we are analyzing the annually averaged output of centennial simulations, changepoints will necessarily be shifts between two multi-annual to multi-decadal average states, such as a switch between alternate stable states centered on different mean values, switches between increasing and decreasing trends, or the emergence of a multidecadal trend from a statistically steady baseline. We use the EnvCpt changepoint package in R (Killick et al., 2020), which relies on the Pruned Exact Linear Time algorithm (Killick et al., 2012) and selects among statistical models with changepoints relative to each other. We considered eight statistical models: (a) a constant mean, (b) a linear trend, (c) constant means with changepoints, (d) linear trends with changepoints, (e-h) AR(1) autoregressive versions of 1-4. We thus in essence fit each time series for each grid point, variable, and century simulation with a piecewise combination of constant values or linear trends, and formally define a changepoint as a point in a given time series where the best-fitting model switches between different segments of this piecewise fit. Qualitatively this corresponds to points in time where approximately interannual or (multi-)decadal shifts in the baseline value or baseline trend are identifiable. Figure  S1 in Supporting Information S1 shows examples of each of these fit to COBALT time series, Figure 1 shows  Table  S1 in Supporting Information S1 shows a pedagogical example wherein the number of changepoints detected is robust to the superposition of a trend.
Both the presence and number of changepoints detected are robust across the models allowing abrupt changes (i.e., models 3, 4, 7, and 8 above detect the same number of changepoints for almost all grid cells for all three centuries). As we are interested in the quantitative characterization of the incidence and prevalence of changepoints, we group the statistical models into those with changepoints and those without, and hereafter focus on the absolute frequency of changepoints detected during each simulation, for each variable, and for each location. We also focus on this robust characteristic of the number of changepoints because the EnvCpt package does not require the piecewise linear models to be continuous; it also does not require changes to be of a certain amplitude for any statistical model. We discard all changepoints that occur in the first two or last 2 yr of any time series to avoid edge effects. Note that the changepoint method does not require ad hoc parameters to be chosen, such as minimum segment length or transition magnitude. To report global changepoint statistics, we weight each grid cell proportional to its area. Note that we analyze each simulation separately. While the sign and occurrence of changepoints is a robust feature of our analysis, the quantification of the amplitude of a given changepoint depends on which of the eight aforementioned models is selected; because which model is selected as best-fitting is sensitive to the choice of fitting metric, we do not consider changepoint amplitude and instead focus on robust aspects of our analysis. In other words, all models 1-8 are fit to each variable, grid cell, and century time series to select the best-fitting model according to the Akaike Information Criterion. Then we focused on the number of changepoints detected by the selected model, which is a robust metric that is almost independent of model choice (i.e., using the Bayesian Information Criterion yields the same number of changepoints almost everywhere but in different models). The overall incidence of changepoints may be more important than their specific magnitude, as recent work suggests that even small-amplitude ecological and biogeochemical changes can have appreciable consequences (Barton et al., 2020;Stock et al., 2014a).

Results and Discussion
The marine ecosystem model we employ exhibits a substantial increase in the prevalence (i.e., spatial extent) of changepoints during the 21st century under climate change conditions (RCP8.5), compared to historical 20th century and pre-industrial era simulations (Figure 2 and Figures S2-S10 in Supporting Information S1). We find fewer changepoints for smaller (phyto-and zoo-) plankton than larger plankton and fewer changepoints for phytoplankton than for zooplankton (Section 3.2), and that changepoints in plankton are more associated with iron and silicate than phosphate and nitrate, and more likely to occur where temporal shifts in latitudinal temperature gradients (a proxy for ocean fronts) are large (Section 3.3). However, the spatial extent of "hotspot" locations, which have more frequent changepoints, reduces from the pre-industrial to the historical simulation and again from the historical to the future simulation (Section 3.4). (As we use the kurtosis of the probability distribution of the number of changepoints to quantify this, it is not specific to e.g., ≥3 or ≥4 changepoints per century. Here by hotspots we mean locations with ≥3 or ≥4 changepoints per century; the number chosen as the minimum per century does not affect our conclusions, but locations with ≥5 changepoints per year are very rare for some variables in some simulations.)

Spatiotemporal Changepoint Patterns
We first quantify and describe the spatial extent and distribution of changepoints in temperature, nutrient concentrations, and plankton populations across PI, 20th, and 21st century simulations (RCP8.5). In terms of changepoint spatial extent, four main features are apparent. The fraction of the upper ocean with changepoints increases for all variables only slightly from the PI to 20th century, but increases dramatically from the 20th to 21st century (Figures 2, 3 and Figures S2-S10 in Supporting Information S1). Temperature has a smaller increase in changepoint extent between the 20th and 21st centuries than the other variables ( Figure 3; the same is true to a lesser extent for iron). Nutrients other than iron have a much higher fraction of ocean area with changepoints that become widespread in the 21st century (>60% of ocean area affected) compared to the other variables. Finally, within each simulation the plankton populations show similar fractions of ocean area with changepoints, and the ranks are relatively stable across centuries and consistent with the scaling of generation time with body size (e.g., Barton et al., 2020). For example, small phytoplankton have changepoints in the smallest fraction of ocean area while large zooplankton have the largest fraction of ocean area with changepoints in each simulation. A plot of changepoint area for the 21st century ( Figure 4) shows that the temporal distribution of these changepoints is somewhat uniform, with a peak in the 2060s ( §3.5) and secondary peaks in the 2030s and 2010s.
In terms of spatial distributions, temperature changepoints are confined mostly to the high latitudes and Northeast Atlantic in the PI century, expanding equatorward in the 20th and 21st centuries ( Figure S6 in Supporting Information S1). All the plankton exhibit variations of this pattern ( Figures  S2-S5 in Supporting Information S1). Temperature changepoints in the PI and 20th centuries are primarily in Southern Ocean locations strongly influenced by interannual variability (Auger et al., 2020;Behrens et al., 2021). This region of the Southern Ocean has strong, approximately zonal fronts in sea surface height that create strong latitudinal gradients in temperature, nutrients, and ecosystems (Chapman et al., 2020). Observations from recent decades suggest that these fronts are moving poleward in response to climate change (Kim & Orsi, 2014;Sallée et al., 2008;Sokolov & Rintoul, 2009). The prevalence of changepoints in the Northeast Atlantic may be due to expansion and contraction of subtropical gyre extent (Bograd et al., 2004;Irwin & Oliver, 2009;Polovina et al., 2008) or variability in the wind-driven gyre circulation in this region, which is particularly variable through time (Häkkinen & Rhines, 2004;Hátún et al., 2005). Nitrate changepoints are, by contrast, concentrated in the subtropics in the PI and 20th centuries before becoming nearly global in the 21st century (RCP8.5; Figure S7 in Supporting Information S1). Nitrate changes primarily occur along the tropical fronts and in the subtropical gyres (associated with very different nutrient regimes; Polovina et al., 2017), suggesting these fronts' interannual movements promote changepoints in upper ocean nitrate concentrations. Strong zonal fronts in environmental and ecological conditions that move meridionally on decadal to centennial timescales appear to be linked to changepoints in the model.

Trophic Changepoint Patterns
We also find an increase in the prevalence of changepoints from small to large phytoplankton, from small and medium to large zooplankton, and from phytoplankton to zooplankton, suggesting that changepoints are more widespread in larger organisms and higher trophic levels ( Figure 5). Environmental variations are filtered through marine food webs, such that longer-lived organisms tend to exhibit more pronounced low-frequency variability than do smaller, shorter-lived organisms (Barton et al., 2020;Di Lorenzo & Ohman, 2013). In other words, the greater area of changepoints in larger organisms and higher trophic levels may be caused, in part, by intrinsic trophic amplification within food webs (Chust et al., 2014;Stock et al., 2014a), which exacerbates extrinsic environmental and climate forcing.

Drivers of Plankton Changepoints
Plankton changepoints in the PI and 20th century can be understood as switches between alternate states associated with movements of frontal positions and co-located shifts with certain nutrients. These relationships break down in the 21st century as these sorts of changepoints and/or the relationships between them are likely eclipsed by the emergence of climate-changedriven trends that appear to emerge differently for different variables and suppress pre-existing patterns of variability (either entirely or to the point of statistical insignificance). In the PI and 20th century, we see that the correlations of different plankton changepoints' locations with each other are fairly strong ( Figure 6).
Spatial correlations in the frequency of changepoints here capture whether changepoints occur in each variable at the same frequency and in the same locations. In other words, the correlation ρ ij of variable i with variable j is calculated by defining a vector for both variables and of the number of changepoints that occur at each grid point k for that variable, and then computing the weighted Pearson correlation for those vectors, weighted by the area of each grid point k. This is preferred to a simple binary association such as a phi coefficient (Cramir, 1946) because variables can have multiple changepoints at a single location across a given time series.
In terms of nutrients, we find that plankton changepoints are most strongly correlated with changepoints in iron, then silicate, then phosphate, then least of all nitrate ( Figure 6). Correlations with temperature changepoints are moderate (Figure 6). These correlations between for example, iron and plankton changepoint locations in the PI and 20th centuries suggest that an appreciable fraction of these plankton changepoints reflect changepoints in these environmental variables, or alternatively similar phenomena specific to these locations are driving changepoints in correlated variables. These correlations also suggest that nitrate's changepoint dynamics are quite different than those for the other variables, because the locations of nitrate changepoints are only very weakly correlated with those of other variables across all simulations. This may in part be due to the greater complexity of the nitrogen cycle than that of other nutrients.
We also find that changepoints in the PI and 20th century tend to be associated with switches between alternate states; Figure 7 shows that in these centuries, most consecutive changepoints at a given location disagree in sign, that is, if one changepoint is associated with an increase in a baseline concentration or trend, the subsequent changepoint is usually associated with a decrease in that baseline concentration or trend. This may be because changepoints in the PI and 20th centuries are strongly associated with shifts in the position of ocean fronts,  Changes between PI and 20th century simulations (left) and the 21st century simulation for all correlations (except nitratephosphate, nitrate-small zooplankton, and small zooplankton-medium zooplankton) are significant (p < 0.01). which meander on interannual to multidecadal timescales. Figure 8 shows the frequency of changepoints, for a given century and variable, as a function of the temporal change in the latitudinal temperature gradient, ∂|∂T/∂y|/∂t. Large latitudinal temperature gradients |∂T/∂y| are commonly a proxy for the positions of ocean fronts; temporal changes in these are then indicative of changes in frontal positions over time. (We evaluated this for decadally averaged temperature values because the model's latitudinal temperature gradients are noisy on annual timescales.) For a given variable and century, a relative probability of changepoint occurrence of, for example, 2 for a given value of ∂|∂T/∂y|/∂t, for example, 0.05°C/°N/decade, means that changepoints are twice as likely to occur, relative to the global average, at places where and times when ∂|∂T/∂y|/∂t = 0.05°C/°N/decade. The strong increase in this relative probability with ∂|∂T/∂y|/∂t in the PI and 20th century for all plankton variables, temperature, silicate, and iron indicates that changepoints in these centuries and variables are associated with shifts in frontal positions.
In the 21st century, however, the correlations between plankton changepoints' locations and those of the environmental variables all decrease to near zero (Figure 6), the association between temporal changes in frontal positions all but disappear (Figure 8), and the sign-disagreement of consecutive changepoints is reversed such that consecutive changepoints tend to be reinforcing in most cases (Figure 7). Altogether this underscores that changepoints in the 21st century are of a qualitatively different nature than those of the PI and 20th century, neither due to changes in frontal positions nor switching between alternate states. This indicates that climate-change-driven changepoints disrupt pre-existing interannual-variability-driven changepoint patterns of plankton ecosystems, raising the question of how best to describe plankton changepoint patterns in the future if patterns of variability are distinct from the past. Furthermore, the lack of co-location in our model of plankton and nutrient changepoints draws into question the extent that nutrients may be indicators of multidecadal changes in plankton ecosystems.

Changepoint Hotspots
In addition to quantifying which locations experience and do not experience ecological changepoints in which variables, we also consider how frequent changepoints are in locations that have them. While the average numbers of plankton changepoints per century (in locations where there are changepoints) remains fairly stable ( Figure S11 in Supporting Information S1), more surprisingly we found that plankton (and temperature) changepoint "hotspots" (locations with ≥3 or ≥4 changepoints per century) tend to disappear over time. In other words, in places where changepoints are frequent, the rate at which changepoints occur slows down. This is most simply quantified by the excess kurtosis (κ = μ 4 /σ 4 − 5.99, where 4( ) = [ ( − ) 4 ] ) is the fourth central moment and σ is the standard deviation) of the changepoint distribution, where 5.99 is used instead of the usual value of 3 because the number of changepoints is non-negative; 5.99 is the kurtosis of an integer-rounded half-normal distribution ( Figure S12 in Supporting Information S1). The kurtosis is a standard measure of a distribution's heavy-tailedness; a positive κ indicates that the distribution has a heavier tail than an integer-rounded half-normal, and the larger the κ, the heavier the tail. In this context, a heavy tail means that changepoints tend to be concentrated in "hotspot" locations where changepoints occur frequently (Figures S12 and S13 in Supporting Information S1). Note that the kurtosis κ is only a measure of the probability distribution of the number of changepoints and is therefore not affected by the total ocean area experiencing changepoints or the total number of changepoints. Figure S13 in Supporting Information S1 also shows the changepoint probability distribution for small phytoplankton for each century to illustrate what a smaller excess kurtosis means in terms of the disappearance of changepoint hotspots. Figure 9 shows the excess kurtosis κ for each variable and simulation. While all variables besides nitrate are somewhat heavy-tailed in the PI century, the most heavy-tailed being small and large zooplankton, we find a decline in excess kurtoses across the plankton types (and for temperature) in the 20th and 21st centuries (RCP8.5). This demonstrates the disappearance of changepoint hotspots for all plankton (and for temperature), even after accounting for the increase in ocean area experiencing changepoints (Figure 3 and Figure S14 in Supporting Information S1). Changepoint hotspots for large and medium zooplankton, and large phytoplankton mostly disappear in the 20th century, whereas small phytoplankton and zooplankton and temperature show declines between both the PI and 20th centuries and between the 20th and 21st centuries. Note that this is not driven by an increase in the area with a low number of changepoints, that is, locations switching from having no changepoints to having low numbers of changepoints over time, because (a) the means of these distributions do not change substantially or consistently over time ( Figure S11 in Supporting Information S1), (b) the total area of locations with 4+ changepoints also decreases with time for all plankton (from between 0.88% and 2.4% in the PI century to between 0.08% and 0.87% in the 21st century; Figure 2 and Figures S2-S10 in Supporting Information S1), and (c) this decrease in excess kurtosis also holds when restricting only to locations that have changepoints in the PI or 20th century ( Figure S14 in Supporting Information S1).
For plankton and temperature, these hotspots all occur in the polar oceans and in the northeast Atlantic (though different variables' hotspots are not always perfectly co-located). Hotspots likely disappear with climate forcing because plankton communities that were previously switching frequently between alternative states no longer do so in the novel environmental conditions into which they are pushed. Temperature changepoint hotspots, concentrated in the high latitudes in the PI and 20th centuries, likely vanish with the poleward recession of seasonal sea ice in the 21st century out of latitudes where temperature has sufficient interannual variability. Nutrients are the exception in Figure 9, with κ values near zero or negative in the PI and 20th centuries, indicating a comparative absence of hotspots. Changes across centuries in nutrients' κ values are smaller, increasing slightly overall from PI to the 21st century. The lack of systematic or substantial change in nutrients' κ values likely reflects that these hotspots' shifts are not eclipsed by climate-change-driven trends. Plankton populations draw nutrient Figure 9. Excess kurtosis of the distribution of the number of changepoints per unit area per century for each variable. High values of excess kurtosis indicate spatial clustering of changepoints into "hotspots". Figure S12 in Supporting Information S1 shows how a reduction in kurtosis corresponds to a reduction in a distribution's tail; Figure S13 in Supporting Information S1 shows an example of a changepoint distribution's change across simulations; Figure S14 in Supporting Information S1 shows that the decrease in kurtosis seen here holds when only considering locations with changepoints in the PI and/or 20th century simulations for each variable. 9 of 11 concentrations down to subsistence concentrations, and therefore nutrient concentrations should have a threshold-like behavior (Tilman, 1982). Food web adjustments that stabilize changes in nutrients may be too small to be identified as changepoints, such that nutrients retain their changepoint hotspots despite climate-change-driven shifts but plankton do not.

Conclusion
We explored the spatial distribution and frequency of changepoints in a plankton community model within an ESM over three centennial simulations representing the pre-industrial period (the final 100 yr of a 1,500 yr control simulation), 20th century , and climate warming conditions in the 21st century (RCP8.5, 2001(RCP8.5, -2100. Anthropogenic forcing in the 21st century results in a substantial increase of ocean area where plankton abundances have changepoints, relative to the preindustrial era and 20th century. However, changepoint hotspots in the Southern Ocean and Northeast Atlantic Ocean, where plankton concentrations or temperature frequently undergo changepoints, largely disappear from the preindustrial period to the 21st century, as plankton communities that were previously switching frequently between alternative states no longer do so in the novel environmental conditions into which they are pushed. In agreement with the hypothesis of trophic amplification, larger plankton have more changepoints than smaller plankton, and zooplankton have more changepoints than phytoplankton; changepoints also often propagate through pelagic food webs rather than being restricted to individual plankton types. Plankton changepoints in the preindustrial and 20th centuries are associated with changes in frontal positions, certain nutrients (especially iron and least of all nitrate), and switches between alternate states. None of this is the case in the 21st century, when instead climate-change-driven changepoints disrupt pre-existing interannual-variability-driven changepoint patterns of plankton ecosystems. These results suggest that globally, plankton populations are susceptible to abrupt changes as a result of anthropogenic climate change, and that population dynamics are important for such changes, with larger organisms and higher trophic levels being more susceptible.
We note though that while latitudinal temperature gradients and nutrient concentrations are critical drivers of plankton communities, they are only some of many environmental factors that influence them; rapid plankton community responses associated with gradual warming may still reflect rapid changes in other environmental factors. We also emphasize that we have only considered one ESM here, and our results are specific to this model and likely differ in particular patterns if not mechanisms for other models. It would be instructive to test whether similar results and mechanisms underlying changepoints hold for other ecosystem models, and even more so whether complex ecosystem models that are able to predict properties such as phytoplankton diversity show similar results, and how these ecosystem properties are projected to change in the future. As the propensity for rapid ecological changes may either increase or decrease with system complexity (McCann, 2000), investigating the susceptibility of model populations to rapid changes across a range of ecosystem complexities will be essential to assessing the implications of model results for marine ecosystems (Cael et al., 2021). The approach we use here may also be suitable for terrestrial systems as well, particularly when considering differences in changepoint extent and frequency across organism size and trophic level at regional and larger spatial scales.

Data Availability Statement
The model outputs used here are available at DOI 10.5281/zenodo.6516343.