From 2020 to 2022, ten seasonal research cruises were conducted in Laizhou Bay. Sampling was performed in May, August, and October 2020; January, March, August, October, and December 2021; and August and October 2022. Each year, more than 40 sampling stations were surveyed, with their locations shown in Figure 1B. In addition, a fixed coastal station in the eastern part of Laizhou Bay (Figure 1B) was monitored monthly from July 2021 to June 2022, yielding 12 high-frequency observations. Although this nearshore station does not fully capture the spatial heterogeneity of the entire bay, its time-series data provide valuable temporal resolution that complements the broader seasonal spatial surveys, thereby facilitating a more detailed understanding of the seasonal dynamics and succession of Synechococcus populations.

Surface seawater (0.5 m depth) temperature, salinity and dissolved oxygen (DO) were measured in situ using a YSI 556 multiparameter probe (Yellow Springs Instruments, United States), and pH was determined onboard using a glass electrode. For nutrient analysis, 300 mL of seawater was filtered through GF/F filters (0.7 μm pore size, Whatman) and stored at −20°C until analysis. Concentrations of nitrate, nitrite, ammonium, dissolved inorganic phosphorus (DIP), and dissolved silica (DSi) were determined using a continuous flow analyzer (Seal QuAAtro, Germany) following Jiang et al. (2022). Dissolved inorganic nitrogen (DIN) was calculated as the sum of nitrate, nitrite, and ammonium. For chemical oxygen demand (COD) analysis, 250 mL of seawater was collected and immediately frozen at −20°C. COD was determined using the standard dichromate oxidation method with reflux digestion and titration with potassium dichromate (Wei et al., 2022). For chlorophyll a (Chl-a) analysis, 1 L of seawater was filtered through GF/F filters and stored in liquid nitrogen. Chl-a concentrations were subsequently quantified using high-performance liquid chromatography.

Seawater samples were pre-filtered through a 20 μm mesh and fixed with 2% paraformaldehyde in 2 mL cryovials, kept in the dark for 15 min, and then rapidly frozen in liquid nitrogen. PE-rich and PC-rich type Synechococcus populations were identified using a Becton Dickinson FACSCalibur flow cytometer equipped with dual lasers (488 nm and 635 nm), following the procedure of Liu et al. (2014). Four fluorescence channels were analyzed: FL1 (530/30 BP), FL2 (585/42 BP), FL3 (670 LP), and FL4 (661/16 BP). PE-rich cells were identified by strong orange fluorescence (FL2), while PC-rich type cells were distinguished by their red fluorescence (FL4).

Environmental DNA was extracted using an enzyme/phenol-chloroform protocol (Riemann et al., 2008), and eluted in TE buffer (10 mM Tris, 1 mM EDTA, pH 8.0). The rpoC1 gene of Synechococcus was amplified via nested Polymerase Chain Reaction (PCR) following the protocol of Mühling et al. (2006). The first round employed primers rpoC1-N5 and rpoC1-C, and the second round used barcoded primers rpoC1-39F and rpoC1-462R, which target a ~ 400 bp fragment of the rpoC1 gene.

To ensure spatial representativeness while maintaining feasibility, 9 to 13 samples were selected from each cruise for rpoC1 analysis. Selection was based on spatial distribution and prevailing weather conditions during each cruise. Successful rpoC1 amplification was achieved for nine of the ten seasonal cruises; no usable sequences were obtained from the December 2021 cruise due to sample degradation.

