In the subarctic North Pacific Ocean, circulation in the upper part of the water column is dominated by two major gyres, the Alaskan Gyre (AG) in the east and the Western Subarctic Gyre (WSG) in the west, outlining the subarctic Pacific (Fukuwaka et al., 2004). A combination of the geostrophic transport, Ekman pumping and vertical mixing within the gyres provides regular winter nutrient enrichment to the upper layer (Capotondi et al., 2005; Nakanowatari et al., 2017). As a result, the subarctic Pacific is one of the three HNLC zones (Harrison et al., 1999; Liu et al., 2004), where phytoplankton dynamics is generally controlled by grazing (Strom and Welschmeyer, 1991) and Fe supply into the epipelagic waters (Sarmiento et al., 2004). Most of the external iron is transported with the subsurface waters from the shelf, shallow enough to be brought up to the surface by winter upwelling and vertical mixing (Lam and Bishop, 2008). Shelf areas also receive Fe with sediment resuspension and land runoff, fuelled mainly by snowmelt in spring and glacier meltwater in summer (Davis et al., 2014; Crusius et al., 2017), primarily where the drainage basins are composed of Fe-bearing volcanic rocks like in the Kuril/Kamchatka region (Lam and Bishop, 2008). The macronutrient and Fe-rich waters critical for primary production (PP) are transported offshore by the mesoscale eddies formed regularly in the Gulf of Alaska (GoA) (Whitney et al., 2005), around the Kamchatka Peninsula (Rogachev et al., 2007), and in the convergence zone of the Kuroshio and Oyashio currents on the Japanese coast (Ueno et al., 2023). Another major source of Fe in the subarctic Pacific is atmospheric deposition, consisting of the Fe-rich dust transported regularly from the Asian deserts (Duce and Tindale, 1991), occasionally emitted volcanic ashes (Hamme et al., 2010), and relatively recently acknowledged Alaskan glacial dust (Crusius, 2021) and excess dust of anthropogenic origin (Pinedo-González et al., 2020; Hunt et al., 2024).

The above processes, coupled with seasonal light availability, shape phytoplankton biomass and phenology in the subarctic Pacific, and remain in a delicate balance susceptible to climate oscillations (Mueter et al., 2004). A key driver of ocean conditions in the subarctic Pacific is the Aleutian Low, a vast low-pressure atmospheric system that dominates the region in winter, controlling the strength of the westerly winds and impacting the winter convective mixing and nutrient replenishment in the euphotic zone (Goes et al., 2004). The strength of subsequent ocean forcing of the Aleutian Low is strongly related to climate indices, including the El Niño–Southern Oscillation (ENSO) (Ortiz-Tánchez et al., 2002; Zhang et al., 2019), the Pacific Decadal Oscillation (PDO) (Bond and Harrison, 2000), and the North Pacific Gyre Oscillation (NPGO) (Di Lorenzo et al., 2008).

