Large diatom bloom off the Antarctic Peninsula during cool conditions associated with the 2015/2016 El Niño

Diatoms play crucial functions in trophic structure and biogeochemical cycles. Due to poleward warming, there has been a substantial decrease in diatom biomass, especially in Antarctic regions that experience strong physical changes. Here we analyze the phytoplankton contents of water samples collected in the spring/summer of 2015/2016 off the North Antarctic Peninsula during the extreme El Niño event and compare them with corresponding satellite chlorophyll-a data. The results suggest a close link between large diatom blooms, upper ocean physical structures and sea ice cover, as a consequence of the El Niño effects. We observed massive concentrations (up to 40 mg m–3 of in situ chlorophyll-a) of diatoms coupled with substantially colder atmospheric and oceanic temperatures and high mean salinity values associated with a lower input of meltwater. We hypothesize that increased meltwater concentration due to continued atmospheric and oceanic warming trends will lead to diatom blooms becoming more episodic and spatially/temporally restricted. Low meltwater input, cold temperatures, and high salinity accompanied a large diatom bloom off the Northern Antarctic Peninsula in 2015/2016 and suggests such blooms may be restricted under warmer conditions, according to direct measurements and remote sensing data.

P oleward shifts in distribution and composition of species are a major evidence of climate warming consequences [1][2][3] . Along the most warming-impacted regions of the Southern Ocean, with a myriad of cascading effects across all higher trophic levels, shifts from large to smaller phytoplankton cells have already been reported in marine ecosystems 1,4,5 . With long-term global climate warming, the West Antarctic Peninsula (WAP) has been experiencing one of the largest atmospheric warming rates on Earth since the 1950s 6 . Together with atmospheric warming, there have been long-term increases in upper ocean temperatures 7 , marine glaciers retreating 8,9 , and significant sea ice season reduction 10 . These long-term physical changes have led to substantial ecosystem alterations, namely in abundance and composition of phytoplankton groups across the WAP, although there have been marked meridional differences 1,4,5 . Whereas the polar conditions are persisting along the mid/southern WAP, with a southward increase in phytoplankton biomass and diatom abundance over time 1,3,11 , in the northern WAP region (hereafter Northern Antarctic Peninsula; NAP) these biological changes are particularly dramatic 1,12 . With its key role at the base of the food web, a downward trend in phytoplankton biomass due to a decrease (increase) in abundance of large diatoms (small-sized cryptophytes) has been observed in situ along the NAP 13,14 .
The NAP marine system, which encompasses the Gerlache and Bransfield Straits, the southernmost sector of the Drake Passage, and the northwestern Weddell Sea, is a diverse region featuring a range of oceanic and coastal ecosystems with complex hydrographic and biogeochemistry dynamics 12 . The surface NAP waters can be split into cold and warm variety. The cold regime, located east of the Peninsula Front and southeast of the Elephant and Clarence islands, is sourced by transition waters from the Weddell Sea, which are generally characterized by low macronutrient concentrations 15 . In turn, the warm variety is mainly derived from the Antarctic Circumpolar Current surface waters, being distributed northwards along the NAP western continental shelves and the west of the Bransfield Front. The warm conditions west of the Peninsula Front are enriched of macronutrients in subsurface waters as warm, nutrient-and carbon-rich Circumpolar Deep Water crosses the western Antarctic Peninsula continental shelves 16 . In addition, there is a northward surface water flow, from the Gerlache Strait to the warmer side of the Bransfield Strait, characterizing transition waters from the Bellingshausen Sea, which is responsible for advecting waters of a highly productive coastal zone to a more open region north of the NAP 12,17,18 .
The trends of atmospheric warming and sea ice declining along the NAP, which have intensified in the mid-twentieth century, have slowed down since the late 1990s 19,20 . This apparent recent short-term plateau of longer-term trends has been associated with regional wind pattern changes coupled with large-scale climate variability, such as Southern Annular Mode (SAM) and El Niño-Southern Oscillation 16 . During negative SAM/El Niño conditions, colder and weaker southerly winds blow across the region, resulting in high sea ice concentrations that reduce wind-mixing in winter and supply meltwater to upper ocean stabilization in spring/summer, enhancing phytoplankton-growth 3,10,21 . The opposite occurs during positive SAM/La Niña. Additionally, the effects of positive SAM may be offset by El Niño conditions, leading to colder winds and increased sea ice extent along the NAP 16 .
The 2015/2016 El Niño was one of the strongest on record, with sea surface temperature anomalies exceeding 2.5°C in the tropical Pacific Ocean 22 . Its effects were associated with high variability in sea ice extent along the Southern Ocean 23 , and an unprecedented increase in the ice sheet accumulation over the Antarctic Peninsula 22 . Moreover, in the same period, an intense bloom of large centric diatoms was registered during the NAP summer 24 , contrasting its long-term trend towards smaller phytoplankton cells 1,12 . However, a clear relationship between the results of the intense diatom bloom and the effects of the extreme 2015/2016 El Niño event was not established. Here, we investigate the link between this first extreme El Niño of the twenty-first century 25 and the atypical seasonal phytoplankton succession dominated by diatoms, with massive concentrations, during the spring/summer of 2015/2016 along the NAP, using collected in situ data and remote sensing measurements. This provides a rare opportunity to understand how diatom succession dynamics can positively respond to the effects of an extreme El Niño event, a counterpoint to drivers of long-term warming trend in a changing NAP.

