INTRODUCTION
The community structure of marine phytoplankton provides essential insights into their functional role within the oceanic carbon cycle, primarily through photosynthesis. As the primary producers in marine ecosystems, phytoplankton contribute 95% of total marine primary production (Sun 2011, Zhou 2014). Their photosynthetic activity generates organic matter that supports marine food webs, forming the foundation of oceanic ecosystems and sustaining broader marine biodiversity (Noman et al. 2019). Due to their short life cycles, phytoplankton are highly sensitive to environmental fluctuations, making shifts in their community composition a valuable indicator of changes in marine ecosystem health, water quality and nutrient dynamics (Jia et al. 2014). Consequently, they are widely recognized as key bioindicators for assessing marine environmental variability and the impacts of climate change (Torrisi and Dell’Uomo 2006). Anthropogenic pressures, such as eutrophication, can drive phytoplankton community shifts toward smaller-sized species, altering nitrogen-to-phosphorus ratio and changing size-structure dynamics, often resulting in increased occurrences of small-size diatom blooms (Chen et al. 2010, Ma 2018). Ocean acidification, characterized by decreasing seawater pH, has been linked to reduction in diatom species diversity (Ma 2018). Additionally, elevated salinity levels significantly inhibit the growth of cryptophytes and diatoms (Van-Meerssche and Pinckney 2017). Global climate-driven environmental changes have further contributed to substantial shifts in phytoplankton community composition, including a transition from diatom-dominated assemblages to communities increasingly dominated by dinoflagellates and cyanobacteria. This shift has been accompanied by a rise in the frequency and intensity of harmful algal blooms (HABs), many of which involve toxigenic species (Anderson et al. 2002, Wu and Kow 2002, Aubry et al. 2012, Liu et al. 2025).
In recent decades, intensified human activities and rapid industrial development have imposed substantial challenges to global environmental stability and climate systems (Hoegh-Guldberg et al. 2018, Sekerci and Ozarslan 2020). Anthropogenic pressures including industrial emissions, fossil fuel combustion, deforestation, and wetland degradation, have exacerbated ocean acidification (Zeebe 2012). Rising global temperatures alter ocean-atmosphere interactions by increasing atmospheric water vapor (Schaffer et al. 2000), which further enhances the greenhouse effect. Additionally, elevated temperatures intensify glacial melting (Najjar et al. 2010) and contribute to sea-level rise, leading to profound transformation in marine ecosystems (Nick et al. 2013). These environmental changes have far-reaching consequences for phytoplankton communities, which form the base of marine food webs. Key consequences include shifts in phytoplankton composition, altering physiological and ecological responses, variations in primary productivity, and disruptions to biogeochemical cycles (Kalhoro et al. 2024, 2025a, 2025b).
The Bohai, Yellow, and East China Sea (BYECS) are ecologically and economically significant marginal seas situated in the northwest Pacific Ocean (Chen 2009). These regions contain some of the world’s most extensive shallow continental shelf systems (Yu et al. 2017). Their environmental dynamics are shaped by strong terrestrial influences and complex hydrodynamic exchange with adjacent oceanic systems. The Bohai Sea (BS) is a semi-enclosed basin influenced by river discharge which delivers high nutrient loads and supports elevated productivity (Wang et al. 2019). The circulation varies seasonally under the combined influence of the East Asian monsoon, tidal mixing, and exchanges with the Yellow Sea (YS) through the Bohai Strait, creating dynamic environmental gradients that shape regional ecosystem (Huang et al. 2012, Zhang et al. 2018, Meng et al. 2020). The East China Sea (ECS), one of the marginal seas of the western north Pacific. The primary circulation in the BYECS is driven by the Taiwan Warm Current, which flows along the Kuroshio-influenced shelf and by nearshore coastal currents (Su and Yuan 2005, Chen 2009, Yoon et al. 2015). Additionally, the large freshwater discharge from the Changjiang (Yangtze) River exerts a profound influence on the ECS hydrography, nutrient distribution, and ecosystem structure (Zhang et al. 2007, Zhao et al. 2019). ECS considered as post productive region with supports fisheries and biological productivity (Furuya et al. 1996). In contrast, the YS is a semi-enclosed shallow shelf sea characterized by distinct hydrodynamic features. Its circulation is driven by a central warm current and peripheral coastal currents. During summer, the central region of the YS is also affected by the presence of cold, stratified seasonal water masses (Zhang et al. 1996). The interplay of these water masses, coupled with significant terrestrial material input, regulates ecological and biogeochemical processes across the BYECS (Zhou et al. 2008, Chen 2009). These water masses exhibit considerable spatial and temporal variability, with distinct hydrographic properties and high nutrient content that support robust phytoplankton growth (Furuya et al. 2003, Chen 2009). Such environmental gradients, both horizontally across the shelf and vertically through the water column promote eutrophic conditions and high primary productivity, making the BYECS an ideal natural laboratory for investigating phytoplankton community dynamics and their response to environmental change.
Phytoplankton ecology in the BYECS has been widely studied, with research covering taxonomy, size structure, floristic distribution, environmental adaptation, and HABs (Zhou et al. 2008, Song 2010, Guo et al. 2011). Previous findings indicate that phytoplankton composition and spatial distribution are primarily governed by marine environmental conditions and oceanographic processes. These include the YS cold water mass, Changjiang (Yangtze) River discharge and the intrusion of the Kuroshio current into shelf regions (Liu et al. 2012, Wang et al. 2014, Song et al. 2017, Zhao et al. 2018). However, most earlier studies have concentrated on the Changjiang River estuary and its adjacent waters (Furuya et al. 2003, Zhou et al. 2008, Guo et al. 2014, Liu et al. 2016). Research has often focussed on seasonal patterns (Yoon et al. 2003, Park et al. 2008, Kim et al. 2020) or interannual variability during summer, when stratification is strongest (Xu et al. 2019). Despite, their efforts, long-term, decadal-scale analysis of phytoplankton community structure across entire BYECS remain limited. Traditional studies have emphasized larger phytoplankton groups such as diatoms and dinoflagellates due to their significant contribution to total chlorophyll-a (TChl-a). These studies consistently highlight diatoms as the dominant phytoplankton group in the BYECS (He et al. 2009, Zhao et al. 2010, Guo et al. 2014, Liu et al. 2016). In contrast, smaller groups, including cryptophytes, cyanobacteria, and prymnesiophytes have often been overlooked because of their relatively lower biomass contributions and the methodological challenges associated with preserving small cells using standard fixation techniques (Mackey et al. 1996).
Satellite remote sensing offers an effective tool for investigating large-scale phytoplankton biomass variations and associated biogeochemical processes (McClain 2009). Global and regional Chl-a distributions and primary production have been previously studied using satellite observations (O’Reilly et al. 2000, Chinta et al. 2024, 2025, Kalhoro et al. 2025a, 2025b, 2025c). However, research on phytoplankton species composition remain constrained due to limited spatial-temporal coverage. Therefore, there is a pressing need to investigate the long-term impacts of the environment on phytoplankton functional groups (PFGs) in the BYECS using remote sensing datasets. To address this gap, the present study examines the spatiotemporal variations in PFGs in the BYECS over the past two decades (2003–2022) using satellite remote sensing observations. This study also evaluates the impact of sea surface temperature (SST) and key nutrients (nitrate [NO3] and phosphate [PO4]) on the seasonal distribution of PFGs. Furthermore, it assesses the environmental drivers regulating phytoplankton community dynamics using a suite of statistical analyses, including the Theil-Sen slope estimator, the Mann-Kendal (MK) trend test, the false discovery rate (FDR) adjustment, correlation analyses, breakpoint detection, principal component analysis (PCA), heatmap correlations, and spatially weighted regression (SWR) analyses maps. Collectively, this study enhances understanding of the interactions between phytoplankton communities and environmental variability and contributes to broader assessments of the ecosystem health and productivity in the BYECS.
MATERIALS AND METHODS
Data collection
The spatiotemporal distribution, seasonal variability, and abundance of PFGs across the BYECS were analyzed using Copernicus Marine Service (CMEMS) reanalysis datasets (Fig. 1). The PFGs included diatoms (DIATO), dinoflagellates (DINO), prokaryotes (PROKAR), Prochlorococcus (PROCHLO), micro-phytoplankton (MICRO), nanoplankton (NANO), green algae (GREEN), Haptophytes (HAPTO), and pico-phytoplankton (PICO). Monthly PFG datasets for 2003–2022 were obtained CMEMS reanalysis, providing multi-sensor ocean-color biogeochemical products at 4 km resolution (https://data.marine.copernicus.eu/product/OCEANCOLOUR_GLO_BGC_L4_MY_009_104/description; product ID: DOI:https://doi.org/10.48670/moi-00281; https://data.marine.copernicus.eu/product/OCEANCOLOUR_GLO_BGC_L4_MY_009_108/download?dataset=c3s_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M_202207; product ID: DOI:https://doi.org/10.48670/moi-00283). Detailed product description, user manuals, and algorithm methodology for the PFGs are documented in Xi et al. (2020, 2021). These studies develop the PFGs retrieval algorithm and provide comprehensive assessments of phytoplankton functional type Chl-a (CHL), offering indirect but meaningful validation of the CMEMS PFGs products. However, satellite-derived datasets, some degree of algorithmic uncertainty and product-specific limitations may influence long-term trend interpretation. SST data with resolution of 0.05° × 0.05° and nutrient data (NO3 and PO4), were retrieved from CMEMS (https://data.marine.copernicus.eu/product/SST_GLO_SST_L4_REP_OBSERVATIONS_010_011/download?dataset=METOFFICE-GLO-SST-L9994-REP-OBS-SST_202003; product ID: DOI:https://doi.org/10.48670/moi-00168 [Worsfold et al. 2024]; https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_BGC_001_029/download?dataset=cmems_mod_glo_bgc_my_0.25deg_P1M-m_202406; product ID: DOI:https://doi.org/10.48670/moi-00019).
Spatial analysis
Monthly CMEMS satellite reanalysis datasets for PFGs, SST, NO3, and PO4 for 2003-2022 were aggregated into seasonal means (DJF, winter; MAM, spring; JJA, summer; SON) to examine large-scale spatiotemporal patterns across the study region. The TChl-a was calculated as the sum of CHL concentrations from all PFGs:
, where, PFGi represent the CHL contribution of each functional group and n is the total number of groups. The fractional contribution (%) of each PFG to the total biomass was computed and followed to Pan et al. (2011).
SWR was applied to quantify the spatial relationships between SST and PFGs, and between nutrients and PFGs, to evaluate environmental controls on phytoplankton variability. Analysis was performed using a geographical weighted regression (GWR) framework (Fotheringham et al. 2002) to account for spatial non-stationary in the relationships between PFGs and SST and NO3. A bi-square kernel with an adaptive bandwidth was applied, and optimal bandwidth (number of nearest neighbours) was selected by minimizing the corrected Akaike Information Criterion (AIC). All spatial coordinates were projected into a Universal Transverse Mercator (UTM) coordinate system, and Euclidean distance was used to compute spatial weights. Predictor variables were standardized using z-score normalized prior to model fitting. The resulting local GWR coefficients represent standardized effect sizes, with positive and negative values indicating the spatially strength and direction of SST and NO3 influence on PFG distributions (Brunsdon et al. 1996).
, where (ui, vi) are spatial coordinates, βk(ui, vi) are location wise specific coefficients, and is error term, and xik are predictor.
Statistical analysis
Annual time series from 2003 to 2022 were analysed to evaluate trends and relationships among PFGs, SST, NO3, PO4, and CHL. Complete observations were used for trend analysis, while pairwise deletion was applied for correlation analysis. All statistical tests were two-sided with a significant level of α = 0.05 and p-values were adjusted using the Benjamini-Hochberg (BH) false discovery rate (FDR) (Benjamini and Hochberg 1995).
MK test
Monotonic trends without distributional assumptions were evaluated using the MK test (Mann 1945, Kendall 1975). For a time-ordered series of size n, the MK statistic S and Kendall’s τ are defined as:
After correcting for ties, the variance of S under the null hypothesis of no trend is given by:
, where tp is the number of tied values for the p−th tie. The standardized test statistic Z is then computed as:
, with the two-sided p-value calculated as: p = 2Φ(−|Z|).
Theil-Sen slope
Trend magnitude was estimated using the Theil-Sen median slope (Sen 1968). For time points ti with i < j, the pairwise slopes are:
Generalized least-squares (GLS) with AR(1) errors
To account for serial autocorrelation, generalized least-square linear trends with AR(1) errors were fitted using the Prais-Winsten transformation (Prais and Winsten 1954). The model was:
The first observation and subsequent observations are transformed as:
The slope is tested with the usual t statistic based on the transformed-model covariance:
The parameter ρ is estimated interactively from ordinary least squares (OLS) residuals and regression is refit until convergence.
Breakpoint analysis
Structural changes in the time series were assessed using single-break (piecewise linear) models (Bai and Perron 2003). Candidate split years tk were restricted to ensure that at least a fraction γ = 0.20 of observations remained on each side. The breakpoint was identified as the year that minimized the total residual sum of squares (RSS).
Continuity at the breakpoint was not enforced. Pre- and post-break slopes were reported together with the estimated breakpoint year and corresponding value. For multiple testing, the BH step-up procedure was applied. Given m p-values sorted in ascending order as p(1) ≤ ⋯ ≤ p(m), the adjusted p-values were computed as:
Adjusted q values are then mapped back to the original hypothesis order.
Correlation and heatmaps
Cross-variable associations were quantified using Pearson’s correlation or Spearman’s rank correlation ρs (Pearson 1895, Spearman 1904). For large samples, the t − statistic:
, with n − 2 degrees of freedom (applied to r or ρs, respectively), and resulting p-values are adjusted by BH-FDR (Student 1908). Variables were hierarchically clustered using a distance 1 − |R| (average linkage), and the reordered correlation matrix R was visualized as a diverging heatmap. Where relevant, significance levels were indicated with symbols for q < 0.05, 0.01, and 0.001.
Standardize linear models
To assess the independent effects of SST, NO3, and PO4 on CHL and all PFGs, linear models were fitted using z-scored predictors:
All variables were standardized before model fitting, ensuring that the reported coefficients represent standardized effects coefficients (per 1 SD change in each predictor) (Cohen et al. 2003). Models were estimated using OLS on complete cases, and standardized coefficients β̂ were reported with 95% confidence intervals (CI). Two-sided p-values were adjusted for multiplicity across all responses-predictors combinations using the Benjamini–Hochberg FDR procedure (q = 0.05) (Benjamini and Hochberg 1995). For interpretability partial R2 values for each predictor were computed using a drop-one approach:
( R f u l l 2 - R r e d u c e d 2 ) / ( 1 - R r e d u c e d 2 )
Principal component analysis
Suitability for PCA was evaluated using the Kaiser-Meyer-Olkin (KMO) and Bartlett’s test of sphericity. Let R be the Pearson correlation matrix and Ω = R−1 the precision matrix. The absolute partial correlation between variables i and j is:
The overall KMO statistic is computed as:
, and the variable-wise measure of sampling adequacy is:
Following standard guidance, KMO > 0.6 (and most MSAi > 0.5) indicates sufficient shared variance for PCA (Kaiser 1974, Kaiser and Rice 1974). Bartlett’s test evaluates whether R departs from the identity. With n observations and p variables, the test statistic is:
, which is approximately chi-square (χ2) with p(p − 1)/2 degrees of freedom. A p-value < 0.05 supports the use of PCA (Bartlett 1950). PCA was then applied to the standardized data matrix Z by eigen-decomposing the correlation matrix:
, where the columns of V are component loadings and Λ = diag(λ1, ..., λp) contains the eigenvalues (component variances). Observation scores were computed as: S = ZV. Percent variance explained as 100λk/∑jλj and present a scree plot (Cattell 1966). The biplot displays (1) observation scores in the PC1–PC2 plane (points 1–20) and (2) variable loadings as plotted arrows, with lengths proportional to variable PC correlations. In correlation-based PCA, these are given by VΛ1/2 (Jolliffe 2002, Legendre and Legendre 2012, Jolliffe and Cadima 2016). All variables, including SST, NO3, and PO4, were included as active variables.
Overall, PCA, heatmaps, summarized Sen’s slope summaries, MK significance tests, GLS trends, and correlation analysis were collectively applied to identify dominant patterns in PFG distributions and their environmental drivers from 2003 to 2022. These methods provide integrated evaluation of long-term trends, seasonal variability, and ecosystem dynamics across the BYECS region. The methodological flowchart (Fig. 2) summarizes the complete analytical workflow for the spatiotemporal analysis of phytoplankton in the BYECS.
RESULTS
Decadal PFGs seasonal climatology
The spatial distribution and seasonal variability of major PFGs (DIATO, DINO, PROKAR, PROCHLO, MICRO, and PICO) across the BYECS (Fig. 3). Seasonal patterns were evaluated for DJF, MAM, JJA, and SON over the two-decade period (2003–2022). DIATO exhibited consistently high concentration (>2 mg m−3) along coastal zones in all the seasons, however offshore waters remained comparatively low (<0.8 mg m−3). Maximum DIATO abundance occurred during DJF and MAM (>2 mg m−3), particularly along nearshore areas and BS coastal regions. Slightly lower concentrations during JJA, likely influenced by localized upwelling and nutrient inputs (Fig. 3). DINO showed uniformly low concentrations (<0.4 mg m−3) across all regions and seasons, indicating minimal spatial or seasonal variability. PROKAR exhibits a strong seasonal signal, with high abundance (>4.8 mg m−3) during DJF and SON, particularly in coastal regions. These periods correspond to cooler temperatures and potentially higher nutrient availability, which may favor PROKAR growth. Higher concentrations were observed at BS and river regions. The abundance of PROKAR decreases significantly (<4.4 mg m−3) during MAM and JJA, reflecting sensitivity to warmer conditions and changing nutrients. PROCHLO, characteristic of oligotrophic waters, maintained a uniform distribution with consistently low abundance (<0.8 mg m−3) across seasons. The lack of significant seasonal variation suggests limited factors such as temperature, light availability and nutrient concentration limits the growth. MICRO abundance remains low across all seasons, with a slight increase in DJF and MAM (>1.2 mg m−3), primarily near the coast. This pattern indicates that larger phytoplankton (diatoms and dinoflagellates) may be limited by nutrient availability or hydrodynamic conditions. PICO shows a pronounced seasonal pattern, with peaks during SON and DJF (>4.8 mg m−3), particularly along the coastal areas (BS), YS margins, and ECS regions, and around river discharge area (Fig. 3). These distributions reflect their adaptability to a wide range of environmental conditions and capacity to thrive under nutrient-poor or stratified conditions. The high abundance during autumn suggests favorable conditions, such as post-monsoon nutrient inputs, which may enhance PICO growth. Overall, decadal-scale patterns show that DIATO, PROKAR and PICO are consistently exhibit higher in BS and YS river-discharge zones, while offshore waters remain relatively oligotrophic (Fig. 3).
Seasonal CHL and SST climatology
The seasonal climatology of CHL and SST in the BYECS (2003–2022) highlights distinct patterns (Fig. 4). CHL concentration peak in winter (DJF: ~1.5–2.5 mg m−3) and spring (MAM: ~1.2–2.0 mg m−3), driven by strong vertical mixing and spring bloom development. Summer (JJA: ~0.8–1.5 mg m−3) exhibits reduced CHL due to profound thermal stratification. While localized high (>5.2 mg m−3) occurred near the river discharge during MAM and JJA because of nutrient-enrich freshwater inputs. Autumn (SON: ~1.0–1.8 mg m−3) often shows a secondary increase associated with post-summer mixing. SST displayed lowest pattern in winter (DJF: ~10–16°C) and highest in summer (JJA: ~24–28°C), while spring (MAM: ~16–22°C) and autumn (SON: ~18–24°C) represent transitional phases (Fig. 4). Spatial analysis shows the inverse CHL-SST relationship observed at higher CHL during cooler, nutrient-rich conditions (winter/spring) particularly along the coast, and lower CHL under warmer, stratified summer conditions. The autumn CHL rise coincides with decreasing SST, reflecting nutrient reintroduction. Overall, these patterns underscore the strong influence of temperature-driven stratification and nutrient dynamics on phytoplankton biomass in the BYECS (Fig. 4).
Nutrients and surface pCO2 distribution
The seasonal distribution of nitrate (NO3), phosphate (PO4), and sea surface partial pressure of CO2 (pCO2) in the BYECS reveals significant variability over the past two decades (Fig. 5). NO3 concentrations peak in DJF and MAM along the coastal regions, particularly near the Yangtze River estuary (20–60 μmol L−1), driven by strong vertical mixing and riverine nutrient inputs. Offshore waters show lower concentrations (<10 μmol L−1). In JJA, coastal NO3 remains high (up to 60 μmol L−1) near the river plume, while offshore regions experience depletion (<0.5 μmol L−1) due to stratification and intense phytoplankton uptake. During SON, offshore concentrations partially recover (10–20 μmol L−1) as mixing increases, whereas coastal waters maintain elevated levels (20–40 μmol L−1) (Fig. 5). PO4 concentrations were generally lower than NO3 but followed similar seasonal patterns. Coastal PO4 ranged from 0.3–0.9 μmol L−1 and increase slightly in spring (0.5–1.0 μmol L−1), supporting seasonal phytoplankton blooms. Summer PO4 remained elevated near the river plume (0.5–0.9 μmol L−1), while offshore waters were depleted (<0.2 μmol L−1). Autumn values stabilized along the coast (0.3–0.8 μmol L−1) with minor offshore increase (up to 0.4 μmol L−1) (Fig. 5). pCO2 exhibits strong seasonal variability from the study area and was meansured in Pascal (Pa). Winter (DJF) values were lowest (10–20 Pa) due to enhanced CO2 solubility under cooler temperatures and strong vertical mixing. Spring values increased moderately (20–30 Pa), particularly offshore, influenced by remineralization of organic matter, riverine CO2 inputs, and warming temperatures, while coastal areas remained lower due to enhanced primary production consuming CO2. Summer exhibited the highest pCO2 offshore (40–50 Pa) as stratification limited gas exchange and reduced CO2 solubility, with coastal pCO2 moderately elevated (30–40 Pa). In autumn, pCO2 declined (20–40 Pa) as cooling and restabilised mixing enhanced CO2 uptake (Fig. 5). The regions with high PFG biomass correspond to lower pCO2, likely reflecting CO2 consumption through photosynthesis under favorable temperature and nutrient conditions. Overall, these patterns emphasize the critical role of riverine inputs, biological processes, and physical mixing in regulating nutrient and carbon dynamics across the BYECS.
SWR of PFGs with SST
Fig. 6 illustrates the SWR between SST and PFGs across the BYECS using monthly data from 2003 to 2022. The pronounced seasonal and regional heterogeneity in the SST-PFG relationship. In BS, DIATO and DINO showed predominantly weak to moderate negative SST associations during DJF and MAM (≈ −1 to −3), particularly in the northern basin. In contrast, strong positive relationships developed in the southern BS during JJA and SON (≈ +3 to +6), indicating an enhanced temperature response in warmer, shallow waters. In the YS, DIATO, and DINO exhibited some of the highest SWR values across the entire region during summer, with central and southern YS reaching +5 to +7. Meanwhile the cold-water mass in the northern YS retained negative coefficients (≈ −2 to −4). PROKAR and PROCHLO showed contrasting temperature responses. PROKAR exhibited strong positive SST relationships along river-influenced margins and southern coastal regions of the ECS (≈ +3 to +6.5), while northern offshore waters maintained negative values (≈ −1.5 to −3). PROCHLO exhibited strong negative SST response in the system during JJA, with extensive area of the YS and offshore ECS showing negative SWR values (≈ −4 to −7), reflecting sensitivity to warm stratified conditions. MICRO and PICO also demonstrated mixed response. MICRO showed positive SST associations along the Jiangsu-coastal Zhejiang coastal belt and the Yangtze River plume during warm months (≈ +2 to +5), but strongly negative patterns in the northern YS during all seasons (≈ −3 to −6). PICO displayed a distinct band of high positive SWR values in the southern YS and coastal ECS in JJA (≈ +4 to +6), while the BS and northern YS remained weakly negative (≈ −1 to −2). CHL exhibited moderate but spatially coherent SST sensitivity. Positive associations dominated in the Yangtze River estuary, nearshore ECS and northern BS (≈ +1 to 3.5), whereas offshore YS regions showed persistent negative patterns (≈ −1 to −3), especially during stratified seasons. Overall, the strongest positive SWR responses (≈ +5 to +7) were concentrated in warm, nutrient-enhanced southern coastal zones, while the strongest negative values (≈ −5 to −7) occurred in cold, stratified offshore regions, especially for PROCHLO and MICRO (Fig. 6). This SWR analysis underscores the complex temperature-driven dynamics of PFG distributions in the BYECS and highlights key seasonal and regional shifts in community structure resulting from SST variability.
SWR of PFGs with NO3
The SWR analysis between NO3 and PFGs revealed clear regional contrasts across the BYECS (Fig. 7). Overall, positive NO3-PFG relationships (≈ +2 to +7) were concentrated in nutrient-enriched coastal and river-influenced zones, while negative relationships (≈ −2 to −7) dominated offshore, stratified, and nutrient-poor regions. In the BS, DIATO, MICRO, and CHL exhibited strong positive SWR values in winter and spring (≈ +3 to +6), consistent with high winter mixing and riverine nutrient inputs. During summer, stratification led to weaker and neutral values (≈ 0 to +2). In the YS, the southern shelf showed moderate positive NO3 effects on DIATO, MICRO, and PICO during MAM and SON (≈ +2 to +4). By contrast, the northern YS displayed negative coefficients (≈ −2 to −5), reflecting nitrate depletion and reduced phytoplankton sensitivity nutrient variability. Across the ECS, the strongest positive NO3 dependence occurred along the Yangtze River plume during JJA (≈ +5 to +7), especially for MICRO, PICO, and PROCHLO, indicating nitrate-driven summer productivity. Conversely, offshore ECS waters consistently showed strong negative SWR values (≈ −4 to −7) for all PFGs, suggesting the low-nutrient, oligotrophic waters weaken the NO3-phytoplankton linkage (Fig. 7). Overall, the SWR patterns confirm that nitrate availability is a dominant driver of coastal phytoplankton variability, with its influence diminishing sharply toward offshore and stratified environments.
Fraction of total CHL for each group
The fractional contribution maps revealed distinct spatial and seasonal dominance patterns among PFGs across BYECS (Fig. 8). DIATO contributed modestly to TChl-a, generally 5–15% across most regions, with slightly higher fractions (15–20%) in the BS during winter and spring. DINO maintained consistently low contributions (<10%) in all seasons, reflecting their minor role in total biomass. PROKAR exhibited strong regional dominance, contributing 25–40% along the southern YS and ECS shelves during spring-autumn, and reaching >45% in localized coastal zones, especially in SON. PROCHLO maintained low fractional contributions (5–15%) throughout the year, with minor enhancements in offshore, warm oligotrophic waters during summer. MICRO displayed substantial contributions (20–35%) across the YS and ECS, peaking at 35–45% in summer, particularly in stratified near-shore waters. PICO was one of the most dominant groups, contributing 30–45% of TChl-a in warm, nutrient-poor offshore regions, especially during SON (Fig. 8). Overall, MICRO and PICO represented the major components of total phytoplankton biomass across the study area, while larger phytoplankton groups (DIATO, DINO) played a secondary but locally important role in colder, nutrient-enriched northern waters.
PFGs heatmap with environmental variables
The Pearson correlation heatmap revealed significant interactions among PFGs and environmental variables across the BYECS (Fig. 9). SST exhibits a strong negative correlation with PROKAR (r = −0.72, q < 0.002), HAPTO (r = −0.70, q < 0.001), PICO (r = −0.70, p < 0.003), and TChl-a (r = −0.65, p < 0.002, q < 0.001), indicating that elevated temperature conditions were associated with reduced abundance of these groups and overall chlorophyll levels. Conversely, SST was moderately positive correlated with DIATO (r = 0.47, p < 0.036, q > 0.369), and MICRO (r = 0.28, p > 0.24, q < 0.001), suggesting that diatoms and microplankton are more tolerant of, or responsive to warmer conditions. Among the PFGs, very strong and highly significant positive correlations were observed between PROKAR and PICO (r = 0.99, p < 0.001), PROKAR and HAPTO (r = 0.94, p < 0.001), and HAPTO and NANO (r = 1.00, p < 0.001), suggesting ecological coupling and similar environmental responses. DIATO exhibited a strong positive relationship with GREEN (r = 0.91, p < 0.001) and with CHL (r = 0.80, p < 0.001), highlighting its key role in driving species average total chlorophyll variability (TChl-a). DINO was moderately correlated with TChl-a (r = 0.61, p < 0.004), further emphasizing the influence of dinoflagellate biomass on total pigment concentration. Nutrient parameters showed a weaker association, such as NO3 displayed no significant correlation with most PFGs (p > 0.34), while PO4 exhibited weak negative correlations with DIATO (r = −0.67, p < 0.001) and MICRO (r = −0.64, p < 0.002), suggesting limited nutrient control compared to thermal effects (Fig. 9). Overall, these correlation patterns indicate that SST was the dominant environmental driver influencing long-term phytoplankton community composition and TChl-a across BYECS. The strong co-variation among several PFGs further reflects shared ecological niches and adaptive strategies under varying environmental conditions (Fig. 9).
Descriptive and long-term trends analysis of PFGs
The long-term time-series of PFGs and environmental variables revealed substantial temporal variability and distinct trends across the BYECS over the past two decades (Table 1, Fig. 10). The mean SST was 17.24 ± 0.39°C ranging from 16.39 to 17.98°C. A significant warming trend was observed (τ > 0.48, p < 0.002), with a Theil-Sen slope of 0.037°C y−1. A distinct breakpoint in 2009 indicated accelerated warming, increasing from 0.02°C y−1 before to 0.11°C y−1 afterward, suggesting intensifying oceanic heat accumulation of the last decade. Among the PFGs, PICO (1.70 ± 0.95 mg m−3) and PROKAR (1.59 ± 0.96 mg m−3) exhibited the highest mean abundances, confirming their dominance in the phytoplankton community. However, both groups showed strongly decreasing trend (τ < −0.83 to −0.87, p < 0.001) with breakpoints around 2015, reflecting a substantial decline in smaller phytoplankton types coinciding with rising SST. Similarly, HAPTO and NANO also showed significant negative trends (τ < −0.52, p < 0.001), with reduced abundance after 2011, suggesting a long-term shift towards lower nanoplankton biomass. In contrast, DIATO (0.41 ± 0.15 mg m−3) showed moderate mean concentrations and weak positive trends, indicating minor but consistent increase in larger phytoplankton groups. DINO averaged 0.20 ± 0.07 mg m−3, with a slight, non-significant decline (τ < −0.27, p > 0.083) (Fig. 10).
Among minor groups, GREEN showed a positive trend (τ > 0.45, p < 0.01), reflecting localized growth under changing nutrient-temperature conditions. The mean TChl-a concentration was 3.38 mg m−3 (±3.74 mg m−3), with values reaching up to 17.70 mg m−3, indicating episodic phytoplankton blooms. Despite these events, TChl-a showed a significant declining trend (τ < −0.57, p < 0.0003; Sen’s slope = −0.358 mg m−3 y−1), suggesting reduced total phytoplankton biomass over time. Conversely, CHL displayed a moderate positive trend (τ > 0.43, p < 0.01; Sen’s slope = 0.011 mg m−3 y−1), which may reflect compositional changes or differential pigment contributions among size fractions. Nutrient levels were relatively stable, with NO3 averaging 4.10 ± 0.28 μmol L−1 and PO4 0.03 μmol L−1 (p < 0.5), showing no significant long-term variations. Overall, the combined descriptive and trend analyses reveal a warming-driven reorganization of phytoplankton communities. The period from 2009 to 2015 marked a major ecological transition, characterized by the decline of small phytoplankton (PROKAR, PICO, HAPTO) and relative increase in larger groups (DIATO, MICRO, CHL). These changes point to a climate-induced shift in phytoplankton productivity and composition in the BYECS during the past two decades (Table 1, Fig. 10, Supplementary Table S1).
Community composition analysis indicated that PICO (28.99%) and PROKAR (27.86%) dominated, undertaking the importance of smaller, bacteria-like organisms. CHL contributed 17.63%, while MICRO (8.44%), DIATO (5.28%), and DINO (3.17%) played key roles in primary production (Fig. 11). The decline of PICO and PROKAR, coupled with increases in DIATO and MICRO, suggests a shifts in ecosystem structure likely influenced by elevated SST and nutrient dynamics.
Correlation analysis
The simple Spearman correlations (Fig. 12A) indicated that SST showed a strong and significant negative relationship with several PFGs, including PROKAR (r = −0.64, q > 0.018), HAPTO (r = −0.73, q < 0.005), NANO (r = −0.73, q > 0.017), and TChl-a (r = 0.64, q > 0.012). Similarly, PICO also showed a negative trend (r = −0.69) though not FDR-significant. These patterns suggest that warmer conditions are associated with reduced abundances of smaller phytoplankton such as PICO, HAPTO, and NANO. In contrast, DIATO (r = 0.38, q < 0.022) and CHL (r = 0.54) exhibited positive associations with SST, implying that diatoms and overall chlorophyll biomass may benefit from moderate temperature increases. Nutrient correlations were generally weak and nonsignificant. However, DIATO (r = −0.65, q > 0.826) and MICRO (r = −0.62) showed negative correlations with PO4, indicating that lower phosphate concentrations may limit larger phytoplankton (Supplementary Table S2).
The partial Spearman correlations (Fig. 12B), controlling for the other two environmental variables, reinforced SST as the primary regulator. Significant negative relationships persisted for PROKAR (r = −0.80, p < 0.003), DINO (r = −0.55, p < 0.0003), NANO (r = −0.84, p < 0.0006), and HAPTO (r = −0.84, p < 0.0004). In contrast, PROCHLO (r = 0.52, p > 0.052) and GREEN (r = 0.48, p > 0.043) showed weak positive associations with SST. Notably, after accounting for nutrient effects, nitrate and phosphate exhibited much weaker or inconsistent relationship, suggesting that temperature-driven changes dominate phytoplankton variability in the BYECS (Supplementary Table S3).
Multiple regression analysis
The multiple regression analysis (Fig. 13) further quantified these interactions. SST consistently emerged as the most influenced predictor across functional groups. Strong and statistically significant standardized coefficients were observed for HAPTO (β = −0.93, p < 0.001, q < 0.0007) NANO (β = −0.93, p < 0.001, q < 0.052), PROKAR (β = −0.87, p < 0.001, q > 0.089), PICO (β = −0.87, p < 0.0001, q < 0.048), and TChl-a (β = −0.86, p < 0.0006), confirming the strong inhibitory effect of warming on smaller phytoplankton and total biomass. Conversely, PROCHLO (β = 0.51, p > 0.062, q < 0.0007) and CHL (β = 0.42, p > 0.103, q < 0.0013) showed positive temperature responses, implying a shift toward certain cyanobacteria and mixed phytoplankton under higher temperatures (Fig. 13, Supplementary Table S4).
Among nutrients predictors, NO3 had positive but weaker effects on PICO (β = 0.35, p < 0.035, q < 0.004), HAPTO (β = 0.38, p < 0.017, q > 0.121), and NANO (β = 0.38, p < 0.017, q > 0.198) indicating that moderate nitrate enrichment may stimulate small-sized phytoplankton growth. PO4 exhibited significant negative effects on DIATO (β = −0.55, p < 0.016, q < 0.090) and MICRO (β = −0.66, p < 0.009, q < 0.05), suggesting nutrient limitation of larger phytoplankton under low phosphate conditions (Fig. 13, Supplementary Table S4).
Overall, both correlation and regression analysis consistently indicate that SST is the key environmental driver shaping phytoplankton composition and productivity in the BYECS. Rising temperatures during the past two decades likely caused a decline in smaller picophytoplankton (PROKAR, PICO, HAPTO, NANO) and relative increase in diatoms, reflecting temperature-induced restricting of the phytoplankton community.
Principal component analysis
The PCA was conducted to explore the relationships between environmental variables (SST, NO3, and PO4), and phytoplankton community composition. The Kaiser-Mayer-Olkin (KMO) measure confirmed the sampling adequacy (KMO > 0.639), and Bartlett’s test of sphericity was highly significant [X2 (91) = 974.49, p < 0.005], confirming strong correlations among the 14 biological and environmental variables and indicating the dataset was suitable for the PCA (Fig. 14). The first two principal components captured 79.5% of the total variances (PC1: 49.3%, PC2: 30.2%), while the remaining components each contributed less than 10%, as confirmed by sharp decline of the scree plot (Fig. 14A & B). PC1 was primary influenced by the SST (SST, loading = −0.295) and representing a temperature-driven gradient. Positive PC1 scores strongly linked to smaller phytoplankton groups, including PROKAR (0.375), PICO (0.374), HAPTO (0.350), NANO (0.350), and TChl-a (0.273), which are typically associated with cooler, more productive waters. In contrast, negative PC1 scores were associated with higher SST and favored PROCHLO (−0.191), DIATO (−0.247), and GREEN (0.244), indicating that warmer waters support a larger-celled taxas such as diatoms and green algae. These results suggest that SST is the principal environmental driver shaping PC1 and regulating the overall size and composition of the phytoplankton community.
PC2 represented a nutrient gradient primarily linked to phosphate availability. PO4 exhibited a strong negative loading (−0.274) and NO3 showed a weak contribution (loading = −0.016), indicating that PC2 was largely structured by phosphate rather than nitrate. Phytoplankton groups responding most strongly to PC2 included DINO (0.433) and MICRO (0.443), which displayed the highest positive loadings, demonstrating that dinoflagellates and microphytoplankton increased markedly under higher phosphate conditions. Conversely, PROCHLO showed a strong negative loading along PC2 (−0.364), indicating its dominance under low-phosphate, oligotrophic conditions. The biplot scores further illustrated this pattern: years with high PC2 scores (Fig. 14, Supplementary Tables S5–S7).
The PCA biplot revealed that the two principal axes, SST (PC1) and phosphate (PC2) act orthogonally to structure phytoplankton communities. High PC1 and low PC2 values dominated by the PICO, while low PC1 and high PC2 scores supported diatoms and dinoflagellates. These findings demonstrate the dual influence of SST and PO4 in shaping phytoplankton dynamics, whereas NO3 plays a comparatively minor role. Overall, the PCA clearly identifies SST as the primary environmental driver shaping the PFGs in the region. Collectively, the PCA provides a robust framework for understanding the environmental regulation of phytoplankton community structure and can inform predictive models of primary productivity in the study region.
DISCUSSION
Phytoplankton composition plays a pivotal role in food web dynamics, fisheries production, and the assessment of HABs (Field et al. 1998). The BYECS is a dynamic marine environment influenced by multiple hydrographic water masses, including Kuroshio water, Taiwan warm current coastal water, and Changjiang River diluted water (CDW), which regulate nutrient distribution and drive phytoplankton community structure (Chen and Wang 1999, Zhang et al. 2007, Guo et al. 2013). This two-decade study (2003–2022) provides a comprehensive evaluation of spatial, seasonal, and long-term ecological variability of PFGs, emphasizing their responses to environmental changes and contributions to marine productivity. Present findings show that smaller phytoplankton, particularly PICO (28.94%) and PROKAR (27.88%), dominated the oligotrophic regions of BYECS, reflecting their ecological importance in nutrient-poor waters. This aligns with previous studies reporting the prevalence of small phytoplankton under low-nutrient conditions (Furuya et al. 2003, Liu et al. 2016).
Seasonal pattern showed that these groups peak during winter months, with higher abundances in coastal and river-influenced zones, reflecting nutrient input from CDW and vertical mixing (Wang et al. 2013, Ndah et al. 2019). DIATO and MICRO, in contrast exhibited moderate increases after 2017, and contributed up to 40% of TChl-a under nutrient-enriched conditions, consistent with previous reports linking diatom blooms to riverine and upwelling nutrient inputs (Zhang et al. 2007, Huang et al. 2019, Pan et al. 2011). Previous research has demonstrated that long-term nutrient variations in BS influenced by riverine inflow, and diffuse sources, with over 40 rivers contributing to the system, among Yellow River is the second largest in China (Wang et al. 2019). These rivers are the dominant factor for nutrients in BS and YS (Xu et al. 2010, Wei et al. 2015, Yang et al. 2018). Present results indicate that nutrient rich coastal zones support phytoplankton growth. However, long-term trends reveal that a warming-driven reorganization of the phytoplankton community. SST increased significantly (τ > 0.48, p < 0.002; Theil-Sen slope = 0.037°C y−1), with a breakpoint around 2009 marking accelerated warming (from 0.02 to 0.11°C y−1). Concurrently, smaller phytoplankton such as PROKAR and PICO decline sharply after 2014 (τ < −0.83 to −0.87, p < 0.001), whereas DIATO and MICRO exhibited weakly positive trends. TChl-a decreased significantly over time (τ < −0.56, Sen’s slope = −0.36 mg m−3 y−1), reflecting overall biomass reduction, suggesting shifts in community composition. These results indicate climate-induced restructuring of the PFGs, consistent with global observations of warming-driven shifts from PICO- and NANO to larger diatom-dominated communities (Behrenfeld et al. 2006, Polovina et al. 2008, Liu et al. 2025).
Similarly, previous studies have identified diatoms as a dominant phytoplankton group in the ECS, South China Sea, and Vietnam (Xu et al. 2001, 2019, Liu and Chen 2012, Ndah et al. 2019, Sun et al. 2022). During inter-monsoon period, high production let DIATO dominated phytoplankton community near Dongha Island (Li et al. 2016). While, in the Jiulong River Estuary, DIATO dominate in spring, and DINO in summer, together making up over 50% of the biomass (Xu et al. 2001). DINO blooms also observed along aquaculture zones Pearl River Estuary in the south China Sea (Qi et al. 2004). The sustained increase in DIATO and MICRO further supports the notion that nutrient availability, particularly PO4 limitation, plays a fundamental role in shaping phytoplankton community composition (Li et al. 2007, Yoon et al. 2015). The period from 2009 to 2015 marked a major ecological transition, characterized by the decline of small phytoplankton and relative increase in larger groups point to a climate-induced shift in phytoplankton productivity and composition in the BYECS during the past two decades. This suggests that the BYECS is predominantly influenced by oligotrophic conditions, favoring smaller, nutrient-efficient phytoplankton taxa. The phytoplankton composition strongly associated with the seasonal variations (Ndah et al. 2019). In addition to nutrients, SST is a key environmental factor influencing changes in phytoplankton composition (Liu et al. 2025). Coastal regions generally exhibit higher phytoplankton abundance compared to offshore waters, likely due to nutrient inputs from terrestrial runoff and upwelling processes (Liu et al. 2022). Wang et al. (2013) also found the high primary production occur in winter compared to the summer due to the mixed layer depth stratification. Comparison with previous studies reveals both consistent and novel findings. The dominance of PICO and PROKAR in oligotrophic waters align with earlier research (Furuya et al. 2003, Liu et al. 2016). Similarly, the increase in DIATO and MICRO under nutrient-enriched conditions corroborates previous observations (Huang et al. 2019). The larger phytoplankton dominated the ECS and adjusted waters during spring (Furuya et al. 2003, Lou and Hu 2014, Liu et al. 2016). The river input likely drives the dominance of MICRO and NANO in coastal waters (Zhang et al. 2018), Jiang et al. (2015) also observed micro-sized DIATO and DINO along the riverside of ECS. They both represent major PFGs capable of causing HAB, which can pose risks to human health and negative impact on fisheries and aquaculture (Moore et al. 2008, Anderson et al. 2012). This study is limited by reliance on satellite and large-scale datasets, which may not capture the fine-scale variability, episodic events, or other ecological factors such as light availability and grazing pressure, including some uncertainty in predicting phytoplankton responses to environmental changes. Future studies should integrate in-situ observations and along with satellite-based data, considering regional, coastal, and offshore zones, and applying isobath-based masks, using advanced machine learning models, to enhanced accuracy and provide a mechanistic understanding of PFGs and their response to environmental change. However, recent studies on phytoplankton have demonstrated that satellite-derived ocean color data are also suitable for this type of analysis from study area (Wang et al. 2015, Sun et al. 2017).
Overall, the present analysis indicates that SST is the primary driver shaping PFGs in the BYECS. Heatmap, correlation, and regression analysis consistently show that warming strongly reduces smaller phytoplankton (PROKAR, PICO, HAPTO, NANO) and TChl-a, while favoring larger groups such as DIATO, MICRO, PROCHLO, and GREEN. Nutrients, particularly PO4 and NO3, exert weaker and less consistent effects (Yoon et al. 2015, Liu et al. 2025). Although shifts in nutrients compositions can significantly influence PFG dynamics (Guo et al. 2014); however, present study indicates that SST was the dominant driver of PFGs variability in the BYECS over the past two decades. PCA further confirmed SST as the main environmental gradient (PC1: 49.30%), with phosphate (PC2: 30.20%) structuring larger taxa, illustrating temperature-driven changes dominate long-term shifts in community compositions. KMO and Bartlett’s test verified the suitability of the dataset for PCA, and PCA loading were interpreted ecologically, linking environmental gradients to PFG distribution. These findings highlight a climate-induced transition from small, oligotrophic-adapted taxa toward larger, nutrient-responsive phytoplankton, providing critical insights into future primary productivity and ecosystem responses under ongoing warming. Furthermore, this research lays a foundation for sustainable marine resource management and conservation efforts in this ecologically and economically significant region.
This study provides a comprehensive analysis of the spatiotemporal variations in phytoplankton composition in the BYECS over the past two decades, utilizing satellite remote sensing data. The findings reveal significant seasonal and long-term trends in phytoplankton community structure, driven primarily by SST, with nutrient availability playing a secondary role. The dominance of smaller phytoplankton groups, such as PICO and PROKAR, underscores their ecological importance in oligotrophic waters, particularly during cooler months. However, their pronounced decline post-2015 highlights their sensitivity to warming SST, indicating temperature-driven reorganization of the phytoplankton community. In contrast, larger phytoplankton groups, such as DIATO and MICRO, exhibited increased abundance in nutrient-rich coastal regions, particularly during periods of upwelling and riverine nutrient inputs. SWR and regression analysis confirmed SST as the dominant environmental driver, with smaller phytoplankton negatively impacted by warming, while larger tax and certain cyanobacteria (PROCHLO, GREEN) responded positively. Nutrient availability, particularly PO4 for MICRO and DINO, and NO3 in coastal river-influenced zones, influenced phytoplankton composition locally, but their effect was weaker than thermal regulation. PICO and PROCHLO thriver in warm, oligotrophic offshore waters, reflecting adaptation to low-nutrient conditions.
Long-term trends revealed a warming SST (~0.037°C y−1) and a decline in the abundance of small phytoplankton (PROKAR, PICO, HAPTO, NANO), coupled with a relative increase in larger groups (DIATO, MICRO), indicating a climate-induced shift in community composition. These changes suggest a potential restructuring of the BYECS ecosystem, likely influenced by climate variability rather than significant nutrient decline. It is important to note that these interpretations are subject to uncertainties associated with satellite-derived PFGs estimates, which may be affected by algorithm limitations, spatial aggregation, and regional weighting. Additionally, local environmental process such as riverine inputs, vertical stratification, and water column mixing can modulate observed chlorophyll distribution, potentially influencing the apparent long-term trends in community composition. The findings align with and expand upon previous research, offering novel insights into the adaptive strategies of phytoplankton communities under warming and nutrient-limited conditions. This study highlights the ecological significance of phytoplankton diversity in sustaining marine productivity and underscores the sensitivity of these communities to environmental shifts. The observed trends have critical implications for understanding the impact of rising temperatures on marine ecosystems, particularly in the context of climate-driven changes in community composition and primary production. By integrating satellite remote sensing with long-term ecological data, this research provides a robust foundation for predicting future changes in phytoplankton dynamics and their cascading effects on marine food webs and biogeochemical cycles. These insights are essential for informing sustainable marine resource management and conservation strategies in the BYECS, a region of immense ecological and economic importance.
PDF Links
PubReader
ePub Link
Full text via DOI
Download Citation 
Supplement
Print 