To describe phytoplankton variability, satellite-derived Chl-a concentrations were retrieved from the 25 km 8-day Globcolour product (ACRI-ST, 2017), based on level‐2 weighted averaging of surface chlorophyll-a (AVW) (O’Reilly et al., 2000; Pauthenet et al., 2024), composed of the SeaWiFS (1997 – 2010), MERIS (2002 – 2012), MODIS-Aqua (2002 – present), VIIRS-NPP (2012 – present), and VIIRS-JPSS-1 (2017 – 2022). The Globcolour-AVW dataset was used since it appeared to be the best-harmonized product out of the several existing merged products for climate studies (Pauthenet et al., 2024). Radiometric observations by Pauthenet et al. (2024) indicated that the Globcolour-GSM dataset had been significantly affected by the drift of the VIIRS sensor, while the OC-CCI dataset showed biases due to higher inputs from MODIS. Furthermore, the PFT validation performed against an independent flow cytometry dataset collected in 2022 confirmed good agreement between the PFT estimates and in situ data (Eisner and Lomas, 2022; Konik et al., 2024). The 25 km pixel size was sufficient for capturing spatial patterns in the subarctic Pacific, and the 8-day temporal resolution was considered optimum to increase cloud-free information required for the phytoplankton phenology analysis based on previous studies (Siswanto et al., 2022; Pramlall et al., 2023; Konik et al., 2024). Phytoplankton phenology metrics in the years 1998–2022 were determined following the methods described in detail by Konik et al. (2024). Generally, the pixels covered by clouds for more than 60% of the time were excluded, and at all the remaining locations, pixel-by-pixel, missing data was interpolated in time on a yearly basis to avoid extrapolation in January and December, the two most clouded months in the study area (Brody et al., 2013). The obtained time series were smoothed with a three-point median filter (Lavigne et al., 2018), and the phytoplankton blooms were identified using the median timelines extracted from a 3x3 pixel window (Konik et al., in prep). The Chl-a threshold method was used since it is preferred for the match/mismatch analysis between the trophic levels (Brody et al., 2013), where the baseline was established 5% above the climatology Chl-a median in the years 1998 – 2022, and bloom was detected only when Chl-a exceeded the baseline for at least two weeks (two consecutive maps) (Soppa et al., 2016; Suchy et al., 2022). Finally, the following phenology metrics were determined ( Figure 1 ): (1) starting date of a bloom (Bstart), the day (Day of Year, DOY) when Chl-a exceeded the baseline for at least two consecutive measurements, with a particular emphasis on the first bloom of the year (FBstart); (2) end date of a bloom (Bend, DOY), which was the following date when Chl-a dropped below the baseline; (3) length of a bloom (BLen, the number of days), derived as the difference between the end and the beginning of a bloom; (4) the maximum Chl-a concentration during a bloom (Cmax, mg m^–3), and (5) the time (Tmax, DOY) when maximum Chl-a concentration was observed (Soppa et al., 2016; Krug et al., 2018); (6) length of the bloom season (SLen, the number of days) was defined as the time elapsed between the start of the first bloom and the end of the last bloom; (7) the number of bloom days during each bloom was summed and given as the total bloom length (TotBLen, the number of days) (Krug et al., 2018); and (8) the Chl-a values during the bloom time were added together to estimate the bloom magnitude (Bmagn, mg m^–3 year^–1) (Soppa et al., 2016).

The phytoplankton composition, or phytoplankton functional types (PFT) (Bracher et al., 2017), from 2002 to 2022, was determined following the methods by Xi et al. (2020), as updated in Xi et al. (2021), but adapted for our regional application as in Konik et al. (2024). The retrieval of the satellite-based PFT time series involved three key steps. First, based on the pigment proportions and composition obtained using High-Performance Liquid Chromatography (HPLC), we determined the phytoplankton composition at the sampling stations with a chemotaxonomic model (CHEMTAX). Next, the identified groups were clustered into six main PFTs, which could be distinguished in the satellite images: diatoms, dinoflagellates, cryptophytes, cyanobacteria, green algae=chlorophytes and prasinophytes, and hapto+pelago=haptophytes and pelagophytes. Finally, an Empirical Orthogonal Function-based model was developed and validated with an independent dataset (Eisner and Lomas, 2022; Weitkamp et al., 2024), allowing for the computation of PFT maps for the entire study area. These satellite-derived PFTs were obtained using remote sensing reflectance (Rrs; sr^–1) products from the GlobColour database across nine spectral bands centered at 412, 443, 469, 490, 531, 555, 645, 670, and 678 nm, and the Sea Surface Temperature (SST) daily maps at 0.25°x0.25° NOAA OISST v2 high-resolution data (Reynolds et al., 2002, 2007) provided by the NOAA/OAR/ESRL PSL, Boulder, Colorado, USA (NOAA, 2007). Among the six PFT, the main focus of this study was the ratio between the diatoms and the small-cell phytoplankton, composed of haptophytes, pelagophytes, green algae, and cyanobacteria. Additionally, the diatom to dinoflagellate Chl-a ratio was included since they are often presented as examples of the opposite survival adaptations to physical and nutritional forcing, captured by Margalef’s mandala (Gilbert and Burford, 2017).

