1 Introduction

Protists have long been recognized as either photoautotroph (typically diatoms and coccolithophores) or heterotrophic predators, which are major contributors to the ocean productivity (Falkowski et al., 2007) and carbon consump tion (Sherr and Sherr, 2002). In addition to the photosynthetic and heterotrophic growth, protists exhibit a range of other trophic modes, including mixotrophy, parasitism, sym biosis, osmotrophy, and saprotrophy (Worden et al., 2015). These complex behavioral strategies and organismal interactions have led to challenges for the comprehensive model of the marine carbon cycle (Worden et al., 2015). Protists are morphologically and genetically diverse and widespread in all types of marine habitats (Massana, 2015). In temperate seas, protist community undergoes strong seasonal changes in composition and abundance, as a result of alterations in both abiotic factors, e.g., irradiance, temperature and nutrient levels, and biotic factors, such as grazing, competition and interaction with pathogens and symbionts (Gran-Stadniczeñko et al., 2018).

Traditional studies of diversity, distribution and dynamics of marine protists were mainly based on microscopical identification and enumeration (Song et al., 2009, 2019; Liu et al., 2017, 2019; Pan et al., 2020). Due to the poor efficiency of the traditional methods, the protist diversity might have been underestimated (Keeling and Del Campo, 2017). The PCR-based approaches on the diversity of 18S rRNA gene have been widely applied in the last decade and revealed massive new lineages of aquatic protists (Moon-van der Staay et al., 2001; Lefranc et al., 2005). Recently, the high throughput sequencing techniques (HTS) has been applied to analyze the hyper-variable regions of the 18S rRNA gene, which has highlighted remarkable seasonal and spatial variations in marine protists (Gong et al., 2015; Sun et al., 2017; Li et al., 2019) and indicated that a large portion of taxa were still uncovered (De Vargas et al., 2015; Gran-Stadniczeñ ko et al., 2018). Compared to the microbiota in pelagic ocean, it is less known about the information of protist community in coastal areas (e.g., estuary and intertidal zone), the involved biotic interactions, and how the seasonal and chemicalphysical factors determine protist community structure and diversity. Coastal areas are potential hotspots of studies on microbial communities interacting with biogeochemical processes and adapting to environmental changes, such as salinity shifts and impacts of human activities (Balzano et al., 2015).

During the seasonal variation, autotrophic phytoplanktons can be affected by fluctuations of light, temperature, and/or nutrient concentration in the coastal area (Montagnes and Franklin, 2001; De la Vega et al., 2011; Vannier et al., 2016). As predators, heterotrophic protozoas are also indirectly infected by the change of chemicalphysical factors (Paffenhö fer, 1998). In addition, the organismal interactions (e.g., parasitism, symbiosis, and mixotrophy) have recently been recognized as major biotic factors contributing to the protist community structure (Caron et al., 2017; Stoecker and Lavrentyev, 2018). To our knowledge, neither abiotic nor biotic factors affecting protist communities have been addressed in the coast of the northern Yellow Sea of China, although implications from the protist community dynamics may provide a glimpse to the inducements of many endemic ecological events, e.g., red tide, jellyfish blooming, hypoxia, which are increasingly occurring in this coastal area (Wei et al., 2008; Wang, 2009; Xu et al., 2013).

In this study, using high-throughput sequencing (HTS) for the V4 region of 18S rDNA gene, year-round and monthly-taken samples were collected in two coastal sites around Yantai, a city along the northern Yellow Sea of China. Seasonal changes of planktonic protist community and the chemical-physical factors were investigated. Here we addressed the following questions: ⅰ) What is the seasonal distribution pattern of the protist community in the studied area? ⅱ) How do the biotic interactions happen among protist groups? ⅲ) What are the main drivers attributing to the seasonal pattern? By answering the above questions, we aim to better understand the dynamic and ecological function of planktonic protist communities in coastal ecosystems.

2 Materials and Methods 2.1 Sample Collection and Processing

Water samples were collected monthly from the surface layer (about 20 cm deep) during the high tide at two sites (X and Y, Fig. 1) in the coast of Yantai from March 2018 to February 2019. As shown in Fig. 1, station X is located at the estuary of Xinan River (37˚25΄N, 121˚37΄E), while station Y is located at an intertidal zone near the Yantai University, about 8 km away from station X (37˚28΄N, 121˚27΄E).