PCR products were gel-purified (Qiagen, Germany), quantified using the Quant-iT picoGreen double stranded DNA (dsDNA) assay kit (Invitrogen, USA), pooled in equimolar concentrations, and sequenced on a Roche GS Junior 454 pyrosequencing platform (Xia et al., 2015). Raw reads were quality-filtered using mothur (v1.35.1; Schloss et al., 2009), retaining sequences between 300 and 500 bp in length with average quality scores > 20. Denoising was performed with shhh.seqs using a sigma value of 0.01, and chimeras were removed via chimera.uchime (Schloss et al., 2011). High-quality sequences were then aligned and compared to a curated reference database of Synechococcus rpoC1 sequences. Each sequence was assigned to a known phylogenetic clade based on BLAST results (identity > 90%, E-value < 0.01). Sequences with < 90% identity were designated as unclassified (Xia et al., 2015). This genotype-level classification enabled detailed analysis of Synechococcus community structure and seasonal succession.

In this study, an exact amplicon sequence variant (ASV) approach was applied to the rpoC1 sequence data using the DADA2 pipeline, enabling single-nucleotide resolution without reliance on arbitrary clustering thresholds (Callahan et al., 2017). For phylogenetic analysis, representative ASV sequences were aligned using a 90% similarity cutoff, a common practice for functional genes with high sequence divergence such as rpoC1 (Turk-Kubo et al., 2015). This approach preserves the high-resolution taxonomic information inherent to ASVs while allowing robust tree construction. All downstream diversity and community composition analyses were conducted at the ASV level to ensure analytical precision and comparability with current microbial ecological standards.

From May 2020 to October 2022, ten surveys revealed pronounced temporal variability in the environmental conditions of Laizhou Bay (Figure 2). Surface water temperature exhibited a clear seasonal cycle, peaking in August (~ 28°C) and declining to around 17°C in May and October. Winter temperatures dropped to near 0°C (Figure 2A). Salinity was notably lower between October 2021 and October 2022, averaging around 25, with similarly low values (~ 27) also observed in October 2020.

DIN and DSi concentrations increased markedly in October 2021, more than doubling compared to earlier cruises (Figure 2B). In contrast, DIP remained consistently low (< 0.5 μmol L⁻¹) throughout the study period, showing irregular fluctuations without a clear seasonal trend. The molar ratios of DIN/DIP and DSi/DIP were substantially elevated, far exceeding the Redfield ratio (N/p = 16), with values reaching ~ 3,000 and ~ 2000, respectively, in October 2021 (Figure 2C). The DSi/DIN ratio peaked in August 2021 (~ 3), suggesting a relative accumulation of silicate. DO concentrations were generally lower in summer, particularly in August, but never dropped to hypoxic levels (Figure 2D). COD showed a slight increase from December 2021 onwards, while pH exhibited notable fluctuations, ranging from a low of 7.8 in August 2020 to a high of 8.7 in October 2021.

The cell abundance of Synechococcus exhibited pronounced seasonal variation, with significantly higher value in summer and autumn than in winter and spring (Figure 3). Peak abundances reached up to ~ 10⁶ cells mL⁻¹ in August, whereas values typically remained below 100 cells mL⁻¹ during colder seasons. Both PC-rich type and PE-rich type Synechococcus displayed similar seasonal trends; however, PE-rich type cells consistently dominated, with abundances approximately one order of magnitude higher than PC-rich type cells.

Horizontal distributions during summer and autumn further illustrated these seasonal dynamics (Figure 4). In summer, Synechococcus abundance decreased northward from the southern bay, where concentrations exceeded 10⁶ cells mL⁻¹, to values below 10⁵ cells mL⁻¹. Notably, a localized low-abundance patch (<10³ cells mL⁻¹) was observed in central to northern Laizhou Bay in August 2021. In autumn, high-abundance zones (>10⁵ cells mL⁻¹) emerged in both the northern and southern regions, while most other areas showed moderate levels (~10⁴ cells mL⁻¹). In contrast, spatial variability was minimal in winter and spring, with nearly all stations recording abundances below 10³ cells mL⁻¹ (Supplementary Figure S1).