The analysis considered two groups of environmental variables: (i) 2D-time series of the environmental parameters obtained from hydrodynamic models, including Sea Surface Temperature (SST), salinity, Sea Surface Height (SSH), Mixed Layer Depth (MLD), surface current speed (CurrS), and satellite-derived Photosynthetically Available Radiation (PAR); and (ii) time series of climate indices.

The statistical significance of the long-term temporal trends in phytoplankton community structure in each bioregion was investigated using the non-parametric Mann-Kendall test (Hollander and Wolfe, 1973; Millard, 2013), and the abrupt changes in the phytoplankton phenology were determined using the changepoint method (Killick and Eckley, 2014; Killick et al., 2022), based on the pruned exact linear time (PELT) algorithm for the regime shift identification (Killick et al., 2012). The modified Bayesian information criterion (MBIC) was the penalty function of the changepoint analysis, guarding against small segments and allowing for any number of changes, including a ‘no change’ result.

To identify the key environmental factors driving phytoplankton phenology within bioregions, the multiple factor analysis (MFA) was performed using the factoextra R package (Kassambara and Mundt, 2020), which is based on the principal component analysis (PCA) (Le et al., 2008). Factors explaining 1% or less of the variability were excluded from the analysis and not presented in the results. The analyses were conducted with the numpy (Harris et al., 2020) and pandas (McKinney, 2010) Python libraries and R software (R Core Team, 2020) using the ggplot2 (Wickham, 2016), tidyverse (Wickham et al., 2019), gridExtra (Auguie and Antonov, 2017), RColorBrewer (Neuwirth and Maindonald, 2022), Kendall (McLeod, 2011), EnvStats (Millard, 2013), and Performance Analytics (Peterson et al., 2020) RCRAN packages.

All the analyses were performed at the scale of the bioregions, which were adapted from the original bioregions defined in Konik et al. (2024). For this, the original bioregions were subdivided into the eastern and western parts due to the impact of large-scale circulation patterns, namely the AG and the WSG, related to the spatial variability in the physical and chemical water properties in the subarctic Pacific (Harrison et al., 1999; Nishioka et al., 2021). Marginal and Connecting bioregions were split into the eastern and western parts as the previous study showed significant contrasts between the Eastern and Western Pacific with respect to hydrology, including the SST, and the zooplankton community structure (Chiba et al., 2015). For example, significant differences in the vertical stratification of the water column between the eastern and western North Pacific have been described in the literature (Old et al., 2019). Vertical stratification controls the mixed layer depth and, therefore, the nutrient concentrations in the upper ocean layer (Goes et al., 2004; Yasunaka et al., 2021); most importantly iron, which is the main nutrient limiting primary productivity within the HNLC zone in the subarctic Pacific (Zhang et al., 2021). The iron loads entering the WSG are significantly higher than those that supply the AG, mainly due to the water exchange with the adjacent seas, the Bering Sea and the Okhotsk Sea (Goes et al., 2004; Nishioka et al., 2020). Furthermore, the western subarctic Pacific is a highly dynamic region characterized by intense eddy activity and strong vertical mixing (Gaube et al., 2015; Siswanto et al., 2022). This is particularly evident in the Kuroshio-Oyashio Extension Region (KOER) near Japan (Itoh and Yasuda, 2010) and along the Kuril-Kamchatka Trench (Itoh and Yasuda, 2010), whereas the central part experiences relatively weaker vertical mixing (Gargett, 1991).

The resulting division consisted of the eight bioregions: Central South (CS), Central Oregon (CO), Central North (CN), Connecting East (ConE), Connecting West (ConW), Marginal East (MarE), Marginal West (MarW), and Archipelago (Arch) ( Figure 2 ).