Results and discussion
Due to contrasts in oceanographic properties along the NAP 24 , the sampling grid was split in two subregions: north and south ( Fig. 1; see "Methods"). The north and south subregions showed from the satellite data a much higher spring/summer (November-February) mean chlorophyll-a (Chl-a) in 2015/2016 than the decadal average time series (2010-2019; Table 1). In agreement with the El Niño effects 10,16 , the sea surface temperature (SST) and the air temperature showed substantially lower mean values during the spring/summer of 2015/2016 along the subregions (Table 1). However, there was an evident spatial/ temporal variability in sea ice concentration/duration between the subregions, with a northward (southward) lower (higher) mean value during 2015/2016 in relation to the decadal average (Table 1). Along the south subregion during the spring/summer of 2015/2016, the increased Chl-a during January followed the decline in the sea ice concentration over the spring and early summer, concurrent with increased SST, which was markedly colder throughout the seasonal phytoplankton succession (Fig. 2a). These results to the south subregion are consistent with previous studies along the WAP, in which years characterized by longer sea ice cover in winter have led to higher phytoplankton biomass in the following summer associated with a more stable water column 11,16,26 . To the north subregion, however, although there was a similar pattern between Chl-a and SST, the increased Chl-a during January was not related with the sea ice retreat (Fig. 2b). Moreover, there was a clear difference between the Chla peaks (the highest Chl-a value reached) along the subregions from the satellite data. The Chl-a peak in the south subregion occurred in early January (10 January 2016, reaching 1.73 mg m -3 ), whereas in the north subregion the Chl-a peak was observed in late January (29 January 2016, reaching 2.23 mg m -3 ).
It has been estimated that drifters entrained in the Gerlache Strait Current and the Bransfield Strait Current exit the Bransfield Strait in 10-20 days 17 , which is consistent with the interval of 19 days between both Chl-a peaks when considering the extreme distance between the subregions (see Fig. 1). These authors also estimated that drifters deployed in the Gerlache Strait Current were quickly advected out of the Gerlache Strait in less than 1 week (i.e., low residence time) 17 , which supports the similar diatom species assemblages identified in our microscopic analysis between stations of the Gerlache Strait and southwestern Bransfield Strait 24 . Therefore, it is plausible that phytoplankton growth in the north of the Gerlache Strait may be laterally advected northward into the Bransfield Strait, explaining the observed concomitant increase of satellite Chl-a data in both subregions from spring, associated with sea ice retreat southward (Fig. 2). In addition, as phytoplankton biomass tends to accumulate northward 17,27,28 , the advection processes could also explain the temporal and intensity differences of the Chl-a peaks along the subregions (see Fig. 2). This suggests that there was a link between the sea ice dynamics, phytoplankton biomass (Chl-a) and advection processes along the NAP during the spring/summer of 2015/2016, in which the sea ice melting first triggered an increase in phytoplankton biomass through water column stratification along the south subregion, and the advection processes led to a subsequent increase northward. Average of chlorophyll-a (Chl-a), sea surface temperature (SST), air temperature, and sea ice concentration during spring/summer (November−February) on a decadal average time series (2010−2019) and on 2015/2016 along the NAP subregions. The sea ice duration was also estimated to the decadal-mean and 2015/2016 (see "Methods"). Diff is the spring/summer difference value between the 2015/2016 and the decadal-mean. The positive (negative) values indicate that parameters were higher (lower) in the 2015/2016 spring/summer compared with the decadal average time series. Significant increase or decrease (-) of the parameters is indicated by asterisks, with a confidence limit of 95% (*) and 99% (**), using Student's t test comparisons. The sea ice duration was not tested (see "Methods"). The satellite Chl-a data require extensive validation with in situ data, especially in polar regions, where cloud cover is ubiquitous and performance is typically poor, due to not properly accurate Chl-a algorithms 12,29 . For that, although the mean Chl-a in 2015/ 2016 from the satellite data was approximately twice as large as the decadal average, there was a severe discrepancy in the mean Chl-a values observed between the in situ and remote sensing data (see Table 1 and Supplementary Table 1). This highlights the importance of the in situ dataset reported here, especially evident during February 2016, when the signal of an intense diatom bloom (> 40 mg m -3 Chl-a) 24 was not captured in the satellite data ( Supplementary Fig. 1), supporting that phytoplankton biomass accumulation during this summer was much higher than recorded by remote sensing observations (see Table 1). In general, the in situ Chl-a achieved its maximum (40 mg m -3 ) and higher mean value (17.4 mg m -3 ) during February comparing to November and January (Supplementary Table 1).
Phytoplankton community structure during the spring/summer of 2015/2016 was assessed through Chemical taxonomy (CHEM-TAX) software, using accessory pigments versus in situ Chl-a concentrations measured via high-performance liquid chromatography (HPLC; see "Methods"). The main phytoplankton group over the season were diatoms, followed by haptophytes (Phaeocystis antarctica), cryptophytes, and dinoflagellates, according to the succession stage (Fig. 3a). Diatoms dominated the phytoplankton community composition in relation to the other groups along the whole in situ sampling period, although their relative biomass (to the total in situ Chl-a) was lower in some stations compared to others in different moments during spring/summer (Fig. 3a). To assess the degree to which the water column structure was a primary driver for development and intensity of diatom growth 3,24 , the mixed layer depth (MLD) and water column stability were calculated as a function of seawater potential density (see "Methods"). There was an inverse polynomial relationship between MLD and mean upper ocean stability (averaged over 5−150 m depth; hereafter referred to as upper ocean stability) (Fig. 3b). The significant positive exponential relationship between the upper ocean stability and diatom absolute concentrations (in situ Chl-a) demonstrates that stability, associated with MLD, was an important driver of diatom dynamics (Fig. 3b). This elucidates the increase in biological production during summer months of 2016, when upper ocean physical structures (MLD and stability) were sufficiently shallow and stable to produce the high phytoplankton biomass (in situ Chl-a) registered here. However, as MLD and stability showed similar values between summer months b Exponential curve (R 2 = 0.57; p < 0.0001; n = 25) between upper ocean stability (Stability) and diatom absolute concentrations on surface. Symbol size indicates phytoplankton size class (pico-, nano-or micro-phytoplankton) that dominated with values >40% the community composition proportion in respect to the total Chl-a (considering the three fractionated size classes). Symbol color indicates the sampling month in respect to November (brown), January (gray), and February (black). The inset shows the polynomial inverse relationship (R 2 = 0.51; p < 0.001; n = 26) between mixed layer depth (MLD) and Stability and indicates that a shallow MLD is not necessarily stable. The subregions are indicated by the underscore gray line below the stations' number (a). The letter N, J, or F (November, January, and February, respectively) denotes the respective month of data sampling according to locations of fixed number stations represented by the map on the lower right corner (c). In (a, b), asterisks indicate decreased diatom proportion associated with increased dinoflagellate proportion during distinct moments of summer months (January and February).
(Supplementary Table 1), only the upper ocean physical structures cannot be accounted for the high differences of in situ Chl-a values observed between diatom blooms in January (maximum of 12 mg m -3 ) and February (maximum of 40 mg m -3 ). Likewise, also not explaining these differences of in situ Chl-a values between summer months, macronutrients were highly abundant throughout the seasonal phytoplankton succession (Supplementary Table 1). Furthermore, although no measurements of dissolved iron, which can be considered as a limiting factor to primary productivity 30 , were carried out here, the Antarctic Peninsula continental shelves have been depicted as a substantial source of this micronutrient to the upper ocean, not limiting phytoplankton growth even during intense blooms 31,32 .
The stability was also seen to influence the dominant phytoplankton size class. While micro-phytoplankton dominated over the whole range of the upper ocean stability, smaller size classes showed a more restricted distribution (Fig. 3b). These differences on biomass and dominant phytoplankton size class could be due to distinct moments of the succession stages linked to interspecific differences. In general, the nano-and pico-phytoplankton were more abundant when: (i) stability was lower in the late spring/early summer (mainly small diatoms and small flagellates) ( Fig. 3b and Supplementary Fig. 2), and when (ii) the large micro-phytoplankton cells left the system in a fast-sinking process, leading to a shift to small dinoflagellates (mainly Gymnodiniales < 20 μm) in the late summer ( Fig. 3b) 24 .
To further examine these phytoplankton composition and biomass differences along the seasonal succession, the physiological condition and grazing pressure were also assessed through indices, using accessory pigments and Chl-a degradation products via HPLC analysis (see "Methods"). Predominantly, the highest proportions of photosynthetic carotenoids were observed in the mid/late summer, indicating high photosynthetic activity from the phytoplankton cells associated with diatom blooms. During the late spring, however, phytoplankton communities accumulated more photoprotective carotenoid pigments (see Supplementary Fig. 3a, b), likely to mitigate the cells photooxidative damage caused by enhanced irradiation, which favors species adapted to light stress 33 . This indicates that the large environmental fluctuations characteristic of the spring, due to the physiological stress, typically favor a more heterogeneous community dominated by nanoflagellates and other smaller phytoplankton cells (Supplementary Fig. 2 and Fig. 3a, b) 26,34 . The grazing and senescence indices, in turn, indicate that top-down control and phytoplankton cell deterioration, respectively, were more intense during January, compared to February, which may explain part of the monthly in situ Chl-a differences (Supplementary Fig. 3c, d).
The top-down control performing on phytoplankton biomass may be exerted by both microheterotrophs (e.g., bacteria, heterotrophic flagellates, and ciliates) and metazoan plankton, such as krill and salps. The fecal pellets of metazoans tend to sink quickly, and hence are improbable to be sampled 35 . However, it is unlikely that in situ Chl-a could have reached the high values recorded here in January and February 2016 under a grazing pressure of metazoans 30,36 . Therefore, it is expected that a grazing pressure associated with microheterotrophs was the major topdown control in this situation. Microheterotrophs differ from metazoans in that they are included together with phytoplankton during sample filtration for HPLC/CHEMTAX analysis and, consequently, pigments associated with their grazing, ingested and/or digested, are extracted and analyzed 35 . Moreover, their small fecal pellets sink slowly and are thus likely to be captured with the sample 35 . To support this case, converging with our grazing index data ( Supplementary Fig. 3c), diatom blooms generally characterized by large centric cells were clearly different between both summer periods. While in February the diatom bloom was mainly composed by the large centric Odontella weissflogii (mostly > 70 µm in length; ref. 24 ), during January a large number (> 2.5 × 10 6 cells L -1 ) of small (< 20 µm in length) pennate Fragilariopsis spp. contributed substantially to the total of diatom bloom biomass. Given that large cell size of phytoplankton can reduce, or prevent, efficiency of zooplankton grazing 37,38 , this intra-group difference between the diatom blooms may be indicating that a selective activity of microheterotrophic grazing occurred on smaller pennate diatoms during January 2016.
Furthermore, the effects of a strong positive SAM trend on sea ice dynamics over the previous years may have caused a severe impact on the survival of larval krill 2 , which was reflected in low krill recruit density during 2016 along the NAP 39 . This further supports a low top-down control from the krill population during this summer, which is also consistent with the grazing index values that were much lower than expected ( Supplementary  Fig. 3c, d) 24 . It has been recorded along the Antarctic regions that these grazing index values are greater than 10% when associated with diatom blooms 40,41 , thus indicating a prevalence of active grazing pressure. Therefore, it is likely that the low top-down control on produced biomass further contributed to the high anomalous values of in situ Chl-a registered over the seasons; even though a grazing pressure associated with microheterotrophic zooplankton assemblages may have driven differences in Chl-a concentrations between the mid/late summer ( Supplementary Fig. 3c, d).
Despite the clear dominance of diatoms throughout the in situ sampling periods in the 2015/2016 spring/summer, there was a notable separation between November and January/February in respect to the biomass contribution (as a proportion to the total Chl-a) of the main phytoplankton functional groups (Supplementary Fig. 4). Using a canonical correspondence analysis (see "Methods"), the diatoms and dinoflagellates were positively associated with temperature, stability and Chl-a, and negatively associated with salinity, MLD and macronutrients to the mid/late summer months. The opposite was observed with the other phytoplankton groups (mostly small nanoflagellates), mainly associated with the stations sampled during November (Supplementary Fig. 4), as discussed above. This indicates that, during summer months, dinoflagellates (mainly Gymnodiniales < 20 μm) were the second most representative taxonomic group, after diatoms.
To determine the statistical relationship between diatoms and dinoflagellates along the NAP seasonal succession, all in situ data collected during February 2016 (ref. 24 ), along with the dataset presented here, were further assessed. In general, there was a positive relationship between absolute concentrations of both groups, suggesting their coexistence (Supplementary Fig. 5a). It is possible to observe that there was a decline in diatom biomass under higher values of upper ocean stability, mainly associated with waters influenced by ice melt coming from the Bellingshausen Sea (Supplementary Fig. 5a) 24 . This was linked to the sedimentation process of larger diatom cells out of the upper layers, due to water column structure dynamics 24 . As a result, following the diatom sinking process, the significant positive relationship between upper ocean stability and dinoflagellates further indicates that these small flagellates were a major contributor for the phytoplankton community composition in the late summer 2016 (Supplementary Figs. 5b and 6). In recent years, it has been observed an increasingly contribution of small dinoflagellates (especially Gymnodiniales < 20 μm) to phytoplankton community composition along the NAP waters during the mid/late summer months [42][43][44] , drawing attention to potential poleward shifts of Antarctic phytoplankton species 12 . These small flagellates have been reported apparently occupying sites/conditions less favorable to diatoms and cryptophytes 24,42 . Furthermore, going into autumn (from March onwards), the decline in SST may still suggest a period of increased wind speed, and thus vertical mixing, which indicates a potential deepening of MLD along the region. This event, together with progressive reduced daylight length 45 , could further explain the decline in Chl-a values from the satellite data to the end of the seasonal phytoplankton succession (Fig. 2), potentially favouring a community dominated by small mixed flagellates and small diatoms (Supplementary Fig. 6) 45 .
Kruskal−Wallis tests were performed (see "Methods"), comparing our February time series data (2008-2017) with the atypical situation of February 2016 (our focus here), to examine differences on water physical parameters and on indices of physiological state and grazing pressure associated with high (> 50%) and low (< 30%) diatom proportions (to the total in situ Chl-a). In general, high diatom proportion (representing values in the time series and in 2016 higher than 50%) showed a strong association with deeper MLD (> 20 m), higher surface salinity (majority > 34), and lower upper ocean stability values (< 4 × 10 -6 m -1 ) relative to low diatom proportion (representing values in the time series lower than 30%; see Fig. 4). These results demonstrate a preferential fundamental niche to diatoms along the NAP (Fig. 4), although there was a marked difference in biomass accumulation during February 2016 (see Supplementary  Fig. 7). In addition, there were also strong differences on indices of physiological state and grazing pressure between high and low diatom proportions. The low proportion of diatoms was associated with higher (lower) ratios of photoprotective (photosynthetic) carotenoids and lower (higher) grazing (senescence) index values, which suggest a poorly productive phytoplankton community relative to high diatom proportion conditions (Supplementary Fig. 8).
The only physical parameter showing a difference between both high diatom proportions was surface temperature, which was colder in 2016. However, this parameter solely cannot explain the much higher diatom biomass observed in the late summer 2016 (mean Chl-a value of 12.4 mg m -3 ) relative to the other years (mean Chl-a values lower than 1 mg m -3 ; see Supplementary Fig. 7). The ratios of photoprotective carotenoids and senescence index showed higher values in the other years (2008−2015, 2017) compared with 2016, pointing for a less favorable condition during these years, even though no statistically significant difference between the ratios of photosynthetic carotenoids and grazing index was found (see Supplementary  Fig. 8). As high diatom proportions showed similar behavior regarding water physical parameters (Fig. 4), differences on both photoprotective and senescence indices may be indicating that less favorable conditions 35,[40][41][42] were preventing diatom biomass accumulation in the other years compared with 2016. However, only through the indices is not possible to disentangle the potential critical factor associated with these unfavorable conditions. Analyses of diatom life cycle, considering interspecific differences, prevailing irradiance, and nutrient availability (mainly those in trace amounts) will be needed to evaluate, and understand, the critical factor (factors) which is (are) hampering diatom bloom formation, even under favorable upper ocean physical structures. The difference in space and time should be considered, assuming the likelihood of sedimentation process (or consumption by metazoan plankton) of large diatom cells before the sampling, that is, low diatom biomass but favorable upper ocean physical structures. Furthermore, although there was no statistical difference between grazing index values associated with high diatom proportions; activity of grazing in the other years cannot be disregarded as a potential driver on diatom biomass, as the grazing index in these years in respect to 2016 showed a wide range of values, which in part exceeded 10% (see Supplementary  Fig. 8c).
Over a longer spatial/temporal scale, the effects of climate warming are progressing southward beyond the NAP 1,3 . The progression of these effects is causing impacts on ocean physics and biological compartments. Accordingly, a decrease in diatom proportion has been pointed northward along the mid/southern WAP, concurrent with an increase in cryptophyte proportion 3 . This shift has also been reported along the NAP as a result of glacier meltwater processes during summer [12][13][14] . Due to meltwater inputs, the cryptophyte fundamental niche has been associated along the NAP with mean values of salinity lower than 34 (refs. 13,14 ), which is consistent with other regions in the mid/ southern WAP 3,11,46 . Moreover, It has also been observed that the peak of cryptophyte abundance is related to warmer air temperature (> 2.3°C), contrasting with diatom peak that was found at lower mean air temperature (~1.2°C) 11 . Such evidence links the reported trends of an increase in cryptophyte proportion and warmer ocean and atmospheric temperatures, which play a pivotal role in glacier retreat rates, driving a substantial glacial melting over the past warming decades along the Antarctic Peninsula 8,9,16 .
During the spring/summer of 2015/2016 reported here, diatom blooms were associated with lower mean air temperatures and high mean surface salinity values (< 0.2°C and > 34, respectively; see Fig. 4, Table 1, and Supplementary Table 1), which could be indicating a lower incoming of meltwater 16 . Likewise, during the moderately strong 2009/2010 El Niño event, large diatom blooms along the NAP coastal region were associated with low air temperature, which caused a delay in the melting of surrounding glaciers, leading to average salinity values greater than 34 (ref. 47 ). These findings are consistent with our results, which reinforce that the effects of the extreme 2015/2016 El Niño on atmospheric/ ocean physics were substantial, triggering anomalous diatom blooms during the seasonal phytoplankton succession (Fig. 4). However, even with large diatom blooms during recent El Niño events, the cryptophyte growth has still been reported along the NAP 13,45 . The emergence of this nanoflagellate was associated with localized inputs of glacial meltwater, decoupled from the diatom blooms, supporting the segregation and/or competition between these two phytoplankton functional groups 11,12 . This differentiation in upper ocean physical structures caused by the influence of glacial meltwater input, which leads to shallower, stabilized, and well-lit upper layers, based on results found here, apparently prevents diatom bloom formation along the NAP due to the fast-sinking process, decreasing biomass accumulation ( Fig. 4d and Supplementary Fig. 5a). It is important to note that, in a first moment, diatom contribution increased concomitantly with the upper ocean stability (Fig. 3), decreasing its biomass on contrary under higher stability values (> 4 × 10 -6 m -1 ; Supplementary Fig. 5a). These results converge with the distribution of low diatom proportion associated with higher (lower) upper ocean stability (surface salinity) and shallower MLD (Fig. 4b−d).
Additionally, other factors such as daily, seasonal and latitudinal fluctuations in solar radiation 41,48 , top-down control 30,49 , interspecific differences 33 , and micronutrient limitations 16 within the upper ocean cannot be disregarded for further understanding diatom dynamics along the NAP, and in other Antarctic coastal regions.
The substantial long-term physical changes along the NAP have been linked to strengthening winds, driven by a positive SAM trend and its co-variability with La Niña 1,10,16 . However, this is opposed to the results reported here, in which the effects of the extreme 2015/2016 El Niño notably superimposed on the persistent positive SAM trend. In this study, we demonstrated a close link between large diatom blooms, upper ocean physical structures and sea ice cover during the spring/summer of 2015/ 2016 along the NAP (Figs. 2-4). Moreover, there was also evidence of low top-down control on the diatom biomass (Supplementary Figs. 3c and 8c), which further supports that bottom-up control is a key mechanism of the marine ecosystems along the WAP 50 . As it is predicted by models that positive SAM frequencies will increase for the next 50 years 51 , coupled with strong winds and warmer temperatures, further glacial retreat and sea ice decline 6 should be expected. Consequently, this would provide a more favorable fundamental niche for cryptophytes 3,12 . This hypothetical future scenario would make large diatom blooms typically episodic events, intensifying and/or establishing what it has been observed over the past warming decades across the NAP 1,[12][13][14]47 . Therefore, despite a seasonal phytoplankton succession with massive concentrations of diatoms during the austral spring/summer of 2015/2016; due to continued warming trends, diatom blooms may become more episodic and restricted at spatial/temporal scales with increasingly meltwater inputs along the NAP.