PE-rich type Synechococcus consistently dominated the community composition across seasons. In summer 2020, PE-rich type cells accounted for over 70% of total Synechococcus abundance throughout most of the bay (Figure 5). In summer 2021, their contribution remained high (> 60%) across much of the area, although it dropped below 40% in parts of the southwestern and central-northern bay. In summer 2022, PE-rich type Synechococcus comprised approximately half of the total abundance, increasing from ~ 20% in the southwest to ~ 80% in the northeast. During autumn, PE-rich type cells remained the dominant pigment type, contributing 88, 94, and 73% to the total abundance in 2020, 2021, and 2022, respectively. They also prevailed in winter and spring, except in May 2020, when their contribution dropped to 29% (Supplementary Figure S2).

A total of 14 Synechococcus clades were identified across the ten cruises and fixed-point time series, including subclusters within clade S5.1 (S5.1. I – S5.1. IX, S5.1. XV), as well as clades S5.2 and S5.3. Among them, clade S5.1. I is typically associated with cold or temperate nutrient-rich coastal waters, while S5.1. II and S5.1. III are thermophilic ecotypes that typically dominate in warm surface layers. Clade S5.1. V has been observed in warm-temperate and estuarine environments. Ecological characteristics of all identified clades are provided in Supplementary Table S1.

In spring 2020, clade S5.1. I overwhelmingly dominated the assemblage, accounting for over 50% of sequences at all stations, and exceeding 80% in central bay areas (Figure 6A). During periods of Synechococcus bloom, genotypic diversity increased substantially. As shown in Figures 6B,C, in August and October 2020, multiple clades such as S5.1. II, III, V, VI, VII, VIII, and IX co-occurred, with the first four clades contributing substantially to the assemblage (> 50%). In summer and autumn of 2021, clades S5.1. VII and S5.1. III increased notably, each contributing approximately 30% to the total composition, while overall genotypic diversity remained high (Figures 6F,G). In August 2022, clade S5.1. I once again became dominant (> 70%) in northern Laizhou Bay, whereas other regions were primarily occupied by S5.1. II and S5.1. III (Figure 6H). By October 2022, the assemblage across the entire bay was again overwhelmingly dominated by S5.1. I (> 80%), accompanied by the reappearance of clade S5.1. IV (Figure 6I). In colder months (January and March 2021), the community structure shifted, with clades S5.1. I, V, VI, and S5.3 representing the main genotypic components (Figures 6D,E).

Monthly observations at the fixed station from July 2021 to June 2022 revealed clear seasonal patterns in environmental conditions, Synechococcus abundance, pigment types, and genotypic composition (Figure 7).

Environmental variables showed strong seasonal variation (Figures 7A,B). Sea surface temperature was above 24°C in July–August, dropped steadily to ~2°C by February, and then increased to 15°C in May and 21°C in June. Salinity also peaked in July–August (> 29), decreased from October onward, and remained below 25 from November through February before recovering to ~27 in late spring. DIN concentrations rose sharply in October (>60 μmol L⁻¹) and stayed elevated (>35 μmol L⁻¹) until February. DIP remained low throughout the year (<0.2 μmol L⁻¹), except for a moderate increase in February. DSi was relatively high (>20 μmol L⁻¹) in August and October, fluctuated around 10 μmol L⁻¹ for most of the year, and dropped below 5 μmol L⁻¹ in May and June.

A pronounced bloom of Synechococcus began in July, with cell densities exceeding 10⁵ cells mL⁻¹ from August to October (Figure 7C). Following the bloom, abundance declined sharply in November and remained below 10³ cells mL⁻¹ through the following May, before increasing again in June (>10⁴ cells m⁻¹). Throughout the year, PE-rich type Synechococcus consistently dominated the community, accounting for over 70% of total abundance, except in November when its relative contribution slightly declined.

Genotypic succession accompanied seasonal shifts in abundance (Figure 7D). In the early Synechococcus bloom stage (July), clades S5.1. V and VI were prominent, each comprising over 20% of the community. However, their proportions declined during peak bloom months (August and October), when S5.1. II and S5.1. III increased markedly, reaching 30 and 20%, respectively. During November and December, S5.1. VIII became dominant (>50%), along with a rise in S5.1. VI (~20%). From January to June, clade S5.1. I gradually increased and became the predominant genotype by June (>60%).