The interannual changes in the main phenology metrics, Bstart and Bend, from 1998 to 2022, generally did not show statistically significant long-term trends for the bioregions. However, a substantial change in the average total bloom length (TotBLen) was observed in recent years. The changepoint analysis revealed that the TotBLen decreased by between six weeks (bioregion CO) and two weeks (Arch) ( Figures 3b, f ). The largest decrease in the TotBLen, of 45 days, was observed in CO, where the Chl-a is generally low (<0.7 mg m^–3). However, an evident drop in the TotBLen was also observed in the other bioregions where the Chl-a and related phytoplankton biomass is much higher (>1 mg m^–3). For example, in the border of the Gulf of Alaska (ConE and MarE bioregions), the TotBLen was shorter by a month after 2014 ( Figures 3d, h ). The only bioregion where the TotBLen increased by almost two weeks was the MarW bioregion, bordering the Kamchatka Peninsula ( Figure 3g ), where changes were observed as early as 2001. This was followed by a decrease in TotBLen in the Arch bioregion in 2004 and the central bioregions between 2010 and 2013, with the latest change in the ConW in 2018 ( Figure 3e ).

We found a relationship between the decrease in TotBLen and Chl-a. Using the Mann-Kendall test, we confirmed the statistical significance of the relationships between TotBLen and Cmax in bioregions CS, CO, ConE, MarE, and Arch ( Figure 4 ). Thus, we may expect further decreases in TotBLen in the future should there be a change in total Chl-a levels, for example, resulting from a shift in phytoplankton composition.

A general decrease in diatom and cyanobacteria abundance in all bioregions and an increase in the hapto-pelago group, except for the CS bioregion, were statistically confirmed by the Mann-Kendall test with the seasonal correction ( Table 1 ). No change in green algae was observed in the connecting bioregions ConE and ConW and marginal bioregions MarE, MarW, and Arch. Only a slight negative trend in dinoflagellates and cryptophytes was found in the ConE and MarE, bordering the Gulf of Alaska. As a result, negative trends in diatom-to-dinoflagellate and diatom-to-small algae Chl-a (haptophytes, pelagophytes, green algae, and cyanobacteria) were observed across the entire study area.

In the central bioregions, such as the CS, CO, and CN, the diatom-to-small-algae ratio anomaly was less pronounced than in the other bioregions due to relatively smaller input from diatoms ( Figures 5a–c ). However, between 2010 and 2013, negative anomalies were evident even in these bioregions ( Figures 5a–c ), which matches the timing of the regime shift in the TotBLen in the respective bioregions ( Figure 3a–c ). This match in time was also observed in the ConE and MarE bioregions ( Figures 5d, h ), where the negative anomalies co-occurred with the decrease in the TotBLen, around 2014 ( Figures 3d, h ), and in 2018 in the ConW bioregion ( Figures 3e , 5e ). The change in bioregions Arch and MarW was not as evident as in the others ( Figures 5f, g ). Still, the positive diatom-to-small-algae ratio anomaly peaks associated with diatom blooms were lower in the last decade, and periods with negative anomalies lasted longer, particularly after 2018.

Comparable results were obtained from the diatom-to-dinoflagellate ratio anomalies ( Figure 6 ), where the negative anomalies in the last five years were evident across all bioregions. Moreover, a similar match in time between the shift towards the predominance of the negative diatom-to-dinoflagellate ratio anomalies ( Figure 6 ) and the TotBLen decrease ( Figure 3 ) was observed, such as a change in 2014 in ConE and MarE ( Figures 3d, h , 6d, h ), and in 2018 in ConW ( Figures 3e , 6e ).