The water temperature, salinity, DO and pH were measured in situ using an electronic sensor (YSI, USA). Water samples for nutrient analysis were collected with Niskin bottles (15 mL) and triplicates were set for each sample. The ammonium (NH₄⁺), nitrate (NO₃⁻), nitrite (NO₂⁻), phosphate (PO₄³⁻) and silicate (SiO₄²⁻) were analyzed using a nutrient Autoanalyzer (Seal, Germany). The total organic carbon (TOC) and total organic nitrogen (TON) of water samples were determined with a Vario Micro Cube Elemental Analyzer (Elementar, Germany).

Protist communities were collected by filtrations of the in situ marine water (500 mL). To remove large plankton, the water samples were prefiltered through a 200 μm mesh. Protist cells were then collected with a peristaltic pump (JINTENG, China) through polycarbonate filters with a diameter of 45 mm a pore size of 0.22 μm (Millipore, Ireland). The filter membranes with protistan cells were immediately put into cryopreservation tubes and stored at −80℃ before DNA extraction.

2.2 DNA Extraction and High-Throughput Sequencing

Extraction and purification of DNA from the filter membranes were carried out using the Fast DNA SPIN Kit for Soil (Q-BIO gene, USA). The DNA concentrations were determined using a Nanodrop 2000c spectrophotometer (Thermo-Fisher, USA).

The V4 region of 18S rDNA has been successfully used to deduce the seasonal dynamics of the Skagerrak plankton communities (Gran-Stadniczenko et al., 2019). Accordingly, in the present research the V4 region was applied to study the protist community. PCR amplification was performed using the primer set 528F (5'-GCGGTAATTCCAGCTCCAA-3') and 706R (5'-AATCCRAGAA TTTCACCTCT-3'). Sequencing was performed bidirectionally on the Ion S5TMXL platform (Thermofisher, USA) in Novogene sequencing company (Tianjin, China).

2.3 Processing Analyses of High-Throughput Sequencing (HTS) Data

Raw data were processed and analyzed using QIIME v.1.8.0 (Caporaso et al., 2010) and Mothur v.1.34.4 (Schloss et al., 2009). Sequences with high quality were obtained after processes of quality controlling, denoising and discarding chimera and singleton.

The clean sequences were clustered to operational taxonomic units (OTU) at a sequence similarity of 97%. The OTUs were blasted, aligned, and classified against the Protist Ribosomal Reference (PR2 database; Guillouet al., 2012). Sequences of non protists (metazoan, macroalgae) and unclassified taxa were omitted (2% of OTUs). Alpha diversity of protists inferred by OTU richness, Shannon, Simpson, and Chao1 were calculated based on the lowest number of the high-quality sequences among all samples. Beta diversity was calculated with Bray-Curtis distances and was visualized using nonmetric multidimensional scaling (nMDS) in PRIMER v.6.0 (Primer-E, UK).

2.4 Network Construction

The co-occurrence network was inferred based on the spearman correlation matrix constructed following the methods of Ma et al. (2016), using the WGCNA package (Langfelder and Horvath, 2012) appended in R v3.5.1 (RCore, 2016) and packages provided by Ma et al. (2016). Only OTUs with relative abundance higher than 1% were subjected to the analyses, and the positive correlations between OTUs were considered. The P-values for multiple testing were adjusted using the Benjamini and Hochberg false discovery rate (FDR) controlling procedure (Benjamini et al., 2006). The direct correlation dependencies were distinguished using the network deconvolution method (Feizi et al., 2013). The co-occurrence networks were constructed based on correlation coefficients and FDR adjusted P-values for correlation. The cutoff of FDR-adjusted P-values was set as 0.05. The cutoff of correlation coefficients was 1.25 for station X and 1.4 for station Y through the random matrix theory-based methods (Luo et al., 2006).

2.5 Statistical Analyses