Redundancy analysis (RDA) based on data from ten cruises revealed the environmental drivers shaping Synechococcus abundance and community composition (Figure 8). Total Synechococcus abundance, as well as PC- and PE-rich type populations, showed strong positive correlations with SST, and negative correlations with DO and pH. However, the relative proportion of PE-rich type cells exhibited a negative association with both SST and COD, suggesting a potential shift in pigment-type dominance under warmer, more eutrophic conditions. Major genotypes such as S5.1. II, VII, VIII, and IX were positively correlated with SST and negatively correlated with DO, indicating temperature-driven niche preferences. In contrast, correlations between Synechococcus dynamics and nutrient concentrations or salinity were generally weak.

These findings were further supported by the monthly time series (Figure 7), where peak Synechococcus blooms consistently coincided with elevated temperature. Notably, although temperature was similar in October and May, Synechococcus abundance was markedly higher in October, which is further discussed in Section 4.2.

Our results demonstrate that temperature plays a pivotal role in shaping both the abundance and community composition of Synechococcus in Laizhou Bay. Throughout the seasonal surveys, Synechococcus consistently reached bloom-level densities from summer to mid-autumn (August and October), which coincided with the high seawater temperatures recorded (Figure 2A). Redundancy analysis further confirmed significant positive correlations between temperature and total Synechococcus abundance, as well as with several dominant pigment types and genotypic clades (Figure 8). These patterns are consistent with previous findings that temperature is a primary environmental factor governing the growth and distribution of picocyanobacteria in marine environments (Zwirglmaier et al., 2008; Mackey et al., 2013; Varkey et al., 2016).

To further elucidate the mechanistic basis of this temperature control, we examined the direct relationship between seawater temperature and Synechococcus abundance. The analysis revealed a strong exponential increase in cell density with rising temperature (Figure 9), consistent with the well-established exponential dependence of biogeochemical and metabolic rates on temperature (Eppley, 1971; Sanz-Lázaro et al., 2011). Notably, our data revealed that Synechococcus abundance began to increase more rapidly around 26°C, suggesting this value may represent an approximate thermal inflection point under the specific environmental conditions of Laizhou Bay (Figure 9). This suggests that 26°C may represent a critical bloom-initiation temperature under the specific environmental conditions of Laizhou Bay. We emphasize that this threshold was derived from statistical fitting and should not be interpreted as a strict Synechococcus bloom-initiation criterion. Indeed, time-series observations at the fixed station (e.g., October 2021) recorded high Synechococcus abundance under slightly cooler temperatures (Figure 7), likely influenced by lag effects or concurrent nutrient availability.

Threshold-type responses of cyanobacteria to environmental factors have been observed in various contexts, including light intensity, depth, and Chl-a concentrations (Zinser et al., 2007; Cram et al., 2015; Yeh and Fuhrman, 2022). Specifically, temperature thresholds triggering Synechococcus blooms have been reported in other marine systems, although the precise inflection points can vary depending on local ecological settings ranging from below 15°C to ~ 25°C in different studies (e.g., Johnson et al., 2006; Larkin et al., 2020). While the identification and mechanistic understanding of such threshold behaviors remain areas for further systematic investigation, our findings provide ecologically relevant evidence of a thermal inflection point under eutrophic semi-enclosed bay conditions. This contributes to the broader understanding of Synechococcus bloom dynamics and may inform predictive models in temperate coastal ecosystems facing future climate warming.

