Widespread changes in Southern Ocean phytoplankton blooms linked to climate drivers

Climate change is expected to elicit widespread alterations to nutrient and light supply, which interact to influence phytoplankton growth and their seasonal cycles. Using 25 years of satellite chlorophyll a data, we show that large regions of the Southern Ocean express significant multi-decadal trends in phenological indices that are typically larger (<50 days decade–1) than previously reported in modelling studies (<10 days decade–1). Although regionally dependent, there is an overall tendency for phytoplankton blooms to increase in amplitude, decline in seasonality, initiate later, terminate earlier and have shorter durations, except in the ice, which initiate earlier and have longer durations. Investigating relationships with prominent climate drivers highlights regional sensitivities and complexities of multiple interacting aspects of a changing climate. Seasonal adjustments of this magnitude at the base of the food web can de-synchronize energy transfer to higher trophic levels, threatening ecosystem services and impacting global climate by altering natural CO2 uptake. Using 25 years of satellite chlorophyll a data, the authors demonstrate significant and widespread changes in the amplitude, timing, duration and seasonality of Southern Ocean phytoplankton blooms. Such changes threaten ecosystem services and can impact global climate by altering natural CO2 uptake.

Extended Data Fig. 1) but that may become so when more years of data are available. Contrary to what has been suggested in modelling studies 24,26 , significant trends in bloom magnitude (maximum, mean and integrated chl-a) were detectable from only 25 years of satellite data. Moreover, bloom magnitude metrics (maximum, mean and integrated chl-a) generally display a larger percentage significance than phenology metrics (initiation, termination, duration and peak timing). This disparity between model predictions and satellite observations may stem from the model's inability to accurately characterize the timing and amplitude of the Southern Ocean's seasonal cycle 30 . The predominance in positive trends in bloom magnitude (Figs. 2a-c and 3 and Extended Data Fig. 1) would arguably be associated with an increase in the amount of bulk carbon being produced during the seasonal bloom (for example, mean of positive significant trends in bloom integrated chl-a = +15 mg m -3 bloom -1 decade -1 , equating to an increase of 2.91% yr -1 ; Supplementary Tables 1 and 2) and thus beneficial to both carbon export and fuelling the food web. However, an increase in chl-a is not necessarily associated with an increase in biomass, particularly in the iron and light co-limited Southern Ocean, where phytoplankton can adjust their cellular chl-a to carbon ratios in response to environmental stresses 16 . This ambiguity in translating trends in chl-a to trends in biomass is highlighted by the negative trends observed in net PP for much of the Southern Ocean (for the majority of satellite models) 16,31 and the negative trends in backscatter (b bp ) (as a proxy for carbon) for the STSS and SPSS biomes (Extended Data Fig. 2a), which accounts for the observed general increasing trend in chl-a/b bp ratios (Extended Data Fig. 2b). Regardless, these trends in chl-a magnitude, whose positive distribution dominates in all biomes ( Fig. 3 and Extended Data Fig. 1), reflect strong adjustments in phytoplankton in response to environmental conditions. Interestingly, the ICE biome has an increasing trend in b bp (coincident with chl-a) in all but the Ross Sea, which is a highly productive and important region for the marine food web, with declining trends instead having the potential to negatively impact higher trophic levels.
Of note when interrogating the phenological trends is their magnitude (typically <50 days decade -1 ; Fig. 2d-f), which is substantially larger than previous estimates from modelling studies (≤10 days decade -1 ) (ref. 26). Indeed, only <6% of significant trends fell within a 10 day decade -1 range. Two example pixels (Extended Data Fig.  3a,b) provide a visual representation of a delayed bloom initiation and advances in maximum peak timing of 47 and 79 days decade -1 , respectively. The STSS and SPSS depict a predominance in the percentage significance of retreating bloom initiations (mean = -33 and -20 days decade -1 , respectively; Fig. 3 and Supplementary Table 1), which is contrary to what was predicted in two modelling studies that instead suggested a more regionally consistent trend (2006-2100) of advancing bloom initiations under IPCC representative concentration pathway 8.5 conditions 26,32 . These characteristics of retreating bloom initiations are coincident with a predominance in the percentage significant distribution of an advance in peak timing and termination, which drives a substantial contraction of bloom duration in the STSS and SPSS (mean = -53 and -20 days decade -1 , respectively; Figs. 2f and 3 and Supplementary Table 1). On the contrary, a predominance in the distribution of advancing bloom initiations in the ICE biome drives extended bloom durations (mean of significant positive trends in duration = +21 days decade -1 ; Figs. 2d,f Fig. 3 and Supplementary Table 1). The sensitivity of these trends to the choice of method for detecting phenological metrics was investigated by comparing these trends with those derived from the cumulative-sum and rate-of-change methods 33 . Although in a mean sense some differences are evident in the regional distribution of seasonal metrics (Extended Data Fig. 4), when it comes to the trends, their distribution (Extended Data Fig. 5) and percentage significance (Extended Data Fig. 6) are similar across all methods.
Research emanating from the past decade has emphasized the important role that subseasonal temporal scales and meso-to the effects of which may be synergistic or antagonistic, making the net impact on phytoplankton diversity, abundance, productivity and export complex and difficult to determine 19 .
The seasonal cycle sets much of the environmental variability in the factors that drive PP 20 . It is also the mode of variability that couples the physical mechanisms of climate forcing to ecosystem response in production, diversity and carbon export 21 . Accordingly, it is expected that any long-term trends in Southern Ocean productivity will be mediated through changes in the characteristics of the seasonal cycle and its interaction with the phenology of the ecosystem. Interrogating long-term alterations in the characteristics of the seasonal cycle is thus expected to provide a sensitive index of climate variability. The most widely used index of phenology is the start of the seasonal bloom, which can significantly impact the success of higher trophic levels that rely on plankton as their principal food source 22 . Other important metrics relate to the bloom amplitude and duration, which dictate the amount of biomass being generated within a season that can be exported to the ocean's interior or transferred to higher trophic levels via the marine food web 23 . Although modelling studies 24,25 have shown that 30-40 years of PP data are necessary to distinguish a climate change trend from background natural variability, others suggest that indicators such as phenology may detect trends more rapidly 24 . Indeed, a follow-up modelling study 26 found that from a global perspective, as little as 20 years of data were required to distinguish a trend in bloom initiation.
Remote sensing of ocean colour is currently the only observational capability that can provide synoptic views of upper-ocean phytoplankton at high spatial and temporal resolution and high temporal extent. In 2011, Thomalla et al. 27 used nine years of satellite chlorophyll a (chl-a) data (available at the time) to characterize the seasonal cycle of Southern Ocean phytoplankton in terms of bloom initiation, amplitude and variability. With an additional 16 years of data, we now have the means to not only characterize regional variability in seasonal indices, but also to investigate trends in those indices over the past 25 years. Understanding the spatial variability in multi-decadal trends in the characteristics of Southern Ocean blooms will provide a unique opportunity to assess regional differences in the response of the biogeochemical system to a changing climate.