One-way ANOVA analysis was performed to examine the differences in environmental factors, α-diversity estimators and proportions of abundant (> 1%) groups at family level among four seasons. Spearman's correlation analysis was conducted to examine the relationships between α-diversity estimators or relative abundance of protistan families and environmental factors. The normality of data was tested and log-transformed when it was necessary. All these analyses were executed using SPSS v. 20.0 (SPSS, Chicago, IL, USA, Kirkpatrick and Feeney, 2012). To statistically test the variations in community structure of planktonic protists in different seasons, ANOSIM (analysis of similarity) was performed using PRIMER V.6 (Clarke and Gorley, 2006). The major environmental factors shaping the community structure of planktonic protists were determined by conducting Redundancy analysis (RDA) using packages in R v3.5.1 (RCore, 2016).

3 Result 3.1 Seasonal Variations of Environmental Factors

Generally, physicochemical properties of marine water (temperature, DO and nutrients) showed clear seasonal patterns (Table 1). The water temperature increased up to 25.4℃ in summer and decreased to a minimum of 2.5℃ in winter (ANOVA, P < 0.01). Ammonium, phosphate, silica and nitrate exhibited the contrast trend to temperature, with the highest values in winter and lowest values in summer (P < 0.02). The concentration of dissolved oxygen (DO) also changed with season significantly. The highest DO concentration was found in February (14.20 mg L⁻¹) and the lowest value was recorded in August and September (5.52 mg L⁻¹) (P = 0.03). Other chemical factors (e.g., the total organic carbon and nitrogen, C:N ratio, and N:P ratio) showed no significant difference with seasonal variation (P > 0.05). No significant regional difference between X and Y stations in environmental factors were detected (Table 1).

3.2 Overview of Sequencing Data and Community Composition

After quality filtering, a total of 1745616 reads were obtained in the 24 samples. Further on, 1231976 protistan sequences were retained after removing reads classified as metazoan and multicellular organisms. Sequence numbers of the 24 samples varied from 8680 to 76258 (Table 2).

By comparing the relative abundance of super-phylum groups between two sites among different seasons, we found that there was no statistically significant difference between two sites (T-test, P > 0.05; Table 3); however, most groups showed significant difference among seasons (ANOVA, P < 0.05). Therefore, we combined the data of two sampling sites and focused on the seasonal variation of the protozoan community.

Generally, infrakingdom Alveolata dominated in the protist communities, accounting for 54% of the total protistan sequences, and was observed in all samples (Fig. 2A). Most of the Alveolata sequences were clustered into phylum Dinoflagellata (74%) and Ciliophora (20%) (Fig. 2B). Higher abundance of Dinoflagellata (40% – 80%) was observed from April to July (spring and summer), while the lowest abundances appeared in autumn (mean 22%). Ciliophora was more abundant in winter than in other seasons (Fig. 2B). Within Dinoflagellata, the most abundant classes were Dinophyceae (50%) and Syndiniales (47%) (Fig. 2C); the most abundant orders were Dino-Group-I (attributing to 34% of Dinoflagellata), Peridiniales (20%) and Dinophyceae_X (16%) (Fig. 2D). In Ciliophora, Spirotrichea was the most representative class (Fig. 2C), in which Tintinnida (31%), Strombidiida (32%) and Choreotrichida (11%) were dominant (Fig. 2D).

Super-group Archaeplastida (22.85% ± 3.81% of the total sequences) was the second abundant group in the protist community (Fig. 2A), which was primarily represented by the phylum Chlorophyta (95%; Fig. 2B). Mamiellophyceae (45%), Trebouxiophyceae (25%) and Ulvophyceae (16%) were the most abundant classes in Chlorophyta (Fig. 2C). At the order level, Mamiellales (97%) was dominant in class Mamiellophyceae, while class Trebouxiophyceae was dominated by order Chlorellales (99%; Fig. 2D). Notably, the order Ulotrichales of Ulvophyceae was markedly dominant in March, which accounted for almost 45% of the sequences of the month (Fig. 2D).

Subkingdom Stramenopiles accounted for 12.57% of the total reads, with Ochrophyta as the major phylum (65%), followed by Sagenista (16%) and Pseudofungi (14%; Figs. 2A, B). Class Bacillariophyta was dominant in phylum Ochrophyta (70%). Labyrinthulomycetes (86%) was the most abundant class in Sagenista, followed by MAST-6 (6.2%), MAST-9 (3.9%) and MAST-7 (2.1%) (Fig. 2C). Oomycota (50%) and MAST-1 (40%) were the most abundant classes in Pseudofungi (Fig. 2C). At the order level, the most abundant groups were Bacillariophyta (46%), Chrysophyceae_X (8.5%) and Thraustochytriales (2.9%; Fig. 2D).

In addition, super-phylum Opisthokonta (4.9% ± 0.89%), Rhizaria (2.6% ± 0.32%), and Hacrobia (2.3% ± 0.56%) exhibited relatively low abundance. Other super-phyla (e.g., Ameobozoa, Apusozoa and Excavata) appeared to be minor ones (less than 1%; Fig. 2A).

3.3 Alpha Diversity of Planktonic Protists

Alpha-diversity estimators of protists varied greatly among the monthly samples, with OTU richness ranging from 548 to 1499, Simpson indices from 0.296 to 0.981, Shannon indices from 1.90 to 7.35, and Chao1 indices from 1348 to 3487 (Table 2). Simpson and Shannon indicators were significantly different among four seasons (P < 0.05, Fig. 3), while the highest value was found in autumn and the lowest value was recorded in spring. Otherwise, Chao1 and Richness indices showed no significant difference among seasons (P > 0.05, Fig. 3).

3.4 Seasonal Distribution of Planktonic Protists

Seasonality of protist community was indicated by the ordination analyses based on the Bray-Curtis dissimilarities. In the nMDS plot (Fig. 4), summer, autumn and winter samples were clearly separated from each other, while spring samples were evenly distributed between the spaces occupied by other season samples. The result of ANOSIM confirmed that samples among four seasons were significantly distinct (P < 0.05, Table 4).

The heterogeneity of protist community among four seasons was further determined by a comparison of relative abundance of the major (relative abundance > 1%) groups at family level (Fig. 5). Significantly, the relative abundances of Chlorellales_X and Suessiaceae were the highest in summer and the lowest in winter (P < 0.05). In contrary, Bathycoccaceae showed higher ratios in autumn (11.6%) and winter (11.4%) than in spring (0.9%) and summer (0.7%) (P < 0.01). Similarly, Mamiellaceae, Ceratiaceae, Cryptomonadales_X, Strombidiidae_R, Tintinnida_TIN_9, Raphid-pennate, and Dino-group-II-Clade-1 had the highest proportions in autumn or winter (P < 0.05). In spite of no statistical significance, the abundant families Dino-Group-I-Clade-1 and Ulotrichales_X had obviously higher proportions in spring and summer.

3.5 Co-Occurrence Network Among Protistan Groups

At both sites, significantly high betweenness centrality scores were observed for OTUs affiliated to two prevailing classes, the ciliate Spirotrichea and the dinoflagellate Syndiniale (Figs. 6, 7). The involved Spirotrichea comprised mainly tintinnids (Favella, Stenosemella, Tintinnopsis, and Tintinnidium) and a few species from Hypotrichia (Bistichella) and Choreotrichida (Parastrombidinopsis shimi). The top frequent co-occurred Syndiniales belonged to Dino-Group-I, with few exceptions from Dino-Group-II. Cooccurrences were also observed among classes Spirotrichea (mainly Strombidiida), Mamiellophycease (genera Micromonas, Ostreococcus, and Bathycoccus), and/or Dinophyceae (e.g., Heterocapsa rotundata) (Figs. 6, 7). Additionally, OTUs belonging to the same class Dinophyceae (blue circles) significantly co-occurred at station X (Fig. 6). These OTUs affiliated to genera Gyrodinium, Heterocapsa and one unidentified Dinophyceae species. Class Bacillariophyta (represented by Thalassiosira, Chaetoceros and one unidentified Polarcentric-Mediophyceae), one species Heterocapsa triquetra from Dictyochophycea, and several OTUs from Mast-groups and Chrysophyceae were positively co-varied (Fig. 6). Additionally, OTUs affiliated to genus Picochlorum in class Trebouxiophycease, showed positive co-variations to two Dinophyceae OTUs, affiliated to genera Symbiodinium and Biecheleriopsis (Fig. 6). At station Y (Fig. 7), we specifically observed significantly high betweenness centrality scores among groups Labyrinthulomycetes (green circles) and two OTUs from Bacillariophyta (Licomophora_paradoxa) and Ulotrichales (Ulva fasciata) (Fig. 7).

3.6 Environmental Drivers of Variation in Protist Diversity and Community Structure

Spearman's correlation coefficients between alpha diversity estimators and environmental factors were given in Table 5. Briefly, Shannon index had the strongest correlation with SiO₄²⁻ (R = 0.56, P = 0.01), followed by PO₄³⁻ (R = 0.54, P = 0.01), TOC (R = 0.52, P = 0.01), NO₂⁻ (R = 0.49, P = 0.02) and pH (R = 0.41, P = 0.05), respectively. Simpson index had the strongest correlation with the C:N ratio (R = 0.56, P = 0.01), followed by pH (R = 0.50, P = 0.01), PO₄³⁻ (R = 0.50, P = 0.01), TOC (R = 0.46, P = 0.03), and SiO₄²⁻ (R = 0.46, P = 0.01), respectively. No significant correlation between OTU richness and Chao1 index with environmental variables was observed, as shown in Table 5.

The redundancy analyses (RDA) plots (Fig. 8) suggested that temperature, PO₄³⁻ and TON were the most significant factors influencing the whole protist community structure among seasons (P < 0.05).

Furthermore, the associations between relative proportion of major families (relative abundance > 1%) and environmental factors were also explored (Table 6). The abundant families Dino-Group-I (Clade-1 and X) and Dinophycease_XX did not show significant relationship with any environmental factors (P > 0.05). Heterocapsaceae only had significant and positive correlation with temperature (R = 0.43, P < 0.05). In contrast, Chlorellales_X was positively related to temperature (R = 0.88, P < 0.01), but negatively related to DO (R = −0.56, P < 0.01) and nutrients including NO₃⁻ (R = −0.78, P < 0.01), NH₄⁺ (R = − 0.71, P < 0.01), PO₄³⁻ (R = −0.46, P < 0.05) and SiO₄²⁻ (R = −0.42, P < 0.05). Similarly, Suessiaceae correlated positively with temperature (R = 0.61, P < 0.01), and negatively with DO (R = − 0.60, P < 0.01), NO₃⁻ (R = −0.42, P < 0.05) and NH₄⁺ (R = −0.44, P < 0.05). The proportion of Bathycoccaceae primarily exhibited significant positive relationship with nutrients PO₄³⁻, SiO₄²⁻ and NO₃⁻ (P < 0.05). Most of the other groups, e.g., Ulotrichales_X, Polar-centric-Mediophyceae, Cryptomonadales_X, Tintinnida_TIN_09, and Dino-Group-II-Clade-1, presented a negative correlation with temperature and a positive correlation with nutrients (P < 0.05).

4 Discussion

Previously, the protist diversity and community composition in coastal water and sediments in the Yellow Sea were mainly studied using morphology-based methods (Song et al., 2009; Zhang and Xu, 2015). In this study, we firstly investigated the seasonal dynamic of protist community in the coastal Yellow Sea based on the HTS technology, which will broaden our insight to the dynamic and environmental response of protists in coastal ecosystems. Using HTS technology, a total of 26038 OTUs from different protistan groups were revealed, which far exceeded the number of taxa detected through morphological observation.

4.1 Dominant Groups in the Surficial Water of the Yantai Coast, Northern Yellow Sea

Alveolata was revealed as the most dominant infrakingdoms in this study, which is comparable with the clone library analyses and HTS studies on microeukaryotes in other water columns (Pernice et al., 2013; Gran-Stadniczeñ ko et al., 2018) and the surficial sediment in the north Yellow Sea (Gong et al., 2015).

Dinoflagellata was the most abundant phylum in Alveolata, dominated by classes Syndiniales Dino-Group-I and Peridiniales. Syndiniales is a group comprised exclusively of parasitic species that infect tintinnid ciliates, crustaceans, dinoflagellates and fish (Guillou et al., 2008). The high abundance of Syndiniales detected in this study corresponds well to the previous findings in several metabarcoding studies on planktonic microeukaryotes (Georges et al., 2014; Christaki et al., 2017), and supports the hypothesis that parasitic interactions could dominate plankton networks (De Vargas et al., 2015; Lima-Mendez et al., 2015). Otherwise, it should be aware that potentially high rRNA gene copy numbers and inactive spores could lead to overestimation on the proportion of Syndiniales (Not et al., 2009; Sassenhagen et al., 2019). Another dominant order Peridiniales is mainly composed of small mixotrophic species (< 200 μm), harboring a number of red tideforming taxa (Iwataki, 2008). In our samples, one OTU related to species Heterocapsa rotundata dominated the Peridiniales (98%). Consistently, H. rotundata was also recorded as one significantly abundant OTU in the surficial sediments of Yellow Sea coast (Gong et al., 2015). These findings suggest a common presence of red tideforming species in the coast of the Yellow Sea.

