Samples (n = 124) were collected from a near-enclosed mariculture bay at the water layer of the surface (0.5 m) and bottom (2–55.5 m) from the Sansha Bay, Fujian Province during four cruises (January, April, July, and October 2019, i.e., winter, spring, summer, and autumn) (Figure 1). Water samples were prefiltered through a 200 μm pore mesh. A portion of water samples were directly filtered onto a 0.22 μm filter membrane (47 mm diameter; Millipore, United States) to obtain bacterial samples. Another portion of water samples were prefiltered through a 20 μm pore mesh and then sequentially filtered to 2 μm and 0.8 μm filter membranes (47 mm diameter; Millipore, United States) to obtain the samples of pico-protists (0.8–2 μm) and nano-protists (2–20 μm). These filters were stored at −80°C until DNA extraction.

Water temperature, salinity and depth were measured in situ with CTD (AML Base X, Canada). Dissolved oxygen (DO) were measured in situ with WTW (Multi3630 IDS, German). Nitrate (NO₃-N), nitrite (NO₂-N), ammonia (NH₄-N), phosphate (PO₄-P), silicic (Dsi), chlorophyll a (Chla), as well as dissolved total nitrogen (TN), total phosphorus (TP), were measured following the standardized method (Xie et al., 2020). Dissolved inorganic nitrogen (DIN) content was calculated by the sum of NO₃-N, NO₂-N, and NH₄-N.

Water samples pre-filtered at 20 μm were fixed with glutaraldehyde at a final concentration of 0.1% (V/V) and the abundance of heterotrophic bacteria (HB), synechococcus (Syn) and photosynthetic picoeukaryotes (PPE) were run on a FACSAria flow cytometer (Becton Dickinson, United States) equipped with laser emitting at 488 nm (Marie et al., 1997, 2000). To measure the abundance of HNF (heterotrophic nanoflagellates) and PNF (pigmented nanoflagellates), water samples pre-filtered at 20 μm were fixed with glutaraldehyde at a final concentration of 0.5% (V/V), then stained with 4′6-diamidino-2-phenylindole (DAPI), and filtered at low pressure (<100 mm Hg) onto 0.8 μm pore-size black nuclear pore filters (25 mm diameter; Millipore, United States). Filters were mounted to glass slides and stored at −20°C in the dark until observed by epifluorescence microscope (Leica Microsystems, Germany). Following Guo et al. (2020) methods, at least 50 randomly selected fields of view were examined for each sample.

Conversion of enumerated microbial abundance to carbon-containing biomass with reference to survey data in the adjacent sea area. The abundance of HB and Syn was converted to the biomass using a factor of 20 fg C cell⁻¹ and 178 fg C cell⁻¹, respectively (Lee and Fuhrman, 1987; Charpy and Blanchot, 1998). The abundance of PPE was converted to biomass using a factor of 1,500 fg C cell⁻¹ (Zubkov et al., 1998). The abundance of HNF and PNF was converted to biomass using a factor of 4,700 fg C cell⁻¹ (Chen et al., 2009).

For protists, DNA was extracted using DNeasy PowerWater kit (Qiagen, United States), the V4 region of the 18S rRNA gene was amplified with primers TAReuk454FWD1 (5′-CCAGCA(G/C)C(C/T)GCGGTAATTCC-3′) and TAReukREV3 (5′-ACTTTCGTTCTTGAT(C/T)(A/G)A-3′) (Stoeck et al., 2010). A quintuple repetition of each sample was amplified as follows: 95°C for 5 min, 10 cycles of 94°C for 30 s, 47°C for 45 s, and 72°C for 1 min, 25 cycles of 94°C for 30 s, 47°C for 45 s, and 72°C for 2 min, and a final extension at 72°C for 2 min. Forward and reverse primers were tagged with 2 bp links and 8 bp barcodes to allow the pooling of multiple samples in one run of sequencing and later differentiation of different samples (Kozich et al., 2013). The PCR products were purified using Agarose gel DNA Recovery Kit (Bioteke, China) and quantified by Nano-200 (Allsheng, China). Samples belonging to two size-fractionated were, respectively, mixed in equimolar concentrations to construct two individual amplicon libraries for sequencing using Illumina Hiseq 2500 platforms with PE250 strategy (Illumina, United States) according to standard protocols. The raw sequence data have been deposited in the NCBI Sequence Read Archive under the BioProject accession number PRJNA1252866.