Trends in phytoplankton seasonal metrics
Using a range of criteria, Fay and McKinley 28 divided the Southern Ocean into three biomes: the subtropical seasonally stratified (STSS), the subpolar seasonally stratified (SPSS) and the ice (ICE). A large degree of spatial variability is evident in the mean distribution of seasonal metrics (bloom timing, amplitude and seasonality) derived from 25 years of chl-a data, which oftentimes differ between ocean basins (for example, shorter blooms of lower amplitude in the Pacific), zonally (for example, associated with the different biomes) and regionally (for example, at continental margins and subantarctic islands) (Fig. 1). This heterogeneity is considered to result primarily from distinct underlying physics that alters the light environment and dominant supply mechanisms of key nutrients, driving subsequent regional variability in phytoplankton growth and variegated bloom characteristics 27 . Determining how large-scale physical processes associated with climate adjustments manifest as ecological responses at smaller scales is critical for interpreting ecosystem variability and deciphering trends and trajectories in the BCP and trophodynamics. However, spatial and temporal irregularity in the characteristics of marine ecosystems makes detection of systematic changes in response to climate pressures challenging 29 .
In this study, significant trends were detectable in all seasonal metrics, with spatial heterogeneity evident as regionally cohesive distributions of significantly positive or negative trends (Figs. 2 and 3 and Supplementary Table 1). This directional dominance in trend distribution is similarly evident in trends that are not significant (P > 0.05; Article https://doi.org/10.1038/s41558-023-01768-4 submeso-spatial scales play in characterizing the seasonal cycle of upper-ocean physics and biogeochemical response in the Southern Ocean [34][35][36][37][38][39][40][41][42][43][44] . The degree of seasonal cycle reproducibility (SCR) (calculated as the correlation between the observed annual time series and the mean climatological seasonal cycle) 27 is a useful metric that can capture these important scales of variability and reflect the characteristics of dominant seasonal versus subseasonal supply mechanisms 27 . An example pixel (Extended Data Fig. 3c) provides a visual representation of a decline in SCR of 29% decade -1 . Trends in SCR in the STSS and ICE biomes (Figs. 2h and 3, Extended Data Fig. 1 and Supplementary Table 1) suggest a predominance in declining SCR (mean of significant negative trends = -12% decade -1 and -20% decade -1 , respectively),  which implies a shift in the characteristics of the seasonal cycle to one that is more variable and less predictable. On the contrary, the SPSS biome suggests a predominance in trends of increasing SCR (mean of significant positive trends = +14% decade -1 ) and adjustments towards a less variable seasonal cycle. Two modified approaches to calculating SCR (as the correlation against a detrended and 5 yr rolling climatology) tested the robustness of the trends and susceptibility of the method to possible bias. Although some discrepancies exist in the climatological distribution of mean SCR (Extended Data Fig. 7b), there is little difference in the proportional distribution of observed trends, which retain Article https://doi.org/10.1038/s41558-023-01768-4 a relative predominance in declining SCR (Extended Data Figs. 7d and 8). These results suggest that the seasonal characteristics of vast areas of the Southern Ocean are evolving to reflect changes in the degree of environmental variability that may lead to alterations in biomass and phytoplankton diversity (for example, from competitive exclusion 45 ). Climate-induced shifts in phytoplankton species composition are key as cell size and elemental stoichiometry impose fundamental constraints on growth rates, food web structure and trophic-level interactions that determine the trajectory of carbon through the food web and the proportion of biomass being exported 46 . An example is the observed decline in mean size and a reduction of diatoms in the western Antarctic 47 associated with a thinning and retreat of the ice sheet over the past decade 48 . These changes are particularly important in the Southern Ocean, where many of the world's water masses are formed and the biogeochemical cycling of organic material dictates surface water nutrient supply 7 .

Linking climate drivers to trends in seasonal metrics
Of the most prominent climate-driven adjustments to Southern Ocean conditions are a general warming 49 , an intensification and southward shift of the westerly winds in association with a more positive phase of the Southern Annular Mode 50 and an increase in stratification and deepening of the summer MLD 15 . All phenological trends reflect the integrated impact of a suite of concurrent physical, chemical and biological processes that may be nonlinear, synergistic or antagonistic and may reflect correlation rather than causation. This makes it challenging to attribute observed trends in the characteristics of the seasonal cycle to a specific mechanistic driver. Nevertheless, with these dominant drivers in mind (Extended Data Fig. 9), we investigate their relationship with bloom duration (which represents both initiation and termination), mean bloom chl-a (the largest percentage significance of the chl-a metrics) and SCR (Fig. 4) (all other correlation maps are available in Extended Data Fig. 10).
Results highlight the regional heterogeneity of driver relationships that do not typically hold across biomes. Trends in mean bloom chl-a were correlated with sea surface temperature (SST) for large regions of the Southern Ocean (Fig. 4a) (similarly noted by ref. 31), which may reflect enhanced metabolic activity in response to warming 51 . The correlations of mean bloom chl-a with summer MLDs (Fig. 4b) and windiness (Fig. 4c) were patchy, with the variegated response in the sign of the correlation probably reflecting regional dependence in the degree   (Fig. 4d,e), which is anticipated to be linked to trends in sea-ice concentration and extent 53 , the onset of ice melt 54 and advancing bloom initiations (Fig. 3c). Negative correlations between SST and trends in SCR typically align with positive correlations in windiness (Fig. 4g,i).
There is some evidence of regional coherence in the trends of windiness and deeper summer MLDs, most notably upstream of the Drake Passage and the Weddell gyre (Extended Data Fig. 9), with a concomitant decline in SCR (Figs. 2h and 4h,i) that may reflect the reported increase in the number and intensity of low-level cyclones around the South Shetland Islands over the past 30 years 27,55 .

Decadal changes to seasonal regimes
More than a decade ago, Thomalla et al. 27 concluded their study with a synthesis schematic that divided the Southern Ocean into a montage of four regions that summarized the varying responses of phytoplankton biomass to different seasonal regimes. We are now able to reproduce that same schematic to explore decadal changes (Fig. 5), with the most notable differences being a decrease in regions of low chl-a and an increase in the spatial extent of regions characterized by high chl-a and low SCR (Fig. 5d,e). There are substantial differences in the spatial distribution of these adjustments, with a strong zonal coherence (Fig. 5c,e) that may reflect the poleward displacement of fronts within the Antarctic Circumpolar Current over the past 20 years 56 . These results suggest an increase in the role of subseasonal temporal scales and meso-to submeso-spatial scales, which is in agreement with a significant increase in Southern Ocean eddy activity attributed to a strengthening of the wind stress since the early 1990s 49 . Physical modifications such as these are likely to alter the intra-seasonal characteristics of light and nutrient supply, thereby impacting phytoplankton growth. This deduction is supported by growing evidence suggesting that in addition to deep winter mixing 57 , storm-driven entrainment

Ecological impacts of widespread trends in seasonal metrics
These trends in the characteristics of the Southern Ocean seasonal cycle reflect changes in environmental conditions that alter the drivers of phytoplankton production, abundance and community composition, with feedbacks that threaten the ecosystem services they provide: sustaining biodiversity, fuelling the food web and mediating global climate through an altered efficiency of the BCP 19 . Support of higher trophic levels is of particular relevance in this region as phytoplankton PP sustains a rich and diverse food web dominated by krill, marine birds, seals and migratory whales. As top predators, Weddell seals integrate information across the entire trophic web with the isotopic composition of their diet being preserved in the archive of their pelts 60 . A stable isotope analysis by Huckstadt et al. 61 found a significant decline in the isotopic characteristics of the diet of Weddell seals in the Ross Sea over the past 100 years, which was linked to widespread changes in upwelled nutrient supply, PP and phytoplankton community structure. These findings emphasize the cascading impact of climate adjustments in environmental conditions through the food web from the trophic base to top predators. The predominance of shorter blooms for the Southern Ocean as a whole (Fig. 3) is expected to have negative implications for carbon drawdown and trophic supply. Indeed, the duration of the seasonal bloom is proposed to be more important for carbon storage and sequestration than the magnitude 47 as longer mealtimes for consumers lead to larger proportions of the bloom being stored in animals as opposed to being reworked through the microbial loop 23 . These phenological trends also have the potential to negatively impact the survival rate of zooplankton and larval fish populations if prey availability decouples from critical life stages 62  species whose timing of spring migrations and breeding phenologies relies heavily on the quantity and quality of marine food sources. For example, seabirds in east Antarctica are delaying their arrival at their colonies by 1.6 days decade -1 and their first egg laying by 0.4 days decade -1 . These delays are linked to a regional decrease in sea-ice extent and an increase in sea-ice duration, which together reduce the quantity and accessibility of food supplies in early spring 65 . Similarly, humpback whales integrate a multitude of environmental signals into their migratory decisions 66,67 with evidence of a one-month advance (1.5 days earlier per year) in departure from their Antarctic peninsula feeding grounds being linked to a reduction in krill abundance in response to the decreasing mass of the autumn ice sheet 68 . The consequences of changes to the timing, quantity and quality of their food source feed back into these higher trophic levels impacting their nutritional stress, reproductive success and survival rates; particularly if they are unable to synchronize their feeding and breeding patterns with that of their food supplies 65,69,70 . In addition to gradual changes engendered by global warming, climate forcing may produce rapid shifts in oceanographic and food-chain dynamics, driving alternative trophic pathways 70 . The absence of any definitive answers from the driver analysis highlights the complexities of multiple interacting aspects of a changing climate. Nonetheless, the trends in seasonal metrics observed here are anticipated to continue or accelerate as climate change drives further adjustments in Southern Ocean physico-chemical conditions. These investigations provide a unique assessment of the response of seasonal bloom characteristics to climate drivers highlighting regional sensitivities that demand our attention in the context of the ecological impacts that such changes are expected to elicit on the BCP and trophic interactions.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-023-01768-4.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate Article https://doi.org/10.1038/s41558-023-01768-4 if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.

Methods
Satellite-derived chl-a concentrations were obtained from the European Space Agency Ocean Colour Climate Change Initiative (https:// esa-oceancolour-cci.org ref. 71) at 4 km and 8 day resolution (v.6.0). Analyses covered the period from 4 September 1997 to 27 December 2022 for the Southern Ocean biomes as defined in ref. 28, the STSS, the SPSS and the ICE. To reduce missing data, chl-a values were first regridded to a regular 25 km grid through bilinear interpolation using the xESMF Python package 72 . The remaining gaps were filled by applying a linear interpolation scheme in sequential steps of longitude, latitude and time 73 using a three-point window. If one of the points bordering the gap along the indicated axis was invalid, it was omitted from the calculation; if two surrounding points were invalid, then the gap was not filled. Finally, the data were smoothed by applying a moving average filter of the previous and next time steps. For more details on this method, see ref. 74.

Seasonal metrics
Phytoplankton blooms typically manifest as a seasonal cycle, with a bloom initiation that identifies the timing of the ramp-up in phytoplankton growth and biomass accumulation followed by bloom peaks within the growing season (which could be multiple) and finally the bloom termination, which defines the end of the growing season 27,[75][76][77][78][79] . The calculation of the phenological indices of bloom initiation, termination and duration followed closely the methods of ref. 26, detailed in the following, but with some additional steps (see Supplementary  Fig. 1 for examples of the phenological indices described).
(1) Bloom slice: each pixel and year's 'bloom slice' was found by locating the timing of the climatological mean bloom maximum then centring the bloom slice per year that starts and ends at the preceding and following 6 months, respectively. (2) Bloom maximum chl-a: the bloom maximum was identified as the local maximum in chl-a that occurs within this bloom slice. (3) Bloom initiation: the bloom initiation date for each bloom slice was calculated by first determining the minimum before the bloom maximum and the range, as the difference in chl-a concentration between the bloom maximum and the preceding minimum. The bloom initiation was then determined as the first date after the pre-peak minimum that the chl-a concentration was greater than the minimum chl-a concentration plus 5% of the chl-a range. (4) Bloom termination: the bloom termination date for each bloom slice was similarly calculated as the first date after the bloom maximum that the chl-a concentration was less than the chl-a concentration at the minimum plus 5% of the range. However, to ensure that seasonal blooms with more than one peak could be accounted for, multiple bloom peaks were defined as a second, third or nth local maxima where the chl-a concentration reached at least 75% of the amplitude of the bloom maximum defined in (2) and was a minimum of 24 days (3 ×8 day time intervals) from the maximum peak. The additional peaks were identified with the Python SciPy 80 function 'find_peaks'. In these instances, the six-month time span for either initiation or termination was calculated from the first or last peak, respectively. (5) Bloom duration: the bloom duration was calculated as the number of days between the bloom initiation and termination dates. (6) Bloom mean chl-a: the mean chl-a over bloom duration was calculated as the mean chl-a concentration between the bloom initiation and termination dates. (7) Bloom integrated chl-a: the seasonally integrated bloom chl-a was calculated using the NumPy 81 trapezoidal function as the chl-a concentration integrated between the bloom initiation and termination dates.
Additional methods for calculating bloom initiation, bloom termination, bloom duration, bloom mean chl-a and bloom integrated chl-a were used to test the sensitivity of the choice in phenological methods.
These included the cumulative-sum method and the rate-of-change method 33 , where the cumulative-sum method is defined as the first moment the cumulative chl-a biomass surpasses 15% of the total biomass, and the rate-of-change method is derived as the first moment the first derivative surpasses 15% of the median derivative. Results from the other methods are presented in Extended Data Figs. 4-6 for comparison.
(8) SCR: the reproducibility of the annual seasonal cycle was calculated by first generating a climatological mean seasonal cycle by resampling the time series to a monthly climatological mean that was then averaged across all years to generate a climatology and subsequently interpolated to 46 data points (the same number of 8-day timesteps per year as in the original chl-a dataset). The SCR is then calculated as the Pearson's correlation coefficient of the annual seasonal cycle against the climatological mean seasonal cycle. A value of 100% is indicative of an annual seasonal cycle that is a perfect replica of the climatological mean, while a value of 0% implies no relationship between the annual seasonal cycle and the climatological mean. Unlike for seasonal metrics (1)- (7), for SCR no 3× time step (3 × 8 day) rolling mean was applied to the regridded and interpolated data. All phenological indices data are publicly available 82 .
Two adjusted methods of calculating SCR were performed to test robustness and susceptibility to possible bias. (1) The time series was detrended before calculating the 25-year climatology and then correlated with each year to determine the percentage SCR. (2) Each year was correlated against a rolling five-year climatology (instead of a 25-year climatology) centred on the year in question to determine percentage SCR. Note that this approach requires the first and second years of the time series to be correlated against a climatology that is generated from only three and four years, respectively (as do the last and second-to-last years of the time series). Results from these approaches are presented in Extended Data Figs. 7 and 8 for comparison.
The cyclical nature of the calendar poses a considerable issue that needs to be addressed when calculating means of phenological indices. For example, we need to avoid a situation where the mean bloom initiation between a year with a bloom in December (for example, day of year = 340) and a year with a bloom in January (for example, day of year = 10) is incorrectly calculated as an average bloom initiation date in June (for example, day of year = 175). To account for this, we used the Python SciPy 80 function 'circmean', which calculates circular means for samples in a range (correct mean = day of year 357).

Trend calculations
The cyclical nature of the calendar similarly poses a problem when examining trends in dates for bloom initiation, termination and the timing of the bloom maximum. To account for this, trends in dates were calculated from the difference in days between each time step and the mean date, that is, the anomaly from the climatological mean phenological date.
Any pixel whose time series had less than 50% of the data available was excluded from any statistical trend analysis. Before linear regressions were performed on the trends, the data were first tested for a normal distribution. If the data were normally distributed, then linear regressions were performed using the Sci-Kit 83 Huber Regressor (ε = 1.35). If the data were not normally distributed, then linear regressions were performed using the non-parametric Mann-Kendall test 84 . All linear regression statistics are reported at the 95% probability, P < 0.05. A sensitivity test of varying ε between 1.2 and 1.5 altered the proportion of trends significant across the Southern Ocean but never by more than 12% (Supplementary Fig. 2).
To examine potential regime shifts across the Southern Ocean, mean chl-a concentrations were calculated using log transformation. The limits for different regimes were a mean chl-a of 0.25 mg m -3 and a mean SCR r 2 of 40%. The four regimes were defined as follows: low chl-a + low SCR, low chl-a + high SCR, high chl-a + low SCR and high chl-a + high SCR.

Driver analysis
To investigate the possible response of phytoplankton phenology to annual changes in the occurrence of high-speed wind events, we include a windiness metric that counts the occurrence of high wind speeds, defined as the number of days per year in which the daily averaged wind speeds are greater than the 25-year 75th percentile of all wind speeds (the percentile is calculated over 1998-2022). This was computed using reanalysis wind speeds between 1998 and 2022 from the Japanese 55-year Reanalysis 85 . For SST, we used the Group for High Resolution Sea Surface Temperature (https://www.ghrsst.org/). MLD data were computed from Hadley EN4.2.2 gridded temperature and salinity profiles 86 after conversion to potential density using a gradient of 0.03 kg m -3 (ref. 87). SST, MLD and wind data were spatially correlated with phenological indices using Pearson correlation coefficients.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
All data used in this study are available publicly 82 . Surface chlorophyll a is available from the Ocean Colour-CCI dataset (v.6.0) at https:// esa-oceancolour-cci.org. Wind data used in this study are from the Japanese 55-year Reanalysis ( JRA-55-do); data are available at https:// esgf-node.llnl.gov/search/input4mips/. To search, select "Target MIP" = "OMIP", "Institution ID" = "MRI" and "Source Version" = "1.4.0" among tabs on the left side. The sea surface temperature data used in this study are available at https://www.ghrsst.org. The gridded temperature and salinity profiles used to derive the mixed-layer depth are available at https://www.metoffice.gov.uk/hadobs/en4/ download-en4-2-2.html. All phenology output data are available at https://doi.org/10.5281/zenodo.8087125 ref. 82. Article https://doi.org/10.1038/s41558-023-01768-4 Extended Data Fig. 10 | Regional distribution of the correlation between physical drivers and the remaining phytoplankton seasonal metrics (selected subset in Fig. 4). Maps depict the spatial correlation of (a,b,c) maximum bloom chlorophyll-a (chl-a), (d,e,f) integrated bloom chl-a, (g,h,i) bloom initiation, ( j,k,l) bloom termination, and (m,n,o) maximum chl-a date against sea surface temperature (SST) on the left, summer mixed layer depth (MLD) in the middle and windiness (that is number of days per year that wind was greater than the 25 year 75th percentile) on the right. Data were excluded if less than 50% of either time series per year was available. Correlations are performed against the corresponding 25 years of SST, MLD and windiness data (1998 -2022), with data only plotted if the seasonal cycle metric trend is significant.