Regarding the environmental drivers, there were clear differences between the bioregions ( Supplementary Figures S2–S7 ). In the central bioregions, CS, CO, CN, and ConW, the spring and summer SSH variability was the primary descriptor of the Chl-a seasonal change. Of secondary importance was the spring salinity change in bioregions CO and ConW, suggesting two distinct water mass mixing ( Figures 7b, d ), and in CS and CN, summer and fall SST and PAR influence was more pronounced ( Figures 7a, c ). The MarE and ConE bioregions, bordering the Gulf of Alaska, were influenced mainly by the winter and spring SST changes, followed by the winter MLD, spring and summer SST, and SSH variability ( Figures 7e, h ). The MarW bioregion was most impacted by the salinity fluctuations and SSH in spring and summer ( Figure 7g ). In contrast, Arch was the only bioregion where the summer CurrS was an important factor ( Figure 7f ), supported by the summer SST and SSH variability.

Considering the climate indices ( Figure 8 ), the global ENSO pattern represented here by the MEI v2 indicated El Niño (MEI>1.5) in 1998, 2009/2010, and 2015/2016, and La Niña (MEI<1.5) in the years 1999 – 2000, 2007 – 2009, 2010 – 2012, and 2020 – 2022. There was also a shift to the negative NPGO at the beginning of 2014, which remained negative until 2022, except for a few months in 2014 and 2016. The NPGO and the PDO did not correlate during the analyzed period ( Supplementary Figure S1 ), where the negative PDO was observed in the years 2008 – 2009, 2010 – 2014, and 2020 – 2022.

The statistical results based on the Mann-Kendall test generally revealed three main patterns: (1) the large impact of the NPGO on the bloom phenology in the central CO and CN and the eastern ConE and MarE bioregions ( Table 2 ), (2) the strong influence of the El Niño-Southern Oscillation (ENSO—MEI v2) on the south-central CS and CO, and (3) the primary role of the Aleutian Low on the western part of the Pacific, including the ConW, MarW, and Arch, which was reflected in significant relationships between the phenology metrics, especially SLen and BLen, and the NPI, the ALBSA, and the PNA.

Our study revealed a general decrease between 2002 and 2022 in the diatom to small-cell phytoplankton Chl-a ratio anomaly ( Figure 5 ) and the diatom to dinoflagellates Chl-a ratio anomaly ( Figure 6 ). The decreasing trend was statistically significant in all bioregions, but the most pronounced anomalies were observed in the connecting and marginal bioregions, ConE and MarE ( Figures 5d, h ), where the biggest changes in the TotBLen were also observed ( Figures 3d, h ). A decrease in bloom-forming diatom contribution and a shift towards smaller phytoplankton could dampen the annual Chl-a variability and likely cause a change in the observed TotBLen. Previous phytoplankton composition analyses in the years 2012–2015 in the GoA and its margin suggested a drop in diatom production and a shift towards smaller phytoplankton groups, such as haptophytes, chlorophytes, and cyanobacteria, mainly during the marine heat wave (MHW) periods (Peña et al., 2019; Taves et al., 2022). The major MHW in the Northeast Pacific, where we observed the most significant change in diatom input, were the ones during the winters of 2013/14 (Bond et al., 2015) and 2014/15 (Di Lorenzo and Mantua, 2016), and summer 2014 and 2019 (Amaya et al., 2020; Song et al., 2023). Lower diatom abundance during the MHW, between 2014 and 2016, was noticed in the Continuous Plankton Recorder (CPRs) measurements in the GoA, among all measurements collected in the years 2000 – 2018 (Batten et al., 2022). Cael et al. (2021) also found the large-cell diatoms and dinoflagellates to be one of the most likely to experience an abrupt shift in biomass in the 2020s, and the more recent models confirmed a significant switch in community composition from diatoms to dinoflagellates in the Pacific (Arteaga and Rousseaux, 2023).

Our findings suggest that the reduction in diatom abundance in the northeast subarctic Pacific is not only a temporary effect of MHW, but may be part of a gradual change in the phytoplankton communities, confirming previous predictions of a decline in the diatom biomass (Boyd et al., 2004; Rousseaux and Gregg, 2015). Consequently, in the eastern bioregions with the highest contribution of diatoms, ConE and MarE, we observed the biggest change towards shorter TotBLen, which may impact the total primary production in the high latitudes of the subarctic Pacific Ocean (Tréguer et al., 2018). It is also important to stress that even a relatively small change in the diatom biomass can have a significant impact on marine export production (Waite et al., 1992) since they are responsible for roughly 40% of the particulate organic carbon export globally due to their fast-sinking silicate frustules and fecal pellets, resulting from the intense diatom grazing in the higher latitudes (> 50°N) (Jin et al., 2006).