For bacteria, 16S rRNA gene sequences were obtained from the BioProject number PRJNA747131 and the Accession number SRP328863 (Zhu J. et al., 2022). Briefly, DNA was extracted using the FastDNA spin kit (MP Biomedicals, USA), the V3-V4 region of the 16S rRNA gene was amplified with primers 341F (5′-CCTAYGGGRBGCASCAG-3′) and 806R (5′-GGACTACNVGGGTWTCTAAT-3′) (Liu et al., 2019). The triplicate PCR products were purified using GeneJET gel purification kit (Thermo Scientific, United States) and quantified by Qubit 2.0 fluorescence quantifier (Thermo Fisher Scientific, Waltham, United States). PCR products were mixed in equal amounts and sequenced on THE Illumina HiSeq platform (Illumina, United States).

Sequencing data processing was performed on Mothur v.1.47.0 following MiSeq standard operating procedure¹ with the steps of sequencing data quality control, amplicon sequence variants (ASV) clustering, and species classification. Specifically, read pairs were aligned, and the tags and primers of reads were removed. To reduce sequencing and PCR error, reads with only one or two sequences were removed. Chimeras were detected using the UCHIME algorithm (Edgar et al., 2011), and if there were flagged chimeras, they were removed from all samples. The remaining high-quality reads with suitable lengths were clustered to distinguish ASVs. The SILVA nr v.138 database (Quast et al., 2012) and PR² protist v4.14 database (Guillou et al., 2012) were used for bacterial and protist taxonomic assignment, respectively. To avoid distortion of the relative abundance of DNA sequences of microbes, non-bacterial or non-protistan sequences (e.g., “unknown,” Archaea, Nucleomorphs, Fungi, Streptophyta, and Metazoa) were removed. To minimize the occurrence of ASVs at both size-fractionated due to cell dislodging or filter clogging, ASVs were further filtered with reference to Deutschmann et al. (2023). For each ASV that appeared at both size-fractionated, the ratio of its relative abundance at the nano-and pico-protist communities was calculated, and when the ratio exceeded 2, the abundance of pico-protists was removed. When the ratio was lower than 0.5, the abundance of nano-protists was removed. Lastly, 124 valid samples were obtained, including 40 winter samples, 30 each of spring and summer samples, and 24 autumn samples. A total of 9,025 bacteria ASVs, 12,274 pico-protist ASVs, and 12,908 nano-protist ASVs were retained.

Referring to previous literature reports (Adl et al., 2012; de Vargas et al., 2015b), the functional groups of pico-and nano-protist ASVs were annotated into photoautotroph, symbiont, parasites, heterotroph, and mixotroph. Here, mixotrophs refer to photosynthetic species that also ingest food by phagocytosis or osmotrophy, while symbiont refers to heterotrophic species that retain prey plastids or symbionts (Adl et al., 2019) (Supplementary Appendix A).

All statistical analyses were performed in R v.4.2.1 software (R Core Team, 2019), and visualized using the “ggplot2” R package (Wickham, 2016). The α-diversity indices, including Shannon, Chao1, and Pielou’s evenness, were calculated using the “vegan” R package (Oksanen et al., 2022). Differences in α-diversity among seasons were analyzed with multiple comparisons (LSD test), using the “agricolae” R package (Mendiburu, 2010).

Prior to multivariate statistical analyses, the ASV tables were Hellinger-transformed and the environmental variables were standardized to zero mean and unit variance by “vegan” R package (Oksanen et al., 2022). The β-diversity was calculated on the Bray-Curtis distance metric and visualized with nonmetric multidimensional scaling (NMDS). The significant differences in microbial communities among water layers were tested by permutation multivariate analysis (PERMANOVA) of variance using the Bray-Curtis distance metric. To explore the correlations between environmental factors and their effects on microbial communities, Spearman correlation analysis and partial Mantel test were performed using the “microeco” R package (Liu et al., 2020). For the partial Mantel test, the distance matrix of microbial communities was calculated using the Bray-Curtis method. Additionally, in order to reveal the biotic interactions between size groups, the affecting biotic factors for each size group did not contain the abundance data of the same size of the organisms. To further identify characteristic taxa, the linear discriminant analysis (LDA) effect size (LEfSe) was used to screen for classes that differed significantly between seasons (LDA > 4.0) and analyze.