Temperature appeared to structure community composition by selecting for genotypes with distinct thermal niches. For example, several PE-rich type lineages (e.g., clades S5.1. II, III, and VII) were positively associated with temperature (Figure 8B), implying thermal differentiation among subpopulations. These observations are in line with the concept of “thermal ecotypes” previously reported in picocyanobacterial assemblages (Johnson et al., 2006; Kling et al., 2023), whereby even closely related genotypes exhibit divergent thermal performance curves. Previous studies have shown that clade S5.1. I typically inhabits cold or temperate waters and can occur in both coastal and oceanic environments (Farrant et al., 2016; Sohm et al., 2016), which supports our observation of its dominance during the spring (May 2020) and autumn (October 2022) cruises. Furthermore, elevated summer temperatures may enhance water column stratification, leading to reduced mixing and prolonged light exposure in surface layers, which benefits phototrophic organisms like Synechococcus, particularly in turbid coastal settings (Behrenfeld et al., 2006; Fischer et al., 2014).

Although the total abundance of Synechococcus increased with temperature, the relative proportion of PE-rich type cells exhibited a negative correlation with sea surface temperature (Figure 8A). This apparent decoupling between absolute growth and pigment-type composition suggests more complex regulatory dynamics. One possible explanation is that different pigment types possess distinct thermal optima or light-harvesting efficiencies, leading to shifts in competitive interactions under warming conditions (Six et al., 2007; Pittera et al., 2018). Alternatively, temperature may interact with other factors such as nutrient stress, pH, or grazing pressure, indirectly modulating pigment-type structure (Zwirglmaier et al., 2007; Grébert et al., 2018). These findings underscore the importance of jointly considering both abundance and community composition when evaluating picocyanobacterial responses to environmental drivers.

Synechococcus exhibited apparent spatial distributions in both cell abundance and community composition during the bloom periods (Figures 4–6), even though seawater temperatures showed no obvious difference across these times. Notably, Synechococcus biomass in October was markedly higher than in May, despite comparable thermal conditions (Figure 7). These divergences highlight that temperature alone cannot fully account for bloom dynamics and underscore the importance of other regulatory factors, such as nutrient availability and ecological interactions, in shaping Synechococcus population responses.

Our results reveal pronounced spatiotemporal variability in Synechococcus community composition in Laizhou Bay, encompassing both pigment types and genotypes. PE-rich types consistently dominated throughout most of the year, accounting for over 70% of total abundance during bloom periods (Figure 5), with their dominance only weakening in early spring (e.g., ~29% in May 2020). This pattern aligns with observations from other coastal systems, where PE-rich type Synechococcus are favored under high-light, stratified conditions due to their efficient light-harvesting antennae and broad ecological plasticity (Liu et al., 2014; Wang et al., 2022b; Li et al., 2024). In contrast, PC-rich type Synechococcus are generally confined to estuarine and low-salinity environments rich in red light, and are rarely found in offshore waters (Wood et al., 1998). However, in Laizhou Bay, frequent Yellow River-derived freshwater and sediment inputs (Liu, 2015) likely create turbid, low-salinity coastal conditions that support PC-rich type populations. In fact, PC-rich type Synechococcus occasionally exceeded PE-rich types in southern coastal areas in summer (Figure 5), highlighting the complex niche partitioning shaped by freshwater influence and local optical properties.

Genotypic data further emphasize the seasonal succession within the Synechococcus assemblage. The S5.1 clade, widely distributed in coastal waters and comprising over 20 known sub-lineages (Li et al., 2024), was the most prevalent in our study. Clade S5.1. I, typically associated with cold or temperate nutrient-rich waters (Zwirglmaier et al., 2008; Sohm et al., 2016), dominated during early spring and winter (Figures 6, 7D). Interestingly, it also re-emerged as the dominant lineage in May 2020 and October 2022 despite relatively warm temperatures, suggesting that temperature alone does not fully control seasonal succession. Other ecological drivers, such as viral pressure, selective grazing, or legacy effects from earlier nutrient pulses, may also influence genotype dynamics. Similar patterns of S5.1. I dominance under warm conditions have been reported in other coastal systems, such as the western North Atlantic (Hunter-Cevera et al., 2016), Southern California Bight (Tai and Palenik, 2009) and the Bohai Sea (Li et al., 2024).