The NPGO index has been widely recognized as strongly linked to the oceanographic conditions in the eastern subarctic Pacific since it describes the water flow intensity in the eastern and central branches of the North Pacific Gyres (Stramma et al., 2020). Here, we observed a positive relationship between the NPGO and the BLen in the eastern bioregions ConE and MarE, and the TotBLen in the central bioregions CO and CN ( Table 2 ). These positive relationships are likely associated with increased nutrient availability in the GoA region and the California Current System during positive NPGO regimes (Di Lorenzo et al., 2008). For these bioregions, the influence of the PDO was of secondary importance, which contrasts with the findings of Bond and Harrison (2000), who identified the PDO as a primary factor impacting the North Pacific. Recent studies have also reported a weaker impact of the PDO on the eastern North Pacific ecosystems (Kilduff et al., 2015; Litzow et al., 2018). The positive PDO phase is typically associated with above-average water temperatures in the eastern subarctic Pacific and cooler-than-average temperatures in the western part (Graham et al., 2021). However, in our study, the PDO was primarily correlated with the TotBLen in the ConW bioregion, located in the Kuroshio bifurcation zone ( Table 2 ). One potential reason for the relationship between BLen and the NPGO, instead of the PDO, is that both PDO and NPGO are influenced by El Niño events, albeit in different ways. Kilduff et al. (2015) noted changes in contemporary El Niño events, suggesting that recent El Niño occurrences are more frequently associated with central Pacific warming, which indirectly modulates the NPGO rather than causing the eastern Pacific warming that regulates the PDO. During the period covered in this paper (2002 – 2022), the relationship between the PDO and the NPGO was not statistically significant ( Supplementary Figure S1 ), confirming the difference in how the two are modulated.

Our findings also showed a direct impact of the ENSO pattern on BLen in the CO bioregion and Cmax in CS, whereas in the western bioregions, ConW, MarW, and Arch, SLen and BLen were more related to the NPI, the PNA, and the ALBSA indices ( Table 2 ). Some degree of intercorrelation between these three indices was observed: proportional between the NPI and the ALBSA and inversely proportional between the PNA and the ALBSA ( Supplementary Figure S1 ). Nonetheless, all three were included in the analysis because they have been widely discussed in the literature as part of a global teleconnection network (Schwing et al., 2010) and were shown to better describe various aspects of the natural variability related to the Aleutian Low together, especially when combined with the AO, than when used separately (Overland et al., 1999). The Aleutian Low — represented by the NPI — affected all the western bioregions, including the Arch, MarW, ConW, and CN. However, the SLen in the MarW bioregion, bordering the Kamchatka Peninsula, was influenced primarily by the ALBSA index, characterizing the ice melt in the Bering Sea region (Cox et al., 2019). Additionally, the relationship between the ALBSA and the AO, which is known to control the sea ice melt in the Okhotsk Sea, the snow melts, and the Amur River discharge around the Kamchatka Peninsula (Tachibana et al., 2008), indicates the critical role of the land runoff in the western Subarctic Pacific (Ogi and Tachibana, 2006), which was reflected in statistically significant relationships with SLen and BLen in the ConW, MarW, and Arch bioregions. These relationships confirm the primary role of salinity among the environmental factors ( Figure 7g ), which is unique to the MarW bioregion. Moreover, the Amur River carries exceptionally high loads of iron (Shamov et al., 2014), likely reflected in the highest observed Chl-a in MarW among all the bioregions ( Figure 4b ).