To assess the environmental adaptability of microbial communities, the community-level niche breadth was calculated, which represents the average ecological tolerance of all taxa in a community. A higher niche breadth indicates a predominance of generalist species, while a lower value reflects a dominance of specialists with narrower ecological preferences (Pandit et al., 2009). This index was calculated as the mean of Levins’ niche breadth values for all taxa within a community (Pandit et al., 2009; Wu et al., 2018). Levins’ niche breadth index (B) was calculated following the description of Pandit et al. (2009) using the “spaa” R package (Zhang, 2016): B j = 1 / ∑ i = 1 N P ij 2 , where B_j represents the habitat niche breadth of ASV j in a metacommunity; N is the total number of communities in the metacommunity; and P_ij is the proportion of ASV j in resource state i, i.e., the abundance of ASV j in community i divided by the abundance of ASV j in the metacommunity. The B value ranges from [1, N] and the higher value indicates that the ASV j is widely and evenly distributed in the metacommunity. Community-level niche breadth was calculated for all seasons, and metacommunities were defined as the set of communities in each season. Differences in community-level niche breadth among communities were analyzed with the LSD test.

All samples were used for network inference (n = 124). To focus on widespread interactions among community members and their relative influence on network properties, ASVs shared by all seasons, presented in more than 1/3 of samples and with an average relative abundance greater than 0.001 were retained (n_Bac = 42, n_Pico-protist = 85, n_Nano-protist = 100). Networks were generated using the “SpiecEasi” R package relying on sparse neighborhood and inverse covariance selection algorithms (Kurtz et al., 2015). Interdomain and intradomain networks were constructed independently using the sparse and low rank (SLR) method by setting parameters (nlambda = 20, minimum lambda ratio = 0.005, pulsar threshold = 0.05, number of representatives = 20) (Shekarriz et al., 2023).

Network topology parameters (i.e., degree, absolute, and betweenness centrality) were calculated using the “igraph” R package, and network modularization and visualization were performed using Gephi software. Referring to the code provided by Shekarriz et al. (2023), we analyzed the efficiency and random attack robustness of networks, which were proposed by Latora and Marchiori (2001) and Iyer et al. (2013), respectively. Normalized robustness (R) of a network was calculated by running iteratively for 10,000 times, using the formula R = 1 N ∑ i = 1 N σ ( 1 N ) , where N is the initial size of the network, σ is the relative size of the largest network component after node removal, and i is the number of vertex or vertices removed from the network. The vulnerability (V) of a network is related to the robustness calculated as V = 0.5 − R . Multiple comparisons of the network properties were conducted using the LSD test. Keystone species were identified as nodes with degree and betweenness centrality measures with values in the top 20 percentiles (Roume et al., 2015).

To further reveal the effects of environmental variables on microbial communities and microbial food web structure, a partial least squares path model (PLS-PM) was modeled using the “plspm” R package (Sanchez et al., 2024). The model included microbial community diversity, community-level niche breadth, and network information transfer efficiency. Network information transfer efficiency index reflects how effectively signals or interactions propagate through the microbial network, and can serve as a proxy for the functional coordination and resilience of the microbial system. Following García-Comas et al. (2016) and Yang et al. (2018), we also included the predator/prey biomass ratio (PPBR) as a proxy for the actual trophic transfer efficiency between predator and prey to assess its ecological relationship to mathematically inferred network connectivity and flow. To improve the reliability of the model, the PPBR was log10 transformed, run using 1,000 bootstraps, and variables with loadings < 0.7 were removed. Based on the path coefficients, the direct and indirect effects of other potential variables on network robustness and information transfer efficiency were calculated. The performance of the model was evaluated using the Goodness-of-Fit (gof) measure. The final PLS-PM model included seven variables: water property (temperature, salinity, and DO), water nutrition (TN and TP), microbial community diversity (only pico-protist and nano-protist), community-level niche breadth (all three groups), PPBR (the biomass ratio of HNF/bacteria and HNF/PPE), network robustness, and network information transfer efficiency.

Environmental conditions varied more in temporal than in spatial variations, with little difference in environmental parameters between two depths in the same season (Supplementary Figure S1). The environment was nitrogen nutrient-rich, with average PO₄-P, DIN, and nitrogen-to-phosphorus ratios (TN: TP) reaching eutrophic levels in autumn and winter, and middling nitrogen-rich nutrient levels in spring and summer (Supplementary Table S1). Temperatures, DO, PO₄-P, and NH₄-N in Sansha Bay were distinct in all seasons (p < 0.05, Supplementary Figure S1). Microbial biomass also fluctuated seasonally (p < 0.05). In winter, the abundance of pico-photosynthetic microorganisms (PPE and Syn) was at a low level, while the abundance of PNF reached its peak throughout the year. The abundance of HB and HNF showed a similar trend of being higher in winter and summer, and lower in spring and autumn (Supplementary Figure S1).

Departing from previous investigations by Ma et al. (2021) and Zhu J. et al. (2022) that focused on bacteria and protists (0.22–200 μm) diversity, here we delved into the diversity of pico and nano size-fractionated protists. We provided a detailed description of changes in community composition and interdomain networks, aiming to elucidate the connections between protist and bacterial communities. The diversity indices of the three groups were significantly different (p < 0.05, Figure 2). The Shannon index and Peilou’s evenness index of the nano-protist community were the highest, followed by the pico-protists and bacteria, while the chao1 index was the opposite. The α-diversity indices of bacterial and protist communities showed greater seasonal than depth variation (Figure 2). The pico-protist and nano-protist community diversity indices were comparable, both fluctuated upward over time and appeared to have complementary seasonal fluctuations to the bacteria (Supplementary Figure S2).

The community structure of bacteria and protists also showed remarkable seasonal variations (Figure 3). The results of the Mantel test showed that temperature, DO and salinity had significant effects on bacterial and protist community structure. Among them, the bacterial and nano-protist community structure has the highest correlation with DO (Mantel’s r_bac = 0.65 and r_nano-protist = 0.74), while the pico-protist community structure has the highest correlation with temperature (Mantel’s r = 0.79). Biological factors also influenced microbial community structure, including PPE and PNF for bacterial communities (Mantel’s r_PPE = 0.26 and r_PNF = 0.19), PNF for pico-protist communities (Mantel’s r = 0.41), and Syn and PPE for nano-protist communities (Mantel’s r_Syn = 0.43 and r_PPE = 0.14). PERMANOVA further revealed significant differences in the communities of all size-fractionated groups in the four seasons (R² > 0.17 and p < 0.001), but there were no significant differences between depths and habitats (R² < 0.02 and p > 0.05, Supplementary Table S2). Based on the above results, subsequent analysis focused on the seasonal dynamics of microbial communities.

The community composition showed pronounced seasonality (Figure 4). The bacterial community was mainly composed of the Alphaproteobacteria class of Proteobacteria phylum (65% ± 10%), followed by Actinobacteria phylum (25% ± 9%). The relative abundance of Alphaproteobacteria peaked in autumn. The Actinobacteria class and the Acidimicrobiia class in the Actinobacteria phylum showed opposite seasonal dynamics, with the relative abundance of the former reaching a minimum in the autumn (2% ± 1%) and the latter reaching a maximum in the autumn (8% ± 2%). Both had opposite responses to seasonal changes in environmental factors such as temperature, salinity, NO₂-N, NH₄-N, PPE biomass, and PNF biomass (Supplementary Figure S3). In addition, the average relative abundance of Cyanobacteria was highest in the summer (18% ± 13%) and below 5% in the rest of the season, showing a positive correlation with temperature, salinity, PNF biomass, and HNF biomass, and a negative correlation with nutrient concentrations (Supplementary Figure S3).

Two size-fractionated protist community composition were similar at supergroup level that both were dominated by Alveolata (26–84%), followed by Hacrobia (4–40%) (Figure 4). In the pico-protist community, two supergroups showed complementary seasonal dynamics, where the relative abundance of Stramenopiles peaked in the winter month (17% ± 8%), while Archaeplastida was higher in the other seasons, with the highest proportion in summer month (24% ± 11%). At the class level, the nano-protist communities contained a higher relative abundance of Dinophyceae (28% ± 11%) and a lower abundance of Syndiniales (27% ± 9%) than the pico-protist communities (8% ± 7 and 45% ± 12% respectively), and the other dominant classes differed slightly between the two size-fractionated communities. The response of dominant taxa to environmental factors was also essentially the same for both pico-protists and nano-protists of the same taxon (Supplementary Figures S4, S5). The winter dominant taxa were negatively correlated with temperature, salinity, NO₂-N, and Syn biomass, and positively correlated with DO, NH₄-N, and PO₄-P. The dominant taxa in other seasons were mainly positively correlated with temperature and negatively correlated with PO₄-P.

The protist communities were trophically divided into photoautotroph, heterotroph, mixotroph, symbiont, parasites, and unknown (Figure 5). From ASV counts, heterotrophs had the highest percentage (31%), parasites and photoautotroph following (24 and 12% respectively), and symbionts were the least abundant (3%). From the relative abundance, parasites were the most abundant in the pico-protist community (46% ± 12%), while mixotrophs were the most abundant in nano-protist communities (31% ± 8%), and symbionts accounted for the lowest percentage in both communities (3% ± 4 and 3% ± 3% respectively). Although the relative abundance of each trophic group differed in the pico-and nano-protist communities, the seasonal trends were generally similar (Supplementary Figures S6, S7). Heterotrophs were most abundant in winter (9–69%), showing negative correlations with temperature, salinity, Chla, and Syn biomass, and positive correlations with PO₄-P, NO₃-N, and NH₄-N. Parasites were abundant in autumn and winter (16–79%), and the correlations with other environmental factors were almost exactly opposite to those of the heterotrophs, except for no significant correlation with PO₄-P. Photoautotrophs, on the other hand, were relatively abundant in the warmer summer and autumn (6–45%) and were negatively correlated with PO₄-P. Seasonal dynamics of the mixotrophs differed between the two size-fractionated communities. Pico-mixotrophs were lowest in winter (4% ± 3%) and showed positive correlations with temperature and salinity, and negative correlations with phosphate and DO. Nano-mixotrophs, in contrast, were most abundant in winter (37% ± 7%), and consisted mainly of the red tide species Heterocapsa rotundata (25% ± 5%), which showed the opposite environmental correlation to the pico-mixotroph.

To evaluate the structural properties of seasonal microbial food webs, bacteria-pico-protist-nano-protist interdomain networks (BPN networks) were constructed using the SpiecEasi method (Figure 6; Table 1). Comparing the networks in the four seasons, the autumn network had the lowest modularity coefficient, average path length, the highest average degree, and clustering coefficient, and was also the most robust (V_mean = 0.085) and nodal efficient in information transfer (μ = 0.41). The summer network was completely opposite, vulnerable (V_mean = 0.100), and inefficient in the transfer of information (μ = 0.37).

There were no significant differences in the specific contributions of bacteria, pico-protists, and nano-protists to the structural properties of the BPN networks, such as absolute edge weight, degree, and betweenness (Supplementary Figure S8). The network had the highest percentage of nano-protists (34–52%), followed by pico-protists (46–28%), and the lowest by bacteria (20–32%). Network keystone species analysis showed that the nano-protist Dinophyceae class were keystone species in all four seasons, and the pico-protist MAST-4 class in winter, spring, and summer (Supplementary Table S3). These results indicated that protists were important foundational components of the interdomain network structure.

The role of bacteria in the network was non-negligible. In BPN networks, although bacteria constituted a small percentage, they were mostly keystone species (33–55% of the keystone species in the networks). In particular, the Rhodobacteraceae family were keystone species in all four seasons, and the Cyanobacteria in spring, summer, and autumn (Supplementary Table S3). It was also worth noting that the network analysis without distinguishing seasons revealed that the inclusion of bacteria improved the robustness and information transfer efficiency of the networks (Supplementary Figure S10). The bacteria-nano-protist interdomain network (BN network) was the most rubust and efficient among the interdomain networks, followed by the BPN network and the bacteria-pico-protist interdomain network (BP network). Therein, the average edge weights and degrees of bacteria were significantly higher than those of protists. In summary, these results highlighted crucial and highly connected keystone roles of bacteria in enhancing the robustness and efficiency of the Sansa Bay protists networks.

To investigate the reasons for the stability and efficiency of the microbial network in autumn, we calculated community-level niche breadth, which can reflect the community’s ability to adapt to the environment, and also used trophic transfer efficiency as represented by the PPBR. The microbial community niche breadth of all three taxa in autumn was the lowest of the year (Figure 7A), indicating that the microbial communities had a narrow environmental tolerance in autumn, which promoted the functional specialization of the community. Furthermore, complex trophic relationships existed between bacteria and nano-protists, and between pico-protists and nano-protists. The results of PPBR showed that the ratio of biomass of HNF to bacteria was significantly lower than that of HNF to PPE, with the ratio of biomass of HNF to PPE in autumn being significantly higher than that in other seasons (Figure 7B). These suggested that trophic transfer efficiency between pico-and nano-sized organisms was particularly high in autumn.