During summer and autumn blooms, community diversity increased markedly and was often dominated by clades S5.1. II, III, and VII (Figure 6). These lineages are known to prefer warm, nutrient-rich surface waters and may have competitive advantages such as rapid growth rates and efficient nutrient uptake mechanisms under stratified conditions (Post et al., 2011; Sohm et al., 2016). Similar seasonal shifts in genotypic composition have been reported in other temperate and subtropical coastal systems. For example, Larkin et al. (2020) and Yeh and Fuhrman (2022) documented positive correlations between clade II and temperature in the western North Atlantic and the San Pedro Channel, respectively. In contrast, Tai and Palenik (2009) observed a summer dominance of clades I and IV, with clades II and III remaining at low abundance, in the Southern California Bight. At our fixed station, a distinct genotypic succession was observed from clades II, III, and VI in July–August to clade VIII in August–October of 2021. These patterns likely reflect intra-clade competition, niche differentiation, or varying responses to shifting light and nutrient regimes during Synechococcus bloom development.

Compared to other eutrophic coastal systems such as Jiaozhou Bay or the Pearl River Estuary, where Synechococcus genotypic shifts are often driven by riverine inputs and salinity gradients (Xia et al., 2015; Li et al., 2024), Laizhou Bay appears to be more strongly influenced by internal processes such as benthic nutrient regeneration, thermal stratification, and biotic interactions. Specifically, factors like viral lysis, allelopathic inhibition, or grazing could differentially affect the fitness of specific genotypes, driving community restructuring during different seasons (Suikkanen et al., 2004; Wang and Chen, 2004; Labrie et al., 2013). Future studies incorporating metatranscriptomic profiling or physiological trait analyses under controlled conditions will be crucial for deciphering the mechanisms underlying clade-specific responses and ecological succession patterns.

The spatiotemporal dynamics of Synechococcus in coastal systems like Laizhou Bay are likely to become increasingly complex under the dual pressures of long-term climate warming and intensifying short-term weather extremes. Rising sea surface temperatures may progressively shift the competitive balance among picocyanobacterial clades, favoring thermophilic lineages (e.g., clade S5.1. II and III) with higher thermal optima and metabolic plasticity (Flombaum et al., 2013; Kling et al., 2023). Moreover, the increasing frequency of extreme precipitation events, driven by anthropogenic climate change (Shiu et al., 2012; Pendergrass, 2018), may trigger abrupt and episodic pulses of freshwater and nutrient influxes, thereby reshaping local salinity, stratification, and nutrient availability. Such stochastic perturbations could disrupt established seasonal succession patterns, favoring opportunistic or fast-growing Synechococcus ecotypes, or even facilitating transitions to alternate bloom regimes.

To anticipate and quantify these future shifts, an integrated approach combining in situ long-term monitoring, culture-based physiological assays, and trait-informed ecosystem modeling is crucial. Recent advances in isolating and culturing diverse Synechococcus clades under controlled conditions have enabled the measurement of key physiological traits such as growth rate, nutrient uptake kinetics, temperature and light optima, and allelopathic potential (Labrie et al., 2013; Pittera et al., 2014; Hunter-Cevera et al., 2016; Browning et al., 2017; Fernández-González et al., 2020). These traits, when accurately characterized across representative pigment types and genotypes, can be incorporated into size- or trait-based ecological models (e.g., NEMURO, CoSiNE model) to simulate clade-specific responses to variable environmental drivers (Xiu and Chai, 2014; Gao et al., 2022). In highly dynamic coastal seas like Laizhou Bay, where seasonal forcing interacts with episodic hydrological inputs, such models offer a powerful tool for disentangling the nonlinear interactions shaping picocyanobacterial community assembly and for guiding coastal management under changing climatic conditions.