The Aleutian Low, captured by the NPI and PNA indices, also affected ConW and CN bioregions ( Table 2 ). Particularly in the ConW, located in the Kuroshio Bifurcation zone, the PNA influence on SLen and BLen was marked, which may be related to the higher frequency of tropical cyclone occurrence in the Kuroshio region (Song and Klotzbach, 2019), reflected in the SSH and SST variability ( Figure 7d ). In the neighboring Arch bioregion, SLen was also impacted by the NPI, in addition to the PNA and the ALBSA indices. Still, the main factors in this bioregion were the CurrS and SST in summer ( Figure 7f ), suggesting the primary role of the ocean currents (Abe and Nakamura, 2013; Nishioka et al., 2020). Positive SLen correlation with the NPGO indicates that the stronger mixing within the gyre likely introduces more nutrients to sustain longer blooms. The negative correlation with the NPI, on the other hand, suggests that stronger storms and diapycnal mixing may weaken water column stratification and lead to deepening of the mixed layer (Suga et al., 2004) below the critical depth described by Sverdrup’s theory (Smetacek and Passow, 1990), and resulting in shorter SLen. Lastly, the NPI and the AO influence on FBstart, SLen, and TotBLen was also noticed in the two largest bioregions, CN and CS ( Table 2 ). These bioregions are governed by the Aleutian Low, controlling the strength of the westerly winds and the wind-related stress over the North Pacific (Goes et al., 2004), which affects the mixed layer depth and the transport of Fe-rich dust from Asian deserts (Duce and Tindale, 1991), especially to the mid-latitudes (between 35° and 45°N). The key factors in the central subarctic Pacific were the SSH and PAR ( Figure 7c ), providing more nutrients and light to sustain primary productivity, which also agreed with the general patterns described by Goes et al. (2004).

Summary characteristics of each bioregion’s dominant environmental factors and climate indices were compiled in Table 3 .

Changes in the length of the TotBLen and diatom contribution will likely affect the entire ecosystem through impacts on zooplankton production, community composition and nutritional quality (Suchy et al., 2022; McLaskey et al., 2024). At a higher trophic level, Piatt et al. (2020) reported high mortality and reproductive failure of common murres (Uria aalge) in the years 2014–2016 linked to the abnormally low phytoplankton biomass and shift in the zooplankton community structure toward less-nutritious taxa. Nielsen et al. (2021) also revealed a significant change in the composition of the ichthyoplankton assemblages during the MHW along the GoA shelf. Nonetheless, the full picture may be more complicated, with Batten et al. (2022) not finding lower zooplankton abundance overall in the GoA during the 2014–2016 MHW, even though a decrease in the diatom-to-small algae ratio was reported. Although before the MHW in 2014, SST had been positively correlated with the abundance of diatoms and zooplankton biomass, during the MHW, zooplankton biomass did not follow the drop in diatom abundance (Batten et al., 2022). On the contrary, the high numbers of zooplankton, with the higher metabolic rates boosted by the elevated SST, were hypothesized to increase the grazing pressure on diatoms. The declines in cold-water zooplankton taxa and general species richness were also noted, indicating that the phytoplankton and zooplankton community structure shift is more probable than a drop in their biomass.

The positive correlation between NPGO and the phytoplankton phenology metrics, particularly BLen in the eastern bioregions, MarE and ConE ( Table 2 ), supports recent research showing a change toward a more pronounced impact of the NPGO on the west coast of North America (Kilduff et al., 2015; Puerta et al., 2019; Suchy et al., 2022) compared to the PDO which has been historically associated with variability of Pacific salmon (Mantua et al., 1997). Furthermore, recently the NPGO combined with the SOI has been shown to be strongly correlated with the recruitment abundance of Japanese sardine (Sardinops melanostictus) in the Kurosio — Oyashio Extension (Yatsu et al., 2021) of the western North Pacific.

All the above findings show certain changes in the food web structure, which appear to be associated mainly with the MHW period. However, the expected increasing frequency of MHW (Song et al., 2023) may change biological landscape permanently in the near future, and the NPGO might be a better predictor of the environmental variability directly affecting the food web in the subarctic Pacific.