1 INTRODUCTION

Qinhuangdao is located to the west of the Bohai Sea (BS) between the Liaodong Gulf and Bohai Gulf. With the development of social economy and acceleration of industrialization over the last two decades, several anthropogenic pollutants have been continuously discharging into the Qinhuangdao coastal waters through overland runoff (Li and Cui, 2012). Meanwhile, with the continuous expansion of mariculture activities, the increasing aquaculture wastewater discharge further aggravates environmental degradation (Zhen et al., 2016; Yao et al., 2019). These nutrient substances of anthropogenic origin lead to severe eutrophication and significant change in the nutrient composition in the coastal waters of Qinhuangdao (Zhen et al., 2016; Ou et al., 2018).

The intensification of coastal eutrophication typically results in a significant increase in the harmful algal blooms (HABs) in terms of frequency, intensity, distribution, and diversity of the causative species (Anderson et al., 2012; Glibert, 2020). A similar phenomenon has occurred in the Qinhuangdao coastal waters. Only five HAB events were recorded in this area in the 1990s, and the dominant causative species was Noctiluca scientillans. However, more than 60 HAB events caused by diversified causative species (30 species from 20 genera) have been recorded in this area from 2000 to 2020 (Ministry of Natural Resources (MNR), 1990–2020; Liang, 2012; Song et al., 2016). Corresponding with the diversification of dominant HAB causative species, there has been a tendency toward miniaturization in the cell size of the HAB causative species in the coastal waters of Qinhuangdao, including the small-size dinoflagellates, diatoms, haptophytes, raphidophytes, and pelagophytes (Ministry of Natural Resources (MNR), 1990–2020; Song et al., 2016; Li, 2021). For example, large-scale brown tides of tiny pico-phytoplankton pelagophyte Aureococcus anophagefferens outbroke in this area from 2009 to 2015, with a total blooming area of more than 10 000 km2 (Ministry of Natural Resources (MNR), 1990–2020; Zhang et al., 2012, 2021). The HAB incidents with raphidophyte Heterosigma akashiwo as the unique or one of the mixed dominant species outbroke in this area in 2007, 2009, 2015, and 2016 (Ministry of Natural Resources (MNR), 1990–2020; Li, 2021). Dinoflagellate Prorocentrum minimum triggered six HAB events in this area during the period from 2006 to 2017 (Ministry of Natural Resources (MNR), 1990–2020; Li, 2021). The diversified and miniaturized HAB causative species reflect the successional tendency of phytoplankton assemblies in the coastal waters of Qinhuangdao; however, there is still little understanding of the composition and succession of the nano- and pico-eukaryotic phytoplankton (NPEP) community in this area under the background of eutrophication.

Green tide, a type of macroalgal bloom, has been recurrently outbreaking in the Qinhuangdao coastal waters since 2015 (Han et al., 2022). Unlike the largescale green tide of macroalgae Ulva prolifera in the Yellow Sea (YS), the green tide in the Qinhuangdao coastal waters is a local bloom distributed along the beaches in the Jinmeng Bay and Beidaihe (Song et al., 2019a, b). Several suspended macroalgal species, including Ulva pertusa, Bryopsis plumosa, Gracilaria lemaneiformis, and U. prolifera, have seasonal succession from April to September (Han et al., 2019; Song et al., 2019a, b). Macroalgal blooms can lead to adverse impacts, such as shading, hypoxia, changes in composition and content of nutrients, and marine biocoenosis structure and biodiversity, on the marine environment (Lapointe et al., 1994, 2005; Valiela et al., 1997; Lyons et al., 2014). Several field studies and simulated experiments have revealed that the decomposing macroalgae have diversified effects on the different taxon of phytoplankton and change the phytoplankton assemblages (Wang et al., 2012; Zhao et al., 2022). Additionally, the phytoplankton assemblages and floating/suspended macroalgae green algae can inter-affect through phenomena such as allelopathy and nutrient and niche competition (Smith and Horne, 1988; Jin and Dong, 2003; Nan et al., 2004; Tang and Gobler, 2011; Zhang et al., 2013). Although the formation mechanism of green tides has almost been elucidated (Han et al., 2019; Song et al., 2019a, b), the potential impacts of macroalgal blooms on the phytoplankton in the Qinhuangdao coastal area need to be clarified.

During recent decades, the coastal waters of Qinhuangdao have been hitting by intensive microalgal blooms, particularly the NPEP blooms (Song et al., 2016; Li, 2021), which are typically neglected in previous studies. Thus, there is an urgent need to determine the composition of NPEP assemblages and their variations in the Qinhuangdao coastal waters with the increasing diversity and tendency of miniaturization in the cell size of the HAB species. Moreover, the frequent outbreak of green tide can certainly alter physical and chemical properties of seawater and further impact the composition of the NPEP assemblages. Therefore, an investigation was conducted in the Jinmeng Bay, where the green tides have been hitting regularly since 2015 (Han et al., 2019; Song et al., 2019a, b). The aims of the study are as follows: 1) analyzing the composition the NPEP assemblages and their variations via high-throughput sequencing (HTS) assay using the eukaryotic 18S rDNA V4 region as the target gene, 2) exploring the relationships between the variations in NPEP assemblages and the environmental factors, 3) evaluating the potential impacts of green tide on the NPEP assemblages.