Based on the previous results, key water property factors (temperature, DO, and salinity), water nutrition (TN and TP), and related biological factors (diversity, community-level niche breadth, PPBR) that shaped the microbial community structure were selected to construct a PLS-PM model. This model was used to investigate the direct and indirect drivers influencing microbial food web stability and the efficiency of information transfer within Sansha Bay (gof = 0.722, Figure 8A; Supplementary Table S4). Water property was found to significantly influence not only microbial community diversity (path coefficient = 0.50, p < 0.01) and community-level niche breadth (path coefficient = −0.45, p < 0.01), but also microbial network robustness (path coefficient = 0.83, p < 0.01) and information transfer efficiency (path coefficient = −0.76, p < 0.01). The influence of water nutrition on microbial networks was weaker than that of water property and was completely opposite. Considering the total effects, the drivers demonstrated opposing impacts on microbial network robustness and information transfer efficiency (Figure 8B). Among all factors, community-level niche had the strongest influence on microbial network structure, followed by PPBR, and then water nutrition.

The α-diversity of pico-protist and nano-protist communities showed similar seasonal trends and structure (Figure 2; Supplementary Figure S9), whereas there were differences in composition (Figure 4), in particular the proportion of parasites and mixotrophs in the communities. Similar to the results of the Tara Oceans expedition for the oceans (de Vargas et al., 2015a), the pico-protist communities had a higher percentage of parasites compared to the nano-protist communities (Figure 5). On the one hand, this may be inflated due to the higher rDNA copy number in some marine alveolate lineages, specifically the clonal abundance of Syndiniales, the dominant class of the pico-community, overestimates the cellular abundance (Not et al., 2009; Massana, 2011). On the other hand, the parasites released hundreds of small, non-phagocytic dinospores into the water column after killing their host (Guillou et al., 2008; Siano et al., 2011). However, because the parasite hosts encompass a wide range of organisms from flagellates to fish, even studies suggesting that Syndiniales may not always be unequivocally host-specific, it was difficult to find a clear association between the parasites and the potential host and to elucidate the drivers of seasonal dynamics of parasites (Guillou et al., 2008; Käse et al., 2021). In addition, the trophic modes of protists may be flexible and influenced by environmental factors (Flynn et al., 2019). We acknowledged that categorizing taxa into discrete functional groups may oversimplify the ecological complexity and plasticity inherent in protists.

Mixotrophs are important components of microbial food webs (Stoecker et al., 2017). Previous studies confirmed that mixotrophic protists could obtain nutrients through direct ingestion of bacteria or algae (Flynn et al., 2019), and coastal PNFs were important grazers of Synechococcus populations (Tsai et al., 2007; Chan et al., 2009). This trophic relationship likely accounts for the significant correlations observed between nano-protist communities and Synechococcus biomass, and between bacterial communities and PNF biomass (Figure 3D). Some mixotrophic protists may also pose risks to aquaculture systems (Flynn et al., 2018). For example, the mixotrophic flagellate Heterocapsa pygmaea (Millette et al., 2017), identified as a keystone species in the spring, summer, and autumn interdomain networks (Supplementary Table S3), was the dominant species in a harmful algal bloom near Sansha Bay (Wu et al., 2022). Another mixotroph, Karlodinium veneficum, a keystone species in the autumn interdomain network (Supplementary Table S3), is known to cause fish kills during its blooms (Place et al., 2012). However, no algal blooms or fish kills occurred during our sampling, despite the relative abundance of Heterocapsa genus exceeded 20% in winter and the relative abundance Karlodinium veneficum exceeded 10% at several mariculture stations. The macroalga Gracilaria lemaneiformis widely cultured in Sansha Bay, has been shown to effectively inhibit harmful algal bloom formation, likely through trophic competition, allelopathic interactions (Yang et al., 2015)., and shading effects that suppress microalgal growth (Xie et al., 2021).

It is generally assumed that network efficiency varied inversely with robustness because edges that contain more connections between nodes make the network more resilient to attacks, but less efficient at creating the shortest paths between two nodes (Rodrigue, 2020). However, in our study, the autumn microbial network showed the highest robustness and the highest efficiency (Figure 6). The PLS-PM results revealed that the overwhelming effects of community-level niche breadth on the network’s robustness and efficiency than environmental factors (Figure 8). Given that protists had the narrowest ecological niche breadth in autumn (Figure 7A), we speculated that the remarkable specialization of protists in autumn was an important driver of rubust and highly efficient interdomain microbial networks.

Multiple environmental ecology studies have demonstrated that microbial specialists play a crucial role in biological interactions (Finke and Snyder, 2008; Yan et al., 2022; Li et al., 2024; Zhou et al., 2024). Synthetic community experiments further revealed the mechanism of microbial specialists’ roles. Specialists established a clearer functional division of labor by utilizing specific resources, such as specific carbon sources, nitrogen sources, or microhabitats (Huang et al., 2019). Specialization reduced the competitive pressure caused by overlapping resources, which made the species interaction more reciprocal or neutral, thus enhancing the stability of the network structure (Bai et al., 2015). When environmental stress increased, specialists played a more crucial role in stabilizing ecological networks (Li et al., 2023). Moreover, specialization also reduced redundant connections and improved the efficiency of information transfer (Tsoi et al., 2019; Wang et al., 2024).

Differently from protists, bacteria contributed to BPN networks not only by specialization. On the one hand, season-independent network construction revealed that the inclusion of bacteria contributed to a robust microbial network, evidenced that bacteria worked to connect members of the community and thereby exchange metabolites or were consumed. For example, Rhodobacteraceae, identified as key species in all four seasonal BPN networks (Supplementary Table S3), were known to have multifaceted and mutually infochemical exchanges with phytoplankton, such as diatoms and dinoflagellate, and their interactions dynamically change according to the physiological state of the phytoplankton (Amin et al., 2015; Seymour et al., 2017). On the other hand, the community diversity of more generalist bacteria was excluded from the PLS-PM analysis due to low loadings. This implied that a select few bacteria, which possess specific relationships with phytoplankton and protist communities, emerged as the keystone species of the interdomain microbial network. However, most generalist bacteria contributed minimally to the interdomain network structure, as their increased diversity did not translate into more efficient resource utilization compared to that of specialists (Finke and Snyder, 2008).

Despite the presence of macroalgae farms that absorbed a portion of the excess nutrients from cage cultures, the waters of Sansha Bay remained persistently eutrophic (Xie et al., 2020). PLS-PM results revealed that the total effect of water nutrients (TN and TP) on microbial network structure exceeded water quality property (Figure 8). This anthropogenic eutrophication resulted in significant changes in microbial community structure (Kiersztyn et al., 2019; Chen et al., 2024). PLS-PM results also showed significant effects of HNF predation on network robustness and information transfer efficiency (Figure 8), suggesting that predatory relationships were central to microbial food webs. HNFs as major predators of pico-sized microorganisms were controlled by food supply (Lin et al., 2016), playing a top-down control role in the microbial food web. Such biological interactions were critical for maintaining the stability of bacterial communities (Liu et al., 2022).

The community structure of marine microorganisms was determined by various abiotic and biotic environmental factors such as temperature, salinity, and nutrients (Logares et al., 2009; Sunagawa et al., 2015; Dai et al., 2022). In cultured waters, dissolved organic matter (DOM) released from macroalgae and fish culture was also a key factor influencing the microbial community structure (Yang et al., 2015; Luo et al., 2024; Ou et al., 2024). In Sansha Bay, year-round cultured Gracilaria lemaneiformis provided a steady input of polysaccharides while inhibiting phytoplankton blooms (Yang et al., 2015; Ou et al., 2024). After Saccharina japonica was harvested in May, there was less complex DOM remaining in autumn, and the resources available to microbes tended to be simpler, reducing the microbial community ecological niche (Tanentzap et al., 2019; Huang et al., 2024). Larimichthys crocea also experienced a rapid growth period in autumn (midway through the 8–13 month rearing cycle), and it had a continuous and steady input of metabolites, creating a predictable pattern of nutrient supply (Luo et al., 2024). Moreover, compared to summer, there were barely any typhoons in autumn, the hydrological conditions were stable, and physical perturbations had little effect on the microbial community structure (Nguyen et al., 2021).