Methods
Cruise design and sampling collection. Data for this study were collected during three seasonal oceanographic cruises conducted on board the RV Almirante Maximiano of the Brazilian Navy during the late spring and mid/late summer months (November, January and February, respectively) between 2015 and 2016 along the NAP. The study area covered the northern sector of the Gerlache Strait, which separates the Anvers and Brabant Islands from the Antarctic Peninsula, and the Bransfield Strait, located between the Southern Shetland Islands and the Peninsula (Fig. 1). The NAP is a diverse region featuring a range of oceanic and coastal ecosystems with complex hydrographic and biogeochemistry dynamics, which are strongly influenced by transition waters from the Weddell and Bellingshausen Seas 12 . Therefore, the NAP sampling grid in this study was split in two subregions: north and south (Fig. 1). The sampling grid consists of nine fixed numbering stations (Fig. 3), which overall were systematically sampled during the three seasonal oceanographic cruises. However, during November due to persistent extended sea ice cover some fixed stations were not sampled (Fig. 1). On the other hand, during February two additional stations (Sts. 10 and 11) were performed to better evaluate the massive diatom bloom observed (Figs. 1 and 3) 24 . In addition, the Brazilian High Latitude Oceanography Group (GOAL) time series data (2008-2017) were also used here. These data were compiled over 10 years of oceanographic cruises conducted on board the polar vessels Ary Rongel A combined Sea-Bird CTD (conductivity-temperature-depth instrument)/ Carrousel 911+system® equipped with 24 five-litre Niskin bottles was used to measure the hydrographic data profiles, such as in situ temperature (°C) and salinity. Seawater discrete samples were collected at regular depths of 5, 25, 50, 75, 100, and 150 m using Niskin bottles to characterize the vertical distribution of phytoplankton communities. Surface seawater samples (5 m) were taken in all CTD stations for both dissolved nutrients and phytoplankton pigments analyses. These seawater samples (0.5-2.5 L) were filtered under low vacuum through GF/F filters, which were subsequently frozen in liquid nitrogen for later HPLC/CHEMTAX analysis. Additionally, size fractionation of phytoplankton biomass (mg m -3 Chl-a) was performed through sequential filtering of surface water samples (1000 mL). Each surface sample for size fractionation was firstly filtered using a polycarbonate filter with pore size of 20 µm (GE Water & Process Technologies, 47 mm). The collected filtrate was subsequently filtered into a polycarbonate filter with pore size of 2 µm (GE Water & Process Technologies, 47 mm). Ultimately, the collected filtrate was filtered into GF/F filter (Whatman) to retain the cells between 0.7 and 2 µm. The first two filtrations were performed by gravity and the last under low vacuum pump (< 5 in. Hg). The filtration by gravity minimizes risks of breaking phytoplankton cells. The filters of size fractionation were also frozen in liquid nitrogen for later HPLC/CHEMTAX analysis. This procedure yielded measurements of concentrations of Chl-a (mg m -3 ) and other pigments, as well as phytoplankton functional groups' contributions for the pico-phytoplankton (0.7-2 µm), nano-phytoplankton (2-20 µm), and micro-phytoplankton (> 20 µm).
Physical measurements. Seawater thermodynamic calculations were performed with the Thermodynamic Equation of Seawater-2010 (TEOS-10) 53 . For each CTD cast, the seawater potential density (ρ , kg m -3 ) referenced to 0 dbar was determined based on seawater conservative temperature, absolute salinity and pressure data. To examine the water column physical structure, MLD (m) and water column stability (×10 -6 m -1 ) were calculated from 5-m smoothed profiles of ρ, due to the high surface variability. Although recent studies have estimated MLD through the depth of the maximum water column buoyancy frequency (N 2 ), this method does not consider the strength of stratification, but only the homogeneity of the upper ocean layer 3,21 . Therefore, to determine the variation of MLD apart of the strength of stratification, MLD was calculated as the depth at which ρ deviates from its 10 m depth value by a threshold of Δρ = 0.03 kg m -3 (ref. 54 ). The water column stability was estimated through the N 2 by the local gravity, using vertical ρ variations 24 .
To represent the variation of stability values estimated to each meter within the upper ocean layers, the mean stability between 5 and 150 m (upper ocean stability) was used based on the analysis of ref. 24 . It is important to note that despite of MLD and upper ocean stability are not completely independent, since both are calculated from the ρ, their statistical relationship is not linear; that is, a shallow MLD is not necessarily stable (Fig. 3b).
Dissolved nutrients analysis. Cellulose acetate membrane filters were used to filtrate the surface seawater samples to determine dissolved inorganic nutrients (i.e., dissolved inorganic nitrogen-DIN: nitrate, nitrite and ammonium; phosphate and silicic acid). Samples were immediately frozen at -80°C until the laboratory analysis. Nutrients were analyzed following the spectrophotometric determination methods described in ref. 55 at the Laboratory of Biogeochemistry at the Federal University of Rio de Janeiro, Brazil. In the form of reactive Si, silicic acid measurements were corrected for sea salt interference. Reaction with ammonium molybdate was used to measure orthophosphate, with absorption readings at 885 nm.
HPLC/CHEMTAX analysis. In the laboratory, the filters were placed in a screwcap centrifuge tube with 3 mL of 95% cold-buffered methanol (2% ammonium acetate) containing 0.05 mg L -1 trans-β-apo-8ʹ-carotenal (Fluka) as internal standard. Samples were sonicated for 5 min in an ice-water bath and placed at -20°C for 1 h. Then, samples were centrifuged at 1100 × g for 5 min at 3°C. To separate the extract from remains of filter and cell debris, the supernatants were filtered through Fluoropore PTFE membrane filters (0.2 μm pore size). Immediately prior to injection, 1000 µL of sample was mixed with 400 µL of Milli-Q water in 2.0 mL amber glass sample vials, which then were placed in the HPLC cooling rack (4°C). A Shimadzu HPLC constituted by a solvent distributor module (LC-20AD) with a control system (CBM-20A), a photodiode detector (SPDM20A), and a fluorescence detector (RF-10AXL) was used to analyze the pigment extracts. The chromatographic separation of the pigments was performed using a monomeric C8 column (SunFire; 15 cm long; 4.6 mm in diameter; 3.5 μm particle size) at a constant temperature of 25°C. The mobile phase (solvent) and respective gradient followed the method developed by ref. 56 , which was discussed and optimized 57 , with a flow rate of 1 mL min −1 , injection volume of 100 μL, and 40 min runs. The absorbance spectra and retention times were used to identify the pigments. The concentrations of the pigments were calculated from the signals in the photodiode array detector in comparison with commercial standards obtained from DHI (Institute for Water and Environment, Denmark). LC-Solution software was used to integrate the peaks that were checked manually and corrected when necessary. In order to correct for losses and volume changes, the concentrations of the pigments were normalized to the internal standard. To reduce the uncertainty of pigments found in low concentrations 58 , a quality assurance threshold procedure through application of quantification limit (LOQ) and detection limit (LOD) was applied to the pigment data. The LOQ and LOD procedures were performed according to ref. 57 .
Three types of Chl-a degradation products: chlorophyllide-a (Chlide-a), pheophytin-a (Phytin-a), and pheophorbide-a (Phide-a) were identified and quantified by HPLC analysis. The sum of the Phytin-a and Phide-a was used to calculate the proxy of grazing pressure, while Chlide-a was used to calculate the senescence of phytoplankton cells 40,41,59 . These two indices were divided by the total of degradation products plus Chl-a (sum of Chl-a, Phytin-a, Phide-a and Chlide-a), and represented as a percentage ( Supplementary Fig. 3c, d).
The phytoplankton community faces variable light regimes within only one single day, when mixed up to the surface from depth, experiencing high irradiance on top, which contrasts with low irradiance at the base of MLD 48 . To efficiently maximize ambient light regimes the phytoplankton cells adjust their photosynthetic apparatus, expanding their light-harvesting complexes under low light conditions, or synthesizing photoprotective pigments under high light conditions, which prevent photodamage 60 . Therefore, via HPLC analysis, photopigment indices were used to investigate phytoplankton pigment acclimations in response to different environmental regimes 40,42 . The carotenoids were separated into photosynthetic carotenoids (PSCs) and photoprotective carotenoids (PPCs). The PSC included 19ʹ-butanoyloxyfucoxanthin, 19ʹ-hexanoyloxyfucoxanthin, fucoxanthin, prasinoxanthin and peridinin, whereas the PPC were composed by zeaxanthin, violaxanthin, lutein, alloxanthin, diadinoxanthin, diatoxanthin, β,βcarotene and β,ε-carotene 61 . Thus, two photo-pigment indices were derived and used here, PSC: Chl-a and PPC: Chl-a. These indices show the proportion of each carotenoid pigment group (PSCs and PPCs) to the total biomass (Chl-a). It is important to highlight that phytoplankton pigment composition may be influenced by photoacclimation processes and nutrient availabilities 48,60,62 . Despite phytoplankton responses could take time scales of hours to days, our analyses associated with PSC: Chl-a and PPC: Chl-a include a wide scale of space and time with our time series of 2008−2017 sampled during February months ( Supplementary Fig. 8). Therefore, it is possible to access the state (taking into account the large number of samples in space and time, which allow us obtaining reference values) by which phytoplankton communities were being exposed and subsequently acclimated, thus indicating in part the environmental conditions associated with these cells, under stress or not. The difference in taxonomic composition should also be considered.
The class-specific accessory pigments and total Chl-a were used to calculate the absolute and relative contribution of phytoplankton groups through the Chemical taxonomy (CHEMTAX) software v1.95 (ref. 63 ). CHEMTAX uses a factor analysis and steepest descent algorithm to best fit the data onto an initial matrix of pigment ratios (the ratios between the respective accessory pigments and Chl-a). This software package has been extensively and successfully used in many worldwide investigations, including in the Southern Ocean 3,13,14,24 , to determine the distribution and biomass of phytoplankton functional groups. This approach provides valuable information about the whole phytoplankton community, including small-size species, which are normally difficult to identify by light microscopy. The Chl-a was used in this study as a proxy for biomass 64,65 , once this photosynthetic pigment is common to all autotrophic phytoplankton. Therefore, we use the term Chl-a referring to either absolute biomass or relative biomass attributed to the corresponding taxonomic groups.
Seven algal groups were chosen for CHEMTAX analysis based on initial pigment ratios derived for the NAP (ref. 42 ). The groups resolved include diatoms Type-A, diatoms Type-B, dinoflagellates Type-A, dinoflagellates Type-B, cryptophytes, haptophytes type 8 (Phaeocystis antarctica) and green flagellates. The two types of diatoms (diatoms Type-A and diatoms Type-B) were grouped as a same chemotaxonomic group: the diatoms. The two types of dinoflagellates (dinoflagellates Type-A and dinoflagellates Type-B) were also grouped as a same chemotaxonomic group: the dinoflagellates. It is important to note that diatoms Type-A were not considered in the analysis of CHEMTAX during February 2016, due to the absence of chlorophyll c 1 (ref. 24 ). Similarly, during November 2015, the diatoms Type-A and dinoflagellates Type-B were not considered in the analysis of CHEMTAX, due to the absence of chlorophyll c 1 and gyroxanthin esters, respectively.
The pigment dataset was separated in two distinct subregions: north and south (Fig. 1). These subregion datasets were further split into three bins according to sample depth (0-25 m, 25-50 m and > 50 m), each of which was separately analyzed by CHEMTAX to allow for variation in pigment: Chl-a ratios with depth 66 . Regarding the fractionation data, an additional split of both pigment datasets was performed according to size class (> 20 µm, 2−20 µm, 0.7−2 µm). This procedure assures a homogeneity of the pigment: Chl-a ratios, considering all samples from the same bin, and provides some compensation for changes in these ratios, which are known to vary with light availability 67,68 . Series of 60 pigment ratio matrices were generated for optimization of the input matrix by multiplying each ratio from the initial matrix by a random function as described in ref. 69 . The average of the best six (10%) output matrices (with the lowest residual or mean square root error) was taken as the optimized result to each bin. Lastly, the results of phytoplankton chemotaxonomic groups derived from the CHEMTAX were quality-controlled by microscopic analysis. . In order to integrate measurements from each sensor, data are band-shifted to match the SeaWiFS bands and subsequently biascorrected. The OC-CCI product is validated with extensive in-situ datasets, exhibiting high correlation at a global level (R 2 = 0.73) 71 .
2010−2019 Satellite SST (°C) and sea ice concentration (%) were extracted from the ESA SST CCI and C3S global Sea Surface Temperature Reprocessed product (available at CMEMS) 72 . This product provides gap-free daily-averaged data with a spatial resolution of 5 × 5 km at a global level. SST data were merged from measurements obtained from the Along Track Scanning Radiometer (ATSR), Sea and Land Surface Temperature Radiometer (SLSTR) and the Advanced Very High Resolution Radiometer (AVHRR) sensors. Sea ice concentration data were obtained from the EUMETSAT Ocean and Sea Ice Satellite Application Facility (OSI-SAF) OSI-450 and OSI-430. Sea ice duration (days) was defined based on the analysis of ref. 10 : as an interval between the sea ice advance in autumn and sea ice retreat in spring. Day of advance is identified when sea ice concentration first exceeds 15% (i.e., the approximate ice edge) for at least 5 days, whereas the day of retreat is identified when sea ice concentration remains below 15% until the end of the study period. For all measured variables (except sea ice duration), including air temperature (recorded at local meteorological stations), daily means were calculated for all pixels inside each subregion (north and south; Fig. 1) from the early November to the end of February (defined as spring/summer here) for a period of 10 years (2010−2019). The November−February mean was determined to the austral spring/summer of 2015/2016 and to the 10 years of the time series, which includes every austral spring/summer into this period, characterizing the decadal average (Figs. 1, 2 and Table 1). The estimated sea ice duration, as a difference between the advance and retreat of sea ice from the satellite data, was calculated to each year over the 10 years. Then, sea ice duration for the decadal average represents an average considering all years, whereas the sea ice duration for 2015/2016 was estimated as a sole value between advance (April) and retreat (September) of sea ice of 2015, which reflected during the spring/summer of 2015/2016. These averages were used to assess the intensity of the effects associated with the extreme 2015/ 2016 El Niño 22,25 , with respect to the decadal pattern of the analyzed variables along the NAP, which were mainly influenced by the persistent positive SAM trend 2,16 .
Statistical analysis. Relationships between biomass (in situ Chl-a) of phytoplankton groups and environmental variables at surface (except stability, which was defined as the mean upper ocean stability between 5 and 150 m, and MLD) were determined through canonical correspondence analysis (CCA) using CANOCO for Windows 4.5 software 73 . Biotic variables were represented by the CHEMTAXderived taxonomic groups' relative biomass (as a proportion to the total Chl-a). Environmental variables included upper ocean stability (Stability), MLD, sea surface temperature (Temperature), sea surface salinity (Salinity), dissolved inorganic nitrogen (DIN), phosphate (Phosphate) and silicic acid (Silicic acid). All variables were log-transformed before analysis to reduce the influence of different scales in the datasets. In order to evaluate the significance of the CCA, Monte-Carlo tests were run based on 499 permutations under a reduced model (p = 0.002).
For north and south subregions, differences in the spring/summer average values from the satellite data of Chl-a, SST, air temperature, and sea ice concentration between the decadal-mean and 2015/2016 were analyzed using Student's t test comparisons (Table 1). This statistical analysis was not performed for sea ice duration because this parameter is presented as a sole value, which prevents the use of Student's t test comparisons. Kruskal−Wallis tests were performed, comparing our February time series data (2008-2017) with the atypical situation of February 2016 (our focus here) to examine differences on water physical parameters, physiological state (indices) and grazing index values associated with high (> 50%) and low (< 30%) diatom proportions. This is a nonparametric test. February data (from 2008 to 2015 and 2017) were split in two bins regarding diatom proportion to the total biomass (in situ Chl-a): > 50% and < 30%. Then, these two bins, along with diatom proportion higher than 50% in February 2016, were statistically tested to access potential differences of diatom distribution in respect to environmental settings ( Fig. 4 and Supplementary Fig. 8).