Archaeplastida was the second dominant super-group, mainly comprised of phylum Chlorophyta. The proportion of Chlorophyta detected in the Yellow Sea samples is much higher than that revealed in other water samples by HTS analyses (Pernice et al., 2013; Gran-Stadniczeñ ko et al., 2018), e.g., in the surface water of estuarine in Skagerrak coast of Europe (22% vs. 7.6%) (Gran-Stadniczeñko et al., 2018).

Lower proportion of subkingdom Stramenopiles was detected in our Yellow Sea coastal samples, compared to those reported previously in the estuarine samples in the Skagerrak coast (12% vs. 20%) (Gran-Stadniczeñko et al., 2018), and in the surficial sediments (12% vs. 15%) of the Yellow Sea (Gong et al., 2015). Moreover, our study revealed less proportion of Oomycota (0.8% vs. 3.6%) and Labyrinthulea (0.6% vs. 3.4%) in the surface water compared to that recorded in the sediment (Gong et al., 2015). This difference between surface water and sediment samples might attribute to the distinct ecological niches. Labyrinthulomycetes was reported playing a key role on the initial stage of the microbial chain food as an organic matter degrader (Raghukumar, 2002), while Oomycota was predicted being saprobic or parasitic on terrestrial or aquatic plants and animals (Nigrelli and Thines, 2013).

4.2 Co-Occurrence of Protist Groups

Significant correlations between OTUs in dinoflagellate Syndiniales and ciliate Spirotrichea were indicated by the co-occurrence analysis (Figs. 6, 7). The co-occurrence might reflect a parasite-host relationship. A few genera of Syndiniales Dino-Group-I were reported as common parasites in tintinnids of Spirotrichea (Coats et al., 2010; Bachvaroff et al., 2012). An increase of tintinnid hosts' population might provide breeding ground for the active parasite or dinospores of syndiniales (Coats et al., 2010). Other spirotrichs (e.g., Hypotrichia and Choreotrichida) involved in this relationship may indicate a wider range of host-parasite relationships than previously known, which need further in vivo observations to validate the in silico inferences (Figueroa et al., 2010). Other biotic interactions, such as the parasites favoring the growth of one species through the infection of its grazers or competitors, may also contribute to this ciliate-parasite co-occurrence (Coats and Heisler, 1989).

The co-occurrence of ciliates and Syndiniales in this study is different from the result of Sassenhagen et al. (2019) on Kerguelen Islands, which showed negative relationship between two tintinnid genera (Tontoniidae and Dictyocystidae) and Syndiniales Group-II-clade-12. The contradiction might attribute to the variable efficiencies that Syndiniales infect ciliates. The life cycles of the Syndiniales parasitic on ciliates encompass a bi-flagellated infective stage, and the non-infectious dinospore stage that Syndiniales may only penetrate the host cell membrane or be ingested by the host as preys (Coats et al., 2010; Jephcott et al., 2016). Indeed, Sassenhagen et al. (2019) observed high concentrations of free Syndiniales spores but relatively few amounts of Syndiniales associated with microplankton.

Remarkable co-occurrence was observed among ciliate Spirotrichea (mainly Strombidiida), Mamiellophyceae (Micromonas and Bathycoccus), and Dinophyceae (e.g., Heterocapsa rotundata). This could represent a predation relationship considering the latter two small phytoplankton are well known food resources for phagocytic ciliates (Lynn, 2008). On another hand, the co-variation might associate with the mixotrophic nature of organisms in the three groups. Strombidiids are typical plastidic ciliates capable of mixotrophic by acquiring phototrophy (Stoecker and Lavrentyev, 2018). Other ciliate species involved in this network, e.g., Cyclotrichium sp. and Askenasia sp. (Figs. 6, 7), are also known mixotrophs with plastid-retention (Stoecker and Lavrentyev, 2018). Meanwhile, the photosynthetic algae Micromonas and Heterocapsa are well known mixotrophic cases which sometimes are more bacterivory than heterotrophic flagellates (Jeong et al., 2005; Stoecker and Lavrentyev, 2018). In this sense, the cooccurrence of Spirotrichea, Mamiellophyceae, and Dinophyceae in this study could also be a result of homologous ecological niches with similar requirement of nutrition and environmental conditions (e.g., light, bacteria food).

4.3 Seasonal Variation of the Protist Community and the Environmental Drivers

Our analyses revealed a significant seasonal pattern of protist community structure (Figs. 2, 5). The Shannon and Simpson index also showed significant seasonal fluctuation (Fig. 3). The highest micro-eukaryotic diversity occurred in autumn, while the lowest one appeared in spring (Fig. 3), which is different to the commonly accepted notion that phytoplankton generally produce maxima in spring and autumn, minima in summer and winter (Wetzel, 1983). This divergence could be interpreted by the specific features in community structures and nutrient fluctuations in the studied coast of the Yellow Sea. The lowest protistan diversity in spring might attribute to the significant blooms of several families in this season, including two parasitic families of Syndiniales Dino-group-I and another family Ulotrichales_X (probably in forms of spore with diameter < 200 μm, or retentions from large foliate bodies) (Fig. 5). On another hand, the protist diversity strongly positively correlated with nutrients (Table 6), which may explain the highest diversity of these groups in autumn. Indeed, we observed higher nutrient concentrations in autumn of the year 2018 in this area (Table 1).

The family Ulotrichales_X blooming in spring was dominantly represented by one green alga Monostroma grevillei in our data (Fig. 5), which has been reported usually occurring in cold-water sublittoral in winter and spring (Sumner and Osburn, 1913; Levring et al., 2019). This is consistent to its negative correlation to temperature in the Spearman's analysis (Table 6). Moreover, the positive Spearman's correlation of Ulotrichales_X with NH₄⁺ indicates a stimulation by ammonia nitrogen to its spring blooming.

Bathycoccaceae in order Mamiellales of phylum Chlorophyta was significantly higher in autumn and winter than in spring and summer (Fig. 5). This finding agrees with several other studies showing that members of Bathycoccaceae have advantages under low-light conditions or persist throughout winter darkness in marine water (Joli et al., 2017). In this study, the distribution of this family could be mainly explained by the changes of SiO₄²⁻ and NO₃⁻ in water, but had no correlation with temperature (Table 6). The family was represented by genus Bathycoccus in our data (42%). The rich nutrients in the winter water may favor the growth of Bathycoccus, which has shown high abundance and wide distribution in nutrient-rich waters but relies less on the temperature and light time (Vannier et al., 2016).

The family Strombidiidae_R of phylum Ciliophora showed significantly higher abundance in winter than in other seasons (Fig. 5). The prevailing of Strombidiidae over other heterotrophic microzooplankton in winter may attribute to their advantages in mixotrophy (Stoecker and Lavrentyev, 2018). These mixotrophic ciliates are less restricted by the decreased food resources of small phytoplankton in winter, compared to the strictly heterotrophic microzooplankton (Mitra et al., 2016).

5 Conclusions

In this study, we obtained a total of 26038 protistan OTUs in the northern Yellow Sea coastal samples using HTS. Alveolata, Archaeplastida, and Stramenopiles are the dominant super-phylums in this coastal area. The protist community composition and relative abundance showed dramatic seasonal variation through the year. However, we did not observe significant difference between the estuary and intertidal zones, which may be caused by the close distance of the two sampling sites, and should be considered in the future sampling strategies. Blooms of the chlorophyte family Ulotrichales_X and the red tideforming taxa Heterocapsa rotundata were recorded in some seasons, which may indicate a sign of eutrophication in the coast of Yantai. The co-occurrence analysis method applied in this work turned out to be efficiency and revealed frequent and strong biotic interactions, which could be widely used in the future ecological studies on protists.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Nos. 31672251, 31772413), the Youth Innovation Promotion Association, CAS (No. 2019216), the Strategic Priority Research Program of the Chinese Academy of Sciences (No. XDA23050303), and the Key Research Project of Frontier Science, CAS (No. QYZDBSSW-DQC013-1). We specially thank Prof. Bin Ma from Zhejiang University for his kindly help in the co-occurrence network analyses.