2 MATERIAL AND METHOD 2.1 Sample collection

Sixteen sampling sites covering the green-tide area in Jinmeng Bay were designed to measure the biomass of the suspended macroalgae (Fig. 1). The biomass of the suspended macroalgae indicates that the green tide started in April and lasted until November 2020 (Han et al., 2022). Two representative sampling sites were selected to collect phytoplankton samples. One site, A32, was in the inshore waters, and the other site, A14, was in the offshore waters (Fig. 1). On May 14, August 23, and September 23, 2020, triplicate 1-L surface water was collected from each site and filtered through sieve silk with apertures of 20 μm, which was then re-filtered on polycarbonate membranes with apertures of 0.4 μm (HTTP, Millipore, USA) under the conditions of negative pressure less than 0.05 MPa. The filtered membranes were then stored in a freezer at -80 ℃ for DNA extraction.

2.2 Determination of environmental factors

Multiple environmental factors, including temperature, salinity, dissolved oxygen (DO), chlorophyll a (chl a), ammonium, nitrate, dissolved inorganic nitrogen (DIN), and phosphate (DIP), dissolved organic nitrogen (DON), dissolved organic phosphorus (DOP), and dissolved organic carbon (DOC), were determined at the two sampling sites. A portable multi-parameter water quality detector (YSI ProQuatro, Yellow Springs, OH, USA) was used to measure the surface seawater temperature and salinity in situ. The mass concentration of DO was measured via the Winkler method (McNeil and D'Asaro, 2014). The chl-a concentrations were measured using a calibrated 10-AU-005-CE fluorometer (Turner Designs, San Jose, Ca, USA), and the specific measurement procedure is described in Han et al. (2022). The ammonium, nitrate, and DIP concentrations were analyzed using an AutoAnalyzer (AA3; BRAN and LUEBBE, Norderstedt, Germany). Total dissolved nitrogen (TDN) and phosphate (TDP) were determined through persulfate oxidation using continuous flow auto analyzer systems (AA3; BRAN and LUEBBE, Norderstedt, Germany). The DON and DOP concentrations were calculated as the differences between TDN and DIN, and TDP and DIP, respectively (Grasshoff et al., 1999). The DOC content was determined via the high-temperature catalytic oxidation method following Yu et al. (2012). Briefly, the seawater samples were thawed at 20 ℃, acidified with 2-mol/L hydrochloric acid, and then purged with high purity oxygen to remove the inorganic carbon in the samples. Finally, DOC concentration in the water samples was measured using Multi N/C 3100 Analyzer (Analytik Jena, Germany). The biomass of suspended algae (SA) was collected, and its biomass was calculated using the method described by Han et al. (2022).

2.3 eDNA extraction, amplification, and Illumina HiSeq HTS

The eDNA was extracted using the cetyltrimethylammonium bromide (CTAB) method described by Winnepenninckx et al. (1993) with the following modifications: incubated with CTAB buffer at 65 ℃ for 1 h; centrifugal separation conditions were 4 ℃, 12 500×g; isopropyl alcohol in equal volume precipitated DNA. The extracted eDNA was stored at -80 ℃. The hypervariable V4 region of 18S rDNA from extracted eDNA was amplified with primer pair D514 (5′-TCCAGCTCCAATAGCGTA-3′) and B706R (5′-AATCCRAGAATTTCACCTCT-3′) (Cheung et al., 2010; Zimmermann et al., 2011). Polymerase Chain Reaction (PCR) were performed in a 20-μL mixture containing 10 μL of 2×Taq PCR MasterMix (TIANGEN), 1 μL of prepared eDNA as a template, 0.8 μL of each primer (10 μmol/L), and 7.4 μL of sterilized deionized water. PCR amplifications were performed with an initial denaturation temperature at 95 ℃ for 3 min, followed by 35 cycles of 95 ℃ for 20 s, 50 ℃ for 30 s, and 72 ℃ for 100 s, and final elongation at 72 ℃ for 7 min. The final reaction products were detected via 1% agarose gel electrophoresis, and the gel imaging system (BIO-RAD Gel Doc EZ Imager, Beijing Bioresources Biotechnology Co., China) images were obtained. The reaction products with targeted DNA bands were purified via quick DNA extraction using SteadyPure PCR DNA Purification Kit AG21003 (AG, China). Sequencing was performed by Novogene Bioinformatics Technology Co. Ltd. (Beijing, China) on a Novaseq 6000 platform using the Miseq reagent kit V3. Raw sequence data were deposited in the NCBI Sequence Read Archive (SRA) under accession number PRJNA809593.

2.4 Data processing and analysis

The adaptor sequence and barcode sequence were removed from both ends of the paired-end reads (raw reads), and raw reads were split into 18 samples based on unique barcodes and merged using FLASH (Magoč and Salzberg, 2011). QIIME2 (Caporaso et al., 2010) software was used to filter the merged raw tags containing more N bases or low-quality bases and chimeric sequences. High-quality sequences were finally obtained for the subsequent analysis.

Subsequently, QIIME2 software was used to conduct two analyses for the effective sequences obtained after quality control. One was DADA2 denoising analysis (Callahan et al., 2016), and the other was clustering analysis (Rognes et al., 2016). Diversity analysis was conducted on the operational taxonomic units (OTUs) and amplicon sequence variants (ASVs) obtained from corresponding analyses, respectively. The classifier method and Silva database (Quast et al., 2013) annotated the representative sequences of each OTU and ASV. The results of their species annotation effect statistics and diversity analysis were compared to select an analysis method that could effectively reflect the species diversity of the study area. The processing results of the selected analysis method were used for the following analysis.

The processed OTUs and ASVs were screened, and the OTUs and ASVs annotated as zooplankton, fungi, and macroalgae were removed from the dataset. The screened dataset was then homogenized. Rarefaction curve for sampling depth verification was mapped to evaluate whether the sequencing volume was sufficient to cover all groups in the samples. Subsequently, the OTUs and ASVs of screened samples were statistically analyzed according to sampling sites and sampling time, and alpha diversity metrics of the total species abundance index (Chao1) and species diversity index (Shannon) were calculated using QIIME2. The Kruskal-Wallis test was conducted for the differences in alpha diversity indexes between samples, which was considered significant at P < 0.05. The Spearman test was performed on the correlation analysis between alpha diversity indexes and environmental and biological factors. Additionally, community structure statistics were conducted at each taxonomic level to obtain the number of the OTUs and ASVs annotated to each taxonomic level of each sample and their proportion and the corresponding species information on species-based abundance distribution. Detrended correspondence analysis (DCA) was performed using the species information to select the appropriate model for correlation analysis among environmental factors, samples, and population and determine the critical environmental drivers affecting sample distribution.

3 RESULT 3.1 Primary data statistics and quality control

Raw sequences of the sequenced samples ranged from 180 076 to 199 788, with an average of 188 728 raw sequences for each sample. After quality control, the sequenced sample obtained an average of 178 771 high-quality and effective sequences with an average length of approximately 300 bp. The values of Q20 and Q30 were greater than 98% and 96%, respectively, implying that the sequencing quality was generally high and stable, and the sequencing data were accurate and reliable (Supplementary Table S1).

3.2 Comparison of two analysis methods

After quality control, the obtained effective sequences of each sample were analyzed via clustering and DADA2 denoising methods to obtain OTUs and ASVs, respectively. The representative sequences of OTUs and ASVs were annotated, and their annotation proportion was similar (Supplementary Table S2). Shannon indexes of the obtained data using the two analysis methods were compared (Table 1), and the indexes of ASVs were greater than those of OTUs, except for the sample A32_S, indicating that the DADA2 denoising method was more suitable to reflect the species diversity of NPEP community and was selected for the subsequent analysis.

3.3 Rarefaction analysis of samples

A total of 2 554 ASVs were obtained from the effective sequences, and 1 236 ASVs representing pico- and nano-phytoplankton with an average of 201 ASVs per sample were retained through screening to remove the zooplankton and fungi and macroalgae species information. Rarefaction analysis was conducted on the screened ASVs, and the rarefaction curves of ASVs were prepared for all the samples (Supplementary Fig.S1). The number of ASVs in each sample flattened with the increase in sequencing quantity, indicating that the sequencing depth was sufficient to reflect the species diversity in the samples.

3.4 Composition of phytoplankton assemblages and their variations

The number of ASVs ranged from 157 to 257 in the samples from the sites A32 and A14, respectively, and there was a significant difference (P < 0.05) in the alpha diversity indexes (Shannon and Chao1 indexes) between every two samples (Supplementary Table S3). Top 40 ASVs with over 80% total relative abundance were selected to compare the spatial and temporal distribution of ASVs among samples (Fig. 2; Supplementary Table S4). The results show that the ASV composition from the two sites was similar within a month, except for that in September. In detail, ASV396, ASV1167, ASV1168, and ASV253 dominated in two samples of May with a total relative abundance of more than 53.7%; ASV161, ASV275, ASV190, and ASV1221 dominated in two samples of August with total relative abundance more than 64.7%. However, ASV134 and ASV669 dominated only in sample A32_S with total relative abundance more than 62.0%, and no ASVs dominated the sample A14_S, with a maximum relative abundance of 10.3% for ASV134. These results showed that the NPEP composition had significantly seasonal variations, particularly the succession of dominant ASV in the survey region.

The Silva classifier method annotated the representative sequences of ASVs, including class, order, family, genus, and species (Supplementary Table S5). A total of 25 classes, 56 orders, 77 families, 110 genera, and 97 species were annotated. The proportions of ASVs annotated at genus and species levels were 44.7% and 17.8%, respectively.

Diatoms, dinoflagellates, chlorophytes, and cryptophytes controlled the NPEP assemblages, with the total relative abundance ranging from 88.3% to 95.8% (Fig. 3; Supplementary Table S6). Diatoms and dinoflagellates dominated the NPEP assemblage in May, with a relative abundance of 41.4% and 30.9% at site A14 and 63.8% and 16.1% at site A32, and Mediophyceae (A14_M, 21.2%; A32_M, 31.9%) and Coscinodiscophyceae (A14_M, 20.1%; A32_M, 30.9%) were the top two classes of diatoms. In August, dinoflagellates (A14_A, 60.7%; A32_A, 50.5%) and chlorophytes (A14_A, 24.4%; A32_A, 29.9%) replaced diatoms to dominate the NPEP assemblages, and Mamiellophyceae and Trebouxiophyceae were the dominant classes of chlorophytes. In September, the dominance of dinoflagellates further improved to 60.9% and 88.8% in both sites, respectively, followed by Coscinodiscophyceae (diatoms) (14.2%) and cryptophytes (5.2%) as the second dominant groups in the sites A14 and A32. Overall, the NPEP communities varied from diatom-dominated in May to dinoflagellate-dominated in August and September.

A wide range of sequences (22.5%–48.2%) could not be classified into known genera (Fig. 4a–b). Among the annotated genus, Skeletonema (Coscinodiscophyceae, diatoms) with the highest relative abundance (19.2% and 26.9% at sites A14 and A32, respectively) was followed by Pelagodinium (Dinophyceae, dinoflagellates) with relative abundance of 5.1% at site A14 and Thalassiosira (Coscinodiscophyceae, diatoms) with the relative abundance of 3.3% at site A32 in May. In August, Gyrodinium (Dinophyceae, dinoflagellates) with a high relative abundance of 40.7% and 34.3% at sites A14 and A32, respectively, dominated the phytoplankton assemblages, followed by Karlodinium (Dinophyceae, dinoflagellates) (11.8%) and Ostreococcus (Mamiellophyceae, chlorophytes) (7.0%) at site A14 and Ostreococcus and Bathycoccus (Mamiellophyceae, chlorophytes) (6.8% and 6.0%) at site A32. In September, Gonyaulax (Dinophyceae, dinoflagellates) was dominant species (10.0% and 31.1% at sites A14 and A32, respectively), accompanied by Teleaulax (Cryptophyceae) (6.7%) and Gyrodinium (Dinophyceae, dinoflagellates) (6.0%) at site A14 and Pelagodinium (Dinophyceae, dinoflagellates) (6.6%) at sites A14 and A32, respectively. The variation tendency of the NPEP communities was dominated by Skeletonema (diatoms) in May to Gyrodinium (dinoflagellates) in August and finally Gonyaulax (dinoflagellates) in September.

Only 16% of NPEP sequences (122 237 sequences/total 765 037 sequences) were annotated to 97 species at the species level, including 38 diatom species, 18 dinoflagellate species, 16 chlorophyte species, 9 cryptophytes species, and 6 haptophyte species (Supplementary Table S7). Diatoms exhibited higher species diversity than the dinoflagellate species; however, their relative abundance was less than that of dinoflagellates. Ten known harmful and toxic species were detected in NPEP communities, including Alexandrium andersonii, A. hiranoi, A. pacificum, Gonyaulax spinifera, and Karlodinium veneficum (dinoflagellates); Pseudochattonella verruculosa (dictyochophyte); A. anophagefferens (pelagophyte); Chattonella subsalsa and Heterosigma akashiwo (raphidophyte); and Phaeocystisglobosa(haptophyte). To further understand the variations in the relative abundance of the annotated species, 37 species with more than 100 sequences in at least one sample were screened, and their variations in the samples, including 14 diatom species, 8 dinoflagellate species, 5 chlorophyte species, 5 haptophyte species, and 3 cryptophyte species, were shown in Fig. 5. Notably, toxic K. veneficum (dinoflagellate) maintained a high abundance during the investigations, with a maximum relative abundance of 11.6%, and Gonyaulax fragilis (dinoflagellate) proliferated suddenly in September with a maximum relative abundance of 30.9%. The diatom species, such as Cylindrotheca closterium, Haslea feriarum, Chaetoceros muellerii, Cyclotella choctawhatcheeana, Skeletonema menzellii, Thalassiosira concaviuscula, Coscinodiscus radiatus, and Guinardia striata, had low abundance (0.0–3.4%) during the three investigations. The results showed the vast diversity of NPEP species. A few different toxic and harmful NPEP species coexisted in the Qinhuangdao coastal waters.

3.5 Relationships between phytoplankton communities and environmental factors

The temporal variations in hydrological and chemical factors were more evident than the spatial variations in the sampling sites. Sea surface temperature increased from 17.5 ℃ in May to 25.2 ℃ in August and decreased to 22.0 ℃ in September, and the salinity gradually reduced from 31.93 to 31.07. The DO was 10.5 mg/L in May, fell to 6.6 mg/L in August, and quickly increased to 11.6 mg/L in September. The ammonium (from 2.5 to 6.8 μmol/L and then to 3.8 μmol/L) and DIP (from 0.11 to 0.46 μmol/L and then to 0.25 μmol/L) exhibited a similar variation tendency with that of sea surface temperature. The nitrate (from 12.2 to 3.8 μmol/L and to 6.5 μmol/L) and DOP (from 0.30 to 0.14 μmol/L and to 0.21 μmol/L) also exhibited similar variation tendency. The DON (from 20.6 to 15.9 μmol/L and to 25.1 μmol/L) exhibited a different variation tendency with the DO. Both chl a and DOC had similar variations, increasing from 8.3 μg/L and 2.9 mg/L in May to high concentrations of 10.5 μg/L and 10.3 μg/L and then to 3.6 mg/L in August and September, respectively (Supplementary Table S8). The suspended green algae had a patchy distribution and were susceptible to wind and tidal changes, making it difficult to detect their biomass quantitatively in specific locations in Jinmeng Bay. In this study, the average biomass of suspended macroalgae at all sites was used to evaluate green tide dynamics in the survey area. The average biomass of suspended macroalgae was 156, 549, and 428 g/m³ in the seawater in May, August, and September, respectively (Supplementary Table S8).

Spearman test was conducted on the correlation between alpha diversity indexes and environmental factors and the biomass of the suspended macroalgae (Table 2). Evidently, the Shannon indexes exhibited significant positive correlations with DOP (P < 0.01), DON and DO (P < 0.05); Chao1 indexes exhibited significant positively correlations with DOP and DO (P < 0.01), and DON and nitrate (P < 0.05), but significantly negatively correlated with DIP (P < 0.05) and the biomass of suspended macroalgae (P < 0.05). Both indexes had no significant correlation with temperature, salinity, chl a, ammonium, and DOC. The results showed that the NPEP diversity in the Qinhuangdao green-tide area was primarily affected by DON, DOP, and DO. Multiple environmental factors, including DON, DOP, DO, nitrate, DIP, and suspended macroalgae, affected the NPEP abundance.

Redundancy analysis (RDA) was conducted to identify the determinant environmental factors that contributed to the species abundance of the top five genera of each sample (Fig. 6; Supplementary Table S9). The first and second RDA axes explained 46% and 42% of the total variations, respectively. The first axis was significantly positively correlated with temperature, ammonium, DIP, and suspended macroalgae and significantly negatively correlated with DO and nitrate. The second axis was significantly positively correlated with temperature and suspended macroalgae and significantly negatively correlated with salinity, nitrate, and DOP.

Diversified correlations existed within the environmental factors. Notably, the suspended macroalgae had a highly significant positive correlation with temperature, ammonium, and DIP (P < 0.01) and a significant negative correlation with DO, nitrate, and DOP (P < 0.01) (Fig. 6; Supplementary Table S10). The results showed complex interactions among the environmental factors, and the suspended macroalgae may affect the composition and content of nutrients in the Qinhuangdao green-tide area.

The top five genera belonged to five groups, including diatoms, dinoflagellates, chlorophytes, cryptophytes, and haptophytes (Supplementary Table S9). The dinoflagellate group exhibited a significant positive correlation between Gyrodinium and Karlodinium (P < 0.01). A significant positive correlation was observed in the chlorophyte group between Bathycoccus, Ostreococcus, and Micromonas (P < 0.01). Between dinoflagellates and diatoms, Gyrodinium had a highly significant negative correlation with Skeletonema and Thalassiosira (P < 0.01); Chaetoceros had a highly significant positive correlation with Gyrodinium and Karlodinium and a significant negative correlation with Pelagodinium (P < 0.01). Between dinoflagellates and chlorophytes, Gyrodinium and Pelagodinium had significant positive and negative correlations, respectively, with Bathycoccus, Ostreococcus, and Micromonas (P < 0.01) (Fig. 6; Supplementary Table S10). The results indicated the complex interactions among the different groups of phytoplankton communities in the Qinhuangdao green-tide area.

Between the phytoplankton abundance and environmental factors, Skeletonema had a highly significant positive correlation with salinity and nitrate (P < 0.01) and a highly significant negative correlation with temperature, suspended macroalgae, DIP, and DOC (P < 0.01); Gyrodinium and Karlodinium (Dinophyceae, dinoflagellates) had a highly significant positive correlation with temperature, ammonium, DIP, and suspended macroalgae (P < 0.01) and had a highly significant negative correlation with DO and DOP (P < 0.01). Green algae Bathycoccus, Ostreococcus, and Micromonas had (a highly) significant positive correlation with temperature, ammonium, DIP, and the suspended macroalgae (P < 0.01 and P < 0.05) and had (a highly) significant negative correlation with nitrate (P < 0.01 and P < 0.05) (Fig. 6; Supplementary Table S10). These results showed that diatoms preferred nitrate (oxidized form of N) and low temperature, and dinoflagellates and chlorophytes preferred ammonium (reduced form of N) and high temperature. Notably, the suspended macroalgae had multiple effects on different groups of phytoplankton communities, such as inhibiting diatoms and promoting dinoflagellates and chlorophytes.

4 DISCUSSION 4.1 Composition of nano- and pico-eukaryotic phytoplankton assemblages and their variations in the Qinhuangdao green-tide area

The wide application of amplicon HTS technology has greatly promoted the study of plankton diversity, particularly for the indistinguishable nano- and picoplankton with a microscope (Howard et al., 2006; Amaral-Zettler et al., 2010; Wang et al., 2015). The OTU clustering, the common method to reveal the biodiversity analysis and species composition, is a traditional analysis method of HTS data (Westcott and Schloss, 2015). Despite its high false positives and negatives, the OTU clustering method is reported as ineffective in comparing species composition in some samples (Rosen et al., 2012; Callahan et al., 2016, 2017). The ASV analysis method can capture variations in all samples and is proposed to replace OTU clustering in the comparative analysis of species composition of different samples (Callahan et al., 2017). This study used OTU clustering and ASV assays to process the HTS data of phytoplankton diversity in the survey area. The results showed that the abilities of species annotation were similar between the two assays. However, the Shannon diversity index obtained by the ASV assay was greater than that of the OTU clustering method (Table 1), indicating that the ASV assay can efficiently reflect the NPEP diversity in the Qinhuangdao green-tide area. Therefore, the ASV assay was recommended to analyze the HTS data of NPEP in the Qinhuangdao green-tide area.

In this study, based on the ASV assay results, 25 classes, 110 genera, and 97 species of NPEP were identified and annotated in the Qinhuangdao greentide area. In a previous study, 18 groups at the class level were determined in the Qinhuangdao brown-tide area from April to August in 2013 and 2014, based on amplicon HTS assay of the 18s rDNA V4 region (Chen et al., 2019). Compared with Chen's study, this study found additional species of eight classes, namely Prasinodermophyceae, Ulvophyceae, Bolidophyceae, Synurophyceae, Phaeophyceae, Chrysomerophyceae, Pavlovophyceae, and Katablepharidophyceae; however, Pinguiophyceae was missing. According to the composition of phytoplankton assemblages and their variations at the different taxonomic levels during the three investigations, the composition of the NPEP communities exhibited a seasonal succession from diatom-dominated to dinoflagellate-dominated (Figs. 3–4). In August and September, the abundance of dinoflagellates was more than 50% and even reached 88.8% in the sample A32_S (Fig. 3). Actually, the previous studies have been recorded the increase of dinoflagellate abundance in the NPEP community and micro-sized phytoplankton community in the Qinhuangdao coastal waters. Dinophyceae (dinoflagellates) represented approximately 50% of OTU richness among the micro- to picophytoplankton (cell size between 0.68 μm to 200 μm) from April to October 2013 (Xu et al., 2017a). The shift from diatom-dominated to flagellate-dominated in the phytoplankton community has appeared in many coastal regions (Fu et al., 2012; Glibert et al., 2014). For example, during the last decades, the evolution of phytoplankton communities has been characterized by a decreasing dominance of diatoms and an increasing contribution of dinoflagellates in the YS (Li et al., 2021). Multiple environmental pressures and biological factors, including the growing nutrient pollution and aquaculture, alteration in nutrient composition in the coastal waters, climate change, and ecological adaptation of different phyla of phytoplankton have been proposed to be responsible for these phenomena (Anderson et al., 2012; Zhou et al., 2017, Glibert, 2020).

The proportions of ASVs annotated at the levels of genus and species were 44.7% and 17.8% (Supplementary Table S5), implying that massive unknown NPEP species exist in the Qinhuangdao coastal waters, which cannot be determined via the morphological detection method through microscopy. For example, 123 species of phytoplankton (cell size of 0.68–200 μm) were identified via microscopy in the monthly samples in the Qinhuangdao coastal waters from March 2013 to January 2014; however, there were 7 NPEP species in the 34 dominant species (Xu et al., 2017b). The studies suggest the abundant species diversity of NPEP and the massive unknown NPEP species in the Qinhuangdao coastal waters.

4.2 Response of nano- and pico-eukaryotic phytoplankton to the environmental factors in the Qinhuangdao green-tide area

The environmental factors considerably changed from May to September and exhibited complex interactive relationships (Supplementary Table S10). The RDA revealed that the first axis was significantly positively correlated with temperature, ammonium, DIP, and the suspended macroalgae and significantly negatively correlated with DO and nitrate. The second axis was significantly positively correlated with temperature and suspended macroalgae and significantly negatively correlated with salinity, nitrate, and DOP (Fig. 6). The results suggested that they are the key environmental factors that impact the species diversity of NPEP in this area. As the diatoms and dinoflagellates are the two dominant species groups in the NPEP assemblages in the Qinhuangdao coastal waters, their variations with the changes in the environmental factors is further discussed as follows.

Our results show that the variations in the phytoplankton composition reflected the different ecological adaptability of the NPEP groups to the variated environments, particularly to the temperature (Fig. 6; Supplementary Table S10). During the period, diatoms dominated the NPEP assemblages in May when the surface seawater temperature was 17.5 ℃; however, dinoflagellates were dominant groups in August and September when the surface seawater temperature was more than 22 ℃, exhibiting a seasonal succession from diatom-dominated to dinoflagellate-dominated (Figs. 3–4; Supplementary Tables S6 & S8). Diatoms are considered taxa with low temperature optima; however, dinoflagellates prefer warmer temperature conditions (Anderson, 2000; Paerl and Huisman, 2008; Paerl and Scott, 2010). It is a normal succession process thorough which the dominance of diatoms gradually decreases in eukaryotic phytoplankton assemblages with the increase in seasonal temperature. Moreover, seawater temperature increased more than 2 ℃ in the BS in the last decades (Sun et al., 2019). The drastic increase in seawater temperature shortens the optimum growth of diatoms (Gobler, 2020). In addition, the high temperature increases seawater stratification and shallows mixed-layer, reducing the vertical exchange rate of nutrient cycling (Wang et al., 2021) and compressing the niches for no-flagellate species (Hallegraeff, 2010; Xiao et al., 2018; Gobler, 2020).

In addition to seawater temperature, the succession of phytoplankton communities from diatom-dominated to dinoflagellate-dominant is affected by the composition and content of nutrients in the Qinhuangdao coastal waters. Silicate is an essential nutrient for the growth of diatoms. However, this study did not include silicates owing to certain condition limitations. Our results showed that the NPEP diversity (Shannon indexes) had a significant positive correlation with DON and DOP in the Qinhuangdao green-tide area (Table 2), implying that DON and DOP contents are closely related to the species diversity of NPEP in this area. Dinoflagellates, chlorophytes, and cryptophytes were three dominant phyla with relative abundance ranging from 51.5% to 90.1% in the samples, except for the 29.45% in the sample A32_M. They have been demonstrated to have various trophic modes, including strong abilities on the utilization of the reduced form of nitrogen and DOP and mixotrophic capabilities (Whitney and Lomas, 2016; Yoo et al., 2017; Glibert, 2020), which may be the direct reason for the increase in species diversity as a result of the high level of DON and DOP in this area. Meanwhile, the study revealed that the NPEP abundance (Chao1 indexes) had a significant positive correlation with multiple factors, including nitrate, DON, DOP, and DO and significant negatively correlated with DIP and suspended macroalgae. The DIN/DIP ratios were 133, 23, and 41 in May, August, and September, respectively. Except for that, at the bloom peak of the green tide in August, DIN/DIP ratios were considerably greater than the Redfield N/P ratio (16:1) (Redfield et al., 1963), implying that the Qinhuangdao coastal waters are DIN-rich and DIP-limited, consistent with several previous studies (Ou et al., 2018; Yao et al., 2019; Zhang et al., 2021). Thus, the DIP is deduced to be a potential environmental constraint for the NPEP abundance in this area.

In this study, Skeletonema and Thalassiosira, the dominant diatom species, had different positive relationships with nitrate but negative relationships with ammonium (Supplementary Table S10). The results are consistent with the previous understanding that diatoms are considered nitrate opportunists and are less adept at using the reduced N forms than dinoflagellates (Glibert, 2020). Coastal eutrophication has occurred widely on a global scale (Glibert, 2020; Wang et al., 2021), and the continual increase in the content of the reduced forms of nitrogen (ammonium and DON) and inorganic N/P ratio virtually changed the nutrient composition and stoichiometry (Ning et al., 2009; Wang et al., 2021). Evidence has indicated that flagellates are more advantageous in such an environment than diatoms (Rhee and Gotham, 1980; Glibert and Burkholder, 2011; Xiao et al., 2018). During the recent 30 years, the increase in the input of land-based pollutants and rapid mariculture development has resulted in a sharp increase in nutrients and inorganic N/P ratio in the Qinhuangdao coastal waters (Zhen et al., 2016; Zhang et al., 2021). Although the land-based sewage discharge and mariculture scale began to decrease after the occurrence of A. anophagefferens brown tide from 2009 to 2015, the eutrophication featuring a highly reduced form of N and inorganic N/P ratio is still severe in this area (Ou et al., unpublished data).

Briefly, the temperature and inorganic N/P ratio have different impacts on diatoms and dinoflagellates, the two dominant NPEP groups in this area. The high temperature and inorganic N/P ratio seem to favor the dinoflagellates (Glibert, 2020; Zhang et al., 2022), resulting in the increasing dominance of dinoflagellates in phytoplankton community as an ecological response to the eutrophication and ocean warming in the Qinhuangdao coastal waters.

4.3 Potential impacts of green tides on the nano-and pico-eukaryotic phytoplankton in the coastal waters of Qinhuangdao

Macroalgal blooms, in terms of frequency, scale, and diversity of causative species, increase in many coastal waters and have become a marine ecological problem worldwide (Smetacek and Zingone, 2013). Since 2015, the green tides have been recurrently breaking out in Jinmeng Bay and Beidaihe of Qinhuangdao. Several suspended macroalgal species of Ulva, Bryopsis, and Gracilaria have seasonal succession from April to September (Han et al., 2019; Song et al., 2019a, b). Owing to the frequent occurrence and increasing scale, the impacts of macroalgal blooms are of great concern. Reportedly, the occurrence of macroalgal blooms not only has detrimental impacts on aquaculture, fisheries, and the tourism-based economy but also exerts a negative effect on the marine environments, particularly bringing changes in composition and stoichiometry of nutrients and biodiversity of marine benthos and plankton communities (Lyons et al., 2014).

To further demonstrate the potential impacts of the green tide in Jinmeng Bay, the RDA of the suspended macroalgae biomass with other environmental factors was conducted. The results showed that the suspended macroalgae had different positive relationships with ammonium, DIP, and DOC and negative relationships with nitrate, DON, and DO (Fig. 6; Supplementary Table S10), implying the complex interaction between the suspended macroalgae and nutrients. As a primary producer, macroalgae can utilize multiple forms of nitrogen and phosphorus in the seawater, implying the suspended macroalgae can compete for the nutrients with phytoplankton, particularly the potential limited nutrient. For example, the massive floating macroalgae uptake of the increased nutrients accounts for the decrease in phytoplankton biomass in the YS in June based on the observation of longterm chl-a changes (Xing et al., 2015). However, the preference for the nutrient absorption and utilization of different groups of microalgae change the nutrient composition in the seawater. An example is that U. prolifera rapidly absorbs and consumes the DIN, DIP, and DOP in the seawater and releases DON into the seawater (Ding, 2014), causing changes in the composition and stoichiometry of nutrients in seawater. Furthermore, during the decline of the green tides and the succession of different groups of macroalgae, the decomposition of the suspended and floating macroalgae releases various forms of carbon, nitrogen, phosphorus, and other biogenic elements (Qi et al., 2016; Chen et al., 2020), which proliferate the growth of several opportunistic microalgae and make changes in the phytoplankton community (Wang et al., 2012; Zhao et al., 2022). Zhao et al. (2022) found that the relative abundance of Pelagophyceae, Mamiellophyceae, Chloropicophyceae, and Bacillariophyceae species significantly increased in the settling region of massive floating macroalgae in the YS. Moreover, Wang et al. (2012) found that the decomposing green algae promoted the proliferation of the Heterosigma akashiwo (raphidophyte) but inhibited the growth of the diatom Skeletonema costatum. Briefly, the outbreak and decomposition of macroalgae change the composition and content of nutrients and then affect the composition of the NPEP communities in the Qinhuangdao green-tide area.

Meanwhile, the RDA on the suspended macroalgae with the annotated species revealed that the suspended macroalgae had significant positive correlations with Gyrodinium and Karlodinium (dinoflagellates) and Bathycoccus, Ostreococcus, and Micromonas (chlorophytes) but significant negative correlations with Skeletonema (diatom), Pelagodinium (dinoflagellate), Pyramimona (chlorophyte), and Chrysochromulina (haptophyte) (Fig. 6; Supplementary Table S10). Several previous studies generally considered that there are allelopathic interactions between macroalgae and phytoplankton. The macroalgae have adverse effects on diatoms; however, their effects on the dinoflagellates are controversial: some experiments found that macroalgae inhibit the growth of dinoflagellates, and the others found that macroalgae promote the development of dinoflagellates (Smith and Horne, 1988; Nan et al., 2004; Huo et al., 2010; Tang and Gobler, 2011; Wang et al., 2013; Zhang et al., 2013; Liu et al., 2018). In this study, Skeletonema was the dominant genus with a relative abundance of more than 20% in the NPEP assemblages in May, Gyrodinium and Karlodinium became the dominant genera with high relative abundance (7.8% and 37.5%) in August, and Gonyaulax (dinoflagellate) with the high relative abundance of 20.6% dominated the NPEP assemblages in September (Fig. 3). The results implied that the occurrence of green tide might accelerate the succession of dominant species in the NPEP community from diatoms to dinoflagellates.

Briefly, macroalgal blooms in the Jinmeng Bay altered the composition and content of nutrients. It exhibited complex interactions with phytoplankton, suggesting an essential but even more unpredictable impact of green tides on the NPEP community in the coastal waters of Qinhuangdao.

5 CONCLUSION

In this study, the composition and variations in the NPEP assemblages were studied by the ASV assay based on the amplicon HTS data in the Qinhuangdao green-tide area during the outbreak of macroalgae bloom from May to September 2020. Although 25 classes, 110 genera, and 97 species of NPEP, including ten toxic and harmful species, were identified and annotated, the NPEP species diversity was massive, considerably more significant than previously viewed in this area. The composition of the NPEP communities exhibited a seasonal succession from diatom-dominated to dinoflagellate-dominated from May to September, and Skeletonema, Karlodinium, and Gonyaulax were the most dominant genera in May, August, and September, respectively. High species diversity and abundance of the NPEP communities were caused by the high DON and DOP contents but might be inhibited by low DIP content in this area. The occurrence of macroalgal blooms may accelerate the variations in the nutrient composition and stoichiometry and accelerate the change in the NPEP communities from diatom-dominated to dinoflagellate-dominated, under the background of seasonal increase in seawater temperature and eutrophication featuring high DON content and high DIN/DIP ratio in the Qinhuangdao coastal waters.

6 DATA AVAILABILITY STATEMENT

The datasets analyzed during this study are available from the corresponding author on reasonable request.

7 ACKNOWLEDGMENT

We acknowledge all the researchers who worked on survey cruises to collect field samples and Oceanographic Data Center, IOCAS, for molecular data analysis.

Electronic supplementary material

Supplementary material (Supplementary Tables S1–S10 and Fig.S1) is available in the online version of this article at https://doi.org/10.1007/s00343-022-2198-7.