Sampling was conducted at a total of 12 sites in the northern Yellow Sea in August 2011-2020 ( Figure 1 ). The sampling sites varied each summer, with 4 to 10 sites per year except for 2019, when there was only one sampling site ( Supplementary Table S1 ). Taking the 10°C isotherm as the boundary for the YSCWM (Guan, 1963; Zhao, 1985), the sites with bottom water temperatures at or below 10°C were regarded as within the Cold Water Mass Region (CWMR), while the others fell into the Outer Region (OR). Site YS6 was notable because its bottom temperatures occasionally exceeded 10°C during the study period, suggesting it is located on the edge of the YSCWM ( Supplementary Table S1 ). However, since its bottom temperatures were at or below 10°C in most years, it was classified within the CWMR group. According to Lyu et al. (2022), the euphotic layer depth of the entire study area ranged from 20 to 30 m in August. Therefore, the water depths of sites YS4, YS9, and YS10 were also within this range ( Supplementary Table S1 ), indicating that their bottom layers were approximately at the base of the euphotic layer. In contrast, the water depths of the other nine sites exceeded 40 m, indicating that their bottom layers were below the base of the euphotic layer.

The collections of in-situ physical parameters, biological samples, and chemical samples were carried out simultaneously on the same voyage every summer. Temperature, salinity, depth, and pH were determined with a SeaBird CTD system (SeaBird, USA), while water samples from various depths were collected using multi-bottled water samplers. Dissolved oxygen (DO) was determined by iodometric titration according to the National Standard of China (GB/T 17378.4-2007). A subsample (500 mL) from each depth layer was filtered through a 0.7 μm GF/F filter (Whatman, UK) using a vacuum of no more than 50 kPa for subsequent analysis of chlorophyll a (Chl-a). The filtered seawater was transferred to a 100 mL polyethylene terephthalate bottle and frozen at -20°C for subsequent nutrient analysis. Phytoplankton samples were collected by towing a plankton net (76 μm mesh size, 37 cm mouth diameter) vertically from a depth of 2 m above the seafloor to the surface at a speed of 0.5 m/s. The volume of filtered water was calculated by multiplying the rope length by the mouth area of the net (Wang et al., 2024). These samples were then preserved in 1 L polyethylene bottles with formaldehyde (final concentration 5%) for subsequent analysis, specifically targeting micro-phytoplankton species. Zooplankton samples were collected using two different mesh sizes of plankton nets (net I: 505 μm mesh size, 50 cm mouth diameter; net II: 160 μm mesh size, 32 cm mouth diameter), fixed with formaldehyde (final concentration 5%), and stored in 500 mL polyethylene bottles. The volume of filtered water was calculated in the same way as for phytoplankton. Sampling of plankton was performed according to the National Standard of China (GB/T 12763.6-2007).

Chl-a was extracted in 5 mL of 90% acetone in the dark for 24 h at 4°C, and a fluorometer (Turner Designs, USA) was used for measuring fluorescence (Welschmeyer, 1994). Dissolved inorganic nitrogen (DIN: NO₃ ^–N, NO₂ ^–N, and NH₄ ⁺-N), dissolved inorganic phosphorus (DIP: PO₄ ^3–-P), and dissolved silicate (DSi: SiO₃ ^2–-Si) were measured using a continuous flow analyzer (AA3, Seal Analytical, Germany). The element (N/P, N/Si, and Si/P) ratios were the ratios of the molar concentrations of elements. Phytoplankton samples were identified and counted using an Olympus CX31 microscope (Olympus, Tokyo, Japan) following the Utermöhl method (Utermöhl, 1958). Each species was identified to the lowest possible taxonomic level and the validity of taxonomic names was checked on the AlgaeBase website (http://www.algaebase.org). Zooplankton samples were identified and counted under a dissecting microscope (Leica, Germany). The phytoplankton or zooplankton abundance was obtained by dividing the number of phytoplankton cells or zooplankton individuals by the volume of filtered water.

The vertical mean temperature gradient (hereinafter referred to as gradient) was calculated by Equation 1 to characterize the thermocline intensity (Zhou et al., 2009):

where T _b is bottom temperature; T _s is surface temperature; D_max is maximum depth, i.e., water column depth.

The depth-averaged value (V _d) of the physicochemical factors was calculated using Equation 2:

where D_max is as defined above; n is the maximum number of sampling layers; V_i is the value in the ith sampling layer; D_i is the depth of the ith water layer.

The average Niño3.4 index (https://psl.noaa.gov/data/climateindices/list/) for the previous winter (December-February) was used as an index for the El Niño-Southern Oscillation (ENSO).

Nutrient limitations were determined following the minimum thresholds for phytoplankton growth (also known as absolute limitation) and stoichiometry ratios (also known as potential limitation), respectively (Yin et al., 2013). DIN concentrations of <1 μM, DIP concentrations of <0.1 μM, and DSi concentrations of <2 μM were considered limiting to phytoplankton growth. Stoichiometric limitations were determined with atomic ratios as follows: N/P<10 and N/Si<1 for nitrogen (N) limitation; N/P>22 and Si/P>22 for phosphorus (P) limitation; and N/Si>1 and Si/P<10 for silicon (Si) limitation. Due to the absence of DSi data, the summers of 2011-2013 were exceptional, with N/P<10 for N limitation and N/P>22 for P limitation (Song et al., 2020).

The diversity indices, including Margalef’s richness (D) (Margalef, 1968), Shannon-Weiner diversity (H′) (Shannon and Weaver, 1949), and Pielou’s evenness (J) (Pielou, 1969), were calculated for each sample using the vegan package according to Equations 3–5.

where S is total number of species in the sample; N is total number of individuals in the sample; n_i is abundance of the ith species in the sample.

The dominance (Yi) of phytoplankton species/genus for each summer was calculated according to Equation 6, and species/genera with Yi ≥ 0.02 were considered dominant in the studied area (Sun et al., 2004).

where n_i and N are the same as above; f_i is ratio of the number of sites where the ith species appear to the total number of sites.

The Dia/Dino index for each sample was calculated by Equation 7 (Mokrane et al., 2019):

where Ab _Dia is diatom abundance; Ab _Dino is dinoflagellate abundance.

Locally weighted least squares (Loess) regression was conducted with R (R Core Team, 2018) to determine the dynamic trends of environmental variables as well as the relationship between the phytoplankton and environmental variables. Multiple comparisons of phytoplankton abundance and diversity indices among years were conducted using R. Shapiro-Wilk and Levene’s tests were performed using the stats and car packages respectively to assess the normality and homogeneity of variance. Kruskal-Wallis test with posthoc Nemenyi test was employed using the stats package and PMCMRplus packages since the data were non-normally distributed. Mann-Whitney U test was performed to compare environmental factors and dominant species between different regions. Statistical significance was set at p < 0.05 for hypothesis testing of differences.

The depth-averaged environmental variables were used to investigate the correlation between phytoplankton and environmental factors. Their relationship was also explored using constrained ordinations, based on square-root transformed abundance data and log₁₀(x+1) transformed environmental data except for pH data. The missing values of environmental factors were replaced by the overall means. Firstly, detrended correspondence analysis (DCA) was conducted on the biological data and the type of canonical ordination methods was selected according to the length of the first DCA axis. The unimodal model (canonical correspondence analysis, CCA) was used when the length of the first DCA axis length was > 4, and the linear model (redundancy analysis, RDA) was used when the length was < 3 (Lepš and Šmilauer, 2003). Then, variables with variance inflation factors (VIF) > 10 were removed from the ordination analysis to avoid high collinearity. Forward selection was carried out using the packfor package to identify statistically significant explanatory variables. Spearman correlations between phytoplankton taxa and diversity indices and between phytoplankton taxa/diversity and environmental factors were calculated with the Hmisc package. Spearman correlation analysis was also conducted between the site-averaged abundance in August and the Niño3.4 index for the previous winter, in order to explore the relationship between climate change and phytoplankton dynamics.

During the study period, several environmental factors in the NYS during summer showed clear temporal trends. The depth-averaged temperature in 2017 was 2.8°C higher than in 2011. Numerous surveys indicated a long-term rising trend in temperature within the NYS (Li et al., 2022; Park et al., 2015; Pei et al., 2017). However, the upward trend observed during the study period was even more pronounced, aligning with the characteristics noted by Li et al. (2022) for the period from 2011 to 2017. The temperature gradient also exhibited a general upward trend. In the CWMR, the gradient changed little due to the simultaneous rise in both surface and bottom temperatures. In contrast, the gradient in the OR tended to rise as the temperature difference between the surface and bottom increased. Additionally, salinity showed an obvious upward trend during the study period, continuing the long-term increase in salinity observed from the mid-20^th century to the early 21^st century (Ma et al., 2006; Lv, 2008). Previous studies attributed this long-term increase in salinity to the growing disparity between evaporation and precipitation, as well as the reduction of freshwater influx into the NYS (Ma et al., 2006; Lv, 2008). Collectively, the long-term trends in temperature, gradient, and salinity in the NYS were related to global climate change.

The nutrient structure in the NYS also underwent drastic changes during the study period. Notably, the DIN concentration and N/P ratio peaked in the summer of 2011, followed by a marked decline and stabilization in the following years. This shift was primarily due to stringent anthropogenic nitrogen emission regulations implemented by the Chinese government in the early 2010s, which led to a reduction of DIN levels in the NYS, which was predominantly sourced from river inputs and rainfall, thereby reversing a previously continuous upward trend (Yu et al., 2019; Zheng and Zhai, 2021). The nutrient structure in the NYS, particularly the N/P ratio, decreased along with the change in DIN level. In addition, the Chl-a concentration, a proxy for phytoplankton biomass, exhibited a significant decline after reaching its highest levels in the summers of 2011 and 2013, ultimately hitting a low in the summer of 2017. Given that the temporal trends in Chl-a concentration mirrored those of DIN and the N/P ratio, it was presumed that these changes in nutrient levels and structure were primary drivers of the observed variations in Chl-a. A recent study have indicated that enhanced summer stratification in the NYS, which restricted the upward turbulent diffusion of deep waters rich in nutrients and Chl-a, also contributed to the decline in Chl-a concentration (Zhai et al., 2023). Additionally, a downward trend in pH was observed in the summers of 2011 to 2017, followed by a rebound in the summers after 2017. Li et al. (2022) reported a significant decrease in pH in the NYS over a 40-year period (1976-2017), particularly during summer, attributing this decline to increased net respiration in biomes due to eutrophication. In this study, pH exhibited a positive correlation with Chl-a, showing a marked decrease initially, followed by a slight increase after 2017. Therefore, the overall decrease in pH was less likely to be linked to an increase in total phytoplankton biomass and community net respiration. However, the abundances of net I zooplankton and certain dominant phytoplankton species, such as Paralia sulcata and Tripos longipes, were significantly and negatively correlated with pH. This suggested that a decline in pH might be associated with an increase in the biomass of these dominant species and macro/mesozooplankton, and vice versa.

The abundances of the total phytoplankton, Bacillariophyta and Dinoflagellata, the Dia/Dino index, and the diversity indices in the NYS showed no significant increasing or decreasing trend during the study period. Bacillariophyta and Dinoflagellata were the main phytoplankton phyla, and their relative dominances varied greatly, with Bacillariophyta dominating in the summers of 2011 and 2015, and Dinoflagellata dominating in 2014 and 2019. This result was consistent with previous surveys that the dominant groups were not consistent in the summertime of different years (Zhang et al., 2014; Hou et al., 2021; Fu, 2021). Additionally, some previous surveys based on multilayer sampling have indicated that the vertical distribution patterns of dominant groups also varied among different years. Zhang et al. (2014) reported that Dinoflagellata dominated in the phytoplankton communities in the NYS in June 2011, with Prorocentrum minimum being the primary dominant species; Dinoflagellata, similar to the total phytoplankton, exhibited a higher abundance in the upper waters (20m depth), whereas Bacillariophyta, primarily consisting of P. sulcata, showed a higher abundance in the bottom waters. Hou et al. (2021) found that in June and July of 2013, Bacillariophyta dominated in the phytoplankton communities in the NYS and the phytoplankton abundance was higher in the middle and bottom layers. The most important dominant species were P. minimum, Cylindrotheca closterium, P. sulcata and Thalassiosira sp. According to the observations of Lv et al. (2016), Bacillariophyta was dominant in phytoplankton communities in the cold water mass region in the NYS in August 2014, and the main dominant species were Thalassiosira sp. and Thalassiosira pacifica. High values of the total phytoplankton abundance and Chl-a concentration occurred from the subsurface (20m depth) to the bottom. Fu (2021) reported that Bacillariophyta was the dominant group in the NYS in July and August of 2020, and P. sulcata was the first dominant species. The abundances of total phytoplankton and Bacillariophyta did not change significantly with depth, but the abundances of Dinoflagellata were higher at the surface and bottom than in the middle layer.

Whether diatoms or dinoflagellates dominated the phytoplankton community often depends on the interaction of multiple environmental factors and is of great concern in the context of climate change and increased eutrophication (Bi et al., 2021). Therefore, we tried to explain what factors controlled the relative dominance of diatoms and dinoflagellates in the phytoplankton communities in the NYS. Our research showed that the DIN concentration and N/P ratio were the highest in the summer of 2011, and the Dia/Dino index was also the highest in 2011. This result was inconsistent with most studies, which indicate that an increase in the N/P ratio tends to increase the dominance of dinoflagellates (Glibert et al., 2011; Xiao et al., 2018). The high diatom abundance in 2011 was mainly due to the high diatom (primarily Chaetoceros pseudocurvisetus) abundance at the nearshore site YS9, and similar situations also occurred at the site YS9 in the summers of 2015 and 2017 (primarily Chaetoceros curvisetus and Proboscia indica, respectively). The DIN level in the surface layer of site YS9 was higher in 2011 (14.47 μM), and the N/P ratio was higher in 2011 and 2015 (2011: 213.63, 2015: 94.97). The elevated concentrations of DIN in this region were primarily attributed to inputs from coastal rivers, which facilitated the prolific growth and reproduction of diatoms. We hypothesized that the elevated levels of nutrients imported by coastal rivers were responsible for the extensive growth and proliferation of diatoms in this coastal region. The high N/P ratio can be interpreted as a consequence of phosphorus depletion due to diatom utilization, while nitrogen levels persisted at relatively high concentrations. In contrast, the high diatom abundance at this site in 2017 was likely attributable to increased temperatures and stratification of the water column. This was because the surface temperature at site YS9 was higher and the gradient was stronger than in other years, and P. indica occurred in high abundance at site YS9 with a positive correlation with temperature. These results were consistent with a previous study showing that P. indica has adaptive advantages in high-temperature and stratified environments (Gómez and Souissi, 2007). In summary, our results suggest that the variation in the relative dominance of diatoms and dinoflagellates across different years may be jointly controlled by multiple factors, including nutrient levels and the intensity of water column stratification.

The hydrological environment of the NYS is mainly influenced by runoff, coastal currents, seasonal water masses and the Yellow Sea Warm Current, with varying degrees of impact during different periods of the year (Bi et al., 2019; Zhang et al., 2021). The thermocline in the NYS begins to form in spring and reaches its peak strength in summer (August). The YSCWM is the most prominent hydrological phenomenon in the NYS during summer, with its formation, development and weakening nearly coinciding with those of the thermocline. Our data showed that the temperature difference between the surface and bottom waters in the cold center of the YSCWM reached up to 20°C. To assess the effects of the YSCWM on phytoplankton communities, we compared environmental factors and biological indicators between the CWMR and OR.

In terms of environmental features, the CWMR exhibited a stronger gradient compared to the OR. The bottom waters in the CWMR had lower temperatures, pH, N/P and Si/P ratios, and Chl-a levels, while exhibiting higher salinity, DO, and nutrient (DIN, PO₄ ^3–-P, and SiO₃ ^2–-Si) concentrations. Regarding biological indices, the total phytoplankton abundance, the absolute/relative abundances of Bacillariophyta, the Dia/Dino index, and species richness (S) were all significantly lower in the CWMR than in the OR. The lower phytoplankton abundance in the CWMR was primarily due to nutrient limitations in the upper waters caused by stratification, which limited the growth of certain phytoplankton taxa. In the CWMR, P limitation was most pronounced in the upper waters, while PO₄ ^3–-P was generally sufficient in the deep waters. However, the strong stratification prevented the upward transport of nutrients from the deep waters to the surface. The significantly lower diatom abundance and Dia/Dino index in the CWMR could also be largely attributed to the stronger stratification and the distinct nutrient distribution patterns in this region. The vertical migration of dinoflagellates enabled them to access both the deep nutrient pool and near-surface light at different times of the day, giving them a competitive advantage over other phytoplankton groups (Fu et al., 2018; Ji and Franks, 2007; Glibert et al., 2013; Wang et al., 2014). Additionally, the significantly lower phytoplankton species richness (S) in the CWMR was influenced by nutrient distribution patterns caused by stratification, as well as lower temperatures, since phytoplankton diversity is suggested to be positively correlated with temperature within a certain range (Righetti et al., 2019; Segura et al., 2015).

During the study period, the absolute and potential P limitations in the NYS were both significant although the relative P limitation was weakened in the summers of 2017-2020. In comparison, the N and Si limitations were less pronounced than the P limitation. Spearman analysis revealed that the Dia/Dino index was significantly associated with temperature, gradient and PO₄ ^3–-P concentration. The absolute abundance of Bacillariophyta showed a negative correlation with DIN and SiO₃ ^2–-Si, while the relative abundance of Bacillariophyta showed a negative correlation with PO₄ ^3–-P. This was because higher depth-averaged DIN, PO₄ ^3–-P and SiO₃ ^2–-Si levels were mainly observed in the CWMR, where the abundances of total phytoplankton and Bacillariophyta were low despite high nutrient levels in the bottom. In contrast to Bacillariophyta, the absolute abundance of Dinoflagellata did not correlate with nutrients while the relative abundance of Dinoflagellata was positively correlated with PO₄ ^3–-P. This was mainly because the adaptive strategies of dinoflagellates enable them to effectively navigate nutrient limitations in the highly stratified environment in summer. Additionally, RDA indicated that the relative abundance of Bacillariophyta was negatively DO. Such relationship did not indicate a direct impact of DO on phytoplankton groups but was attributable to dominance of dinoflagellates in CWMR compared to the OR and the higher average depth-averaged DO concentration was also higher in the CWMR. The higher DO concentration in the CWMR than the OR in summer was mainly due to the low temperature of the YSCWM, which could increase the solubility of oxygen (Wei et al., 2019).

Many studies have reported that Paralia sulcata dominated the phytoplankton communities in the NYS during the summer months, primarily concentrated at the bottom of the water column (Fu, 2021; Hou et al., 2021; Zhang et al., 2014). In this study, P. sulcata consistently dominated throughout the study period, with the highest abundances recorded in the summers of 2015-2017. The locations of higher P. sulcata abundance were variable, sometimes found in the CWMR (e.g., YS1 in 2017 and YS6 in 2018) and sometimes observed in the OR (e.g., YS8 in 2015 and YS4 in 2016). The absolute abundance of P. sulcata in the CWMR was slightly but not significantly lower than that in OR, while the relative abundances in both regions were nearly equal, which indicated that the distribution of this species was not greatly affected by the cold water mass. P. sulcata’s sustained dominance in the NYS during summer over a long period might be attributed to its broad adaptability to various environmental factors. As a benthic-pelagic species, P. sulcata is well adapted to survive on sediment, a trait which is supported by its high tolerance to low light conditions and high nutrient availability. Although elevated concentrations of nutrients favor the growth of P. sulcata, this species can tolerate slight limited nutrient levels (Gebühr et al., 2009; Gebühr, 2011). The optimal growth temperature for P. sulcata ranges from 10 to 16°C, as reported in Helgoland Roads, but it can also tolerate lower (<5°C) and higher (~20°C) temperatures, which enables it to remain a dominant species in both summer and winter in this region (Gebühr, 2011) as well as in the Yellow Sea (Liu et al., 2015; Liu and Glibert, 2018). Our findings indicated that P. sulcata exhibited lower abundances during the summers of 2014 and 2020, which aligned with the abundances of diatoms and total phytoplankton. Conversely, dinoflagellates showed a higher relative abundance during these periods. This could be attributed to the ecological niche of P. sulcata overlapping with that of certain dinoflagellate species, particularly in phosphorus-deficient waters where P. sulcata is not subjected to severe phosphorus stress (Wang et al., 2014; Yu et al., 2015). Consequently, the dominance of P. sulcata was partially replaced by dinoflagellate species during these times.

The species Tripos muelleri was noteworthy because, in contrast to most dominant species, it exhibited a significantly higher relative abundance in the CWMR compared to the OR. Our results indicated that the absolute abundance of T. muelleri was significantly positively correlated with the average Niño3.4 index for the previous winter as well as the temperature. Current research has a limited understanding of the physiology of this species, but several in-situ studies have demonstrated that T. muelleri and its varieties are widely distributed in temperate and tropical waters, with varying distribution depths depending on the region, predominantly during non-monsoon seasons (Anderson et al., 2022; Chitari et al., 2017; Mikaelyan et al., 2021; Varghese et al., 2022). A recent report of a substantial bloom of T. muelleri in the Gulf of Maine in May attributed its occurrence to a mild winter, reduced spring winds, and warming temperatures in the gulf (Ray, 2023). Our results suggest that this species has a strong adaptability to the highly stratified environment of the CWMR and demonstrates a clear response to ENSO, implying that T. muelleri may be more sensitive to climate change than other species.

Overall, this study demonstrated that the YSCWM had important effects on the phytoplankton community structure in the NYS in summer, with taxa that exhibited strong adaptability to the stratified environment, such as Dinoflagellata and P. sulcata, dominating this area. Although eutrophication and nutrient structure imbalances caused by human activities decreased significantly and phytoplankton diversity remained stable, the environmental changes related to climate change, such as rising water temperatures, salinity, and stratification of the water column, had potential impacts on the stability of phytoplankton community structure. This includes the abundance ratios of diatoms and dinoflagellates, as well as the composition of dominant species. Thus, the effects of climate change on the coastal environment and biome deserve continued attention. Due to the limitations of the sampling method used in this study, namely the vertical trawl, the analysis was based on the phytoplankton abundance in the entire water column rather than at various depths, resulting in the absence of data on the vertical distribution of phytoplankton. We expect to employ a layered sampling method in future studies, which will better demonstrate the spatial patterns of phytoplankton and thus more accurately explain the interactions between phytoplankton and the environment.