Variations in nano- and pico-eukaryotic phytoplankton assemblages in the Qinhuangdao green-tide area*

Weiqian ZHANG , Hongbin HAN , , Limei QIU , Chao LIU ,Qingchun ZHANG , , Guizhong ZHOU ,

1 Qingdao University of Science and Technology, Qingdao 266071, China

2 CAS and Shandong Province Key Laboratory of Experimental Marine Biology, Institute of Oceanology, CAS Center for Ocean Mega-Science, Chinese Academy of Sciences, Qingdao 266071, China

3 Laboratory for Marine Fisheries Science and Food Production Processes, Pilot National Laboratory for Marine Science and Technology (Qingdao), Qingdao 266237, China

4 CAS Key Laboratory of Marine Ecology and Environmental Sciences, Institute of Oceanology, Chinese Academy of Sciences,Qingdao 266071, China

5 Laboratory for Marine Ecology and Environmental Science, Pilot National Laboratory for Marine Science and Technology(Qingdao), Qingdao 266237, China

6 Key Laboratory of Marine Eco-Environmental Science and Technology, First Institute of Oceanography, Ministry of Natural Resources, Qingdao 266061, China

7 University of Chinese Academy of Sciences, Beijing 100049, China

8 Center of Ocean Mega Science, Chinese Academy of Sciences, Qingdao 266071, China

Abstract Qinhuangdao coastal waters have been frequently hitting by nano- and pico-eukaryotic phytoplankton (NPEP) blooms and green tides (macroalgal blooms) in the recent decade. However,understanding about the impacts of environmental factors and the green tides on the NPEP assemblages in this area is limited. In this study, the composition of NPEP assemblages and their variations were analyzed via amplicon sequence variants (ASVs) assay based on amplicon high-throughput sequencing data with the 18S V4 region as a targeted gene in the Qinhuangdao green-tide area during the green tide. Consequently, average NPEP eff ective sequences and ASVs of 178 000 and 200 were obtained from each sample, respectively.Although there were 25 classes, 110 genera, and 97 species of NPEP were identif ied and annotated, the proportions of annotated ASVs at genus and species levels were only 44.7% and 17.8%, respectively. The NPEP communities had a seasonal succession from diatom-dominated to dinof lagellate-dominated. During the three investigations, Skeletonema, Karlodinium, and Gonyaulax were the most dominant genera in May,August, and September, respectively. Species diversity and the abundance of NPEP communities could be increased by a high content of dissolved organic nitrogen (DON) and dissolved organic phosphorus(DOP) but inhibited by low dissolved inorganic phosphorus content. The outbreak of green tides could alter the composition and content of nutrients and accelerate the succession of the NPEP communities from diatom-dominated to dinof lagellate-dominated under the background of a seasonal increase in seawater temperature. These results preliminarily revealed the impacts of the recurrent occurrences of green tides on the NPEP assemblages in the Qinhuangdao green-tide area exhibiting high DON content and dissolved inorganic nitrogen/phosphorus ratio.

Keyword: nano- and pico-eukaryotic phytoplankton; high-throughput sequencing; green tide; eutrophication;Qinhuangdao

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 signif icant change in the nutrient composition in the coastal waters of Qinhuangdao (Zhen et al., 2016; Ou et al.,2018).

The intensif ication of coastal eutrophication typically results in a signif icant 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 f ive HAB events were recorded in this area in the 1990s, and the dominant causative species wasNoctilucascientillans. However, more than 60 HAB events caused by diversif ied 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 diversif ication 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 dinof lagellates, 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 pelagophyteAureococcusanophageff erensoutbroke 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 raphidophyteHeterosigmaakashiwoas 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). Dinof lagellateProrocentrumminimumtriggered six HAB events in this area during the period from 2006 to 2017 (Ministry of Natural Resources(MNR), 1990–2020; Li, 2021). The diversif ied and miniaturized HAB causative species ref lect 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 macroalgaeUlvaproliferain 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, includingUlvapertusa,Bryopsisplumosa,Gracilarialemaneiformis, andU.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 f ield studies and simulated experiments have revealed that the decomposing macroalgae have diversif ied eff ects on the diff erent taxon of phytoplankton and change the phytoplankton assemblages (Wang et al., 2012;Zhao et al., 2022). Additionally, the phytoplankton assemblages and f loating/suspended macroalgae green algae can inter-aff ect 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 clarif ied.

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.

Fig.1 Map of the Bohai Sea (a) and the sampling sites (b)

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 off shore waters (Fig.1). On May 14,August 23, and September 23, 2020, triplicate 1-L surface water was collected from each site and f iltered through sieve silk with apertures of 20 μm, which was then re-f iltered on polycarbonate membranes with apertures of 0.4 μm (HTTP, Millipore, USA)under the conditions of negative pressure less than 0.05 MPa. The f iltered membranes were then stored in a freezer at -80 °C for DNA extraction.

2.2 Determination of environmental factors

Multiple environmental factors, including temperature, salinity, dissolved oxygen (DO),chlorophylla(chla), 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-aconcentrations were measured using a calibrated 10-AU-005-CE f luorometer (Turner Designs, San Jose, Ca, USA),and the specif ic 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 f low auto analyzer systems (AA3;BRAN and LUEBBE, Norderstedt, Germany). The DON and DOP concentrations were calculated as the diff erences 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). Brief ly,the seawater samples were thawed at 20 °C, acidif ied 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, amplif ication, and Illumina HiSeq HTS

The eDNA was extracted using the cetyltrimethylammonium bromide (CTAB) method described by Winnepenninckx et al. (1993) with the following modif ications: incubated with CTAB buff er at 65 °C for 1 h; centrifugal separation conditions were 4 °C, 12 500×g; isopropyl alcohol in equal volume precipitated DNA. The extracted eDNA was stored at -80 °C. The hypervariable V4 region of 18S rDNA from extracted eDNA was amplif ied 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 amplif ications were performed with an initial denaturation temperature at 95 °C for 3 min, followed by 35 cycles of 95 °C for 20 s, 50 °C for 30 s, and 72 °C for 100 s, and f inal elongation at 72 °C for 7 min. The f inal 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 purif ied via quick DNA extraction using SteadyPure PCR DNA Purif ication 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 f ilter the merged raw tags containing more N bases or low-quality bases and chimeric sequences. High-quality sequences were f inally obtained for the subsequent analysis.

Subsequently, QIIME2 software was used to conduct two analyses for the eff ective 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 classif ier method and Silva database (Quast et al., 2013) annotated the representative sequences of each OTU and ASV. The results of their species annotation eff ect statistics and diversity analysis were compared to select an analysis method that could eff ectively ref lect 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 verif ication was mapped to evaluate whether the sequencing volume was suffi cient 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 diff erences in alpha diversity indexes between samples, which was considered signif icant atP<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 aff ecting sample distribution.

Table 1 Shannon indexes of the samples based on OTUs and ASVs data

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 eff ective 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 eff ective 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 ref lect 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 eff ective 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 f lattened with the increase in sequencing quantity, indicating that the sequencing depth was suffi cient to ref lect 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 signif icant diff erence (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 signif icantly seasonal variations,particularly the succession of dominant ASV in the survey region.

The Silva classif ier 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, dinof lagellates, 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 dinof lagellates 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,dinof lagellates (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 dinof lagellates 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 dinof lagellate-dominated in August and September.

A wide range of sequences (22.5%–48.2%) could not be classif ied 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 byPelagodinium(Dinophyceae,dinof lagellates) with relative abundance of 5.1% at site A14 andThalassiosira(Coscinodiscophyceae,diatoms) with the relative abundance of 3.3% at site A32 in May. In August,Gyrodinium(Dinophyceae,dinof lagellates) with a high relative abundance of 40.7% and 34.3% at sites A14 and A32, respectively,dominated the phytoplankton assemblages, followed byKarlodinium(Dinophyceae, dinof lagellates)(11.8%) andOstreococcus(Mamiellophyceae,chlorophytes) (7.0%) at site A14 andOstreococcusandBathycoccus(Mamiellophyceae, chlorophytes)(6.8% and 6.0%) at site A32. In September,Gonyaulax(Dinophyceae, dinof lagellates) was dominant species(10.0% and 31.1% at sites A14 and A32, respectively),accompanied byTeleaulax(Cryptophyceae) (6.7%)andGyrodinium(Dinophyceae, dinof lagellates)(6.0%) at site A14 andPelagodinium(Dinophyceae,dinof lagellates) (6.6%) at sites A14 and A32,respectively. The variation tendency of the NPEP communities was dominated bySkeletonema(diatoms)in May toGyrodinium(dinof lagellates) in August and f inallyGonyaulax(dinof lagellates) in September.

Fig.3 Composition of the six phytoplankton assemblages and their variations based on the relative abundance at the class level

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 dinof lagellate species, 16 chlorophyte species,9 cryptophytes species, and 6 haptophyte species(Supplementary Table S7). Diatoms exhibited higher species diversity than the dinof lagellate species;however, their relative abundance was less than that of dinof lagellates. Ten known harmful and toxic species were detected in NPEP communities, includingAlexandriumandersonii,A.hiranoi,A.pacif icum,Gonyaulaxspinifera, andKarlodiniumvenef icum(dinof lagellates);Pseudochattonellaverruculosa(dictyochophyte);A.anophageff erens(pelagophyte);ChattonellasubsalsaandHeterosigmaakashiwo(raphidophyte); andPhaeocystisglobosa(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 dinof lagellate species,5 chlorophyte species, 5 haptophyte species, and 3 cryptophyte species, were shown in Fig.5. Notably,toxicK.venef icum(dinof lagellate) maintained a high abundance during the investigations, with a maximum relative abundance of 11.6%, andGonyaulaxfragilis(dinof lagellate) proliferated suddenly in September with a maximum relative abundance of 30.9%. The diatom species, such asCylindrothecaclosterium,Hasleaferiarum,Chaetocerosmuellerii,Cyclotellachoctawhatcheeana,Skeletonemamenzellii,Thalassiosiraconcaviuscula,Coscinodiscusradiatus,andGuinardiastriata, had low abundance (0.0–3.4%)during the three investigations. The results showed the vast diversity of NPEP species. A few diff erent toxic and harmful NPEP species coexisted in the Qinhuangdao coastal waters.

Fig.4 Composition of the six phytoplankton assemblages and their variations based on the relative abundance at the genus level

Fig.5 Sequence abundances of the top 37 most abundant annotated species with more than 100 sequences in at least one sample were included

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 °C in May to 25.2 °C in August and decreased to 22.0 °C in September, and the salinity gradually reduced from31.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 diff erent variation tendency with the DO. Both chlaand 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 diffi cult to detect their biomass quantitatively in specif ic 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/m3in the seawater in May,August, and September, respectively (Supplementary Table S8).

Table 2 Correlation relationships between alpha diversity, environmental factors, and the suspended macroalgae

Fig.6 Redundancy analysis ordination plots of the relationships between the top f ive genera,environmental factors, and suspended macroalgae

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 signif icant positive correlations with DOP (P< 0.01),DON and DO (P< 0.05); Chao1 indexes exhibited signif icant positively correlations with DOP and DO (P<0.01), and DON and nitrate (P< 0.05), but signif icantly negatively correlated with DIP (P< 0.05)and the biomass of suspended macroalgae (P< 0.05).Both indexes had no signif icant correlation with temperature, salinity, chla, ammonium, and DOC.The results showed that the NPEP diversity in the Qinhuangdao green-tide area was primarily aff ected by DON, DOP, and DO. Multiple environmental factors, including DON, DOP, DO, nitrate, DIP,and suspended macroalgae, aff ected the NPEP abundance.

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

Diversif ied correlations existed within the environmental factors. Notably, the suspended macroalgae had a highly signif icant positive correlation with temperature, ammonium, and DIP(P< 0.01) and a signif icant 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 aff ect the composition and content of nutrients in the Qinhuangdao green-tide area.

The top f ive genera belonged to f ive groups,including diatoms, dinof lagellates, chlorophytes,cryptophytes, and haptophytes (Supplementary Table S9). The dinof lagellate group exhibited a signif icant positive correlation betweenGyrodiniumandKarlodinium(P< 0.01). A signif icant positive correlation was observed in the chlorophyte group betweenBathycoccus,Ostreococcus, andMicromonas(P< 0.01). Between dinof lagellates and diatoms,Gyrodiniumhad a highly signif icant negative correlation withSkeletonemaandThalassiosira(P< 0.01);Chaetoceroshad a highly signif icant positive correlation withGyrodiniumandKarlodiniumand a signif icant negative correlation withPelagodinium(P< 0.01). Between dinof lagellates and chlorophytes,GyrodiniumandPelagodiniumhad signif icant positive and negative correlations, respectively, withBathycoccus,Ostreococcus, andMicromonas(P< 0.01) (Fig.6;Supplementary Table S10). The results indicated the complex interactions among the diff erent groups of phytoplankton communities in the Qinhuangdao green-tide area.

Between the phytoplankton abundance and environmental factors,Skeletonemahad a highly signif icant positive correlation with salinity and nitrate (P< 0.01) and a highly signif icant negative correlation with temperature, suspended macroalgae, DIP, and DOC (P< 0.01);GyrodiniumandKarlodinium(Dinophyceae, dinof lagellates)had a highly signif icant positive correlation with temperature, ammonium, DIP, and suspended macroalgae (P< 0.01) and had a highly signif icant negative correlation with DO and DOP (P< 0.01).Green algaeBathycoccus,Ostreococcus, andMicromonashad (a highly) signif icant positive correlation with temperature, ammonium, DIP, and the suspended macroalgae (P<0.01 andP< 0.05) and had (a highly) signif icant negative correlation with nitrate (P<0.01 andP< 0.05) (Fig.6; Supplementary Table S10). These results showed that diatoms preferred nitrate (oxidized form of N) and low temperature, and dinof lagellates and chlorophytes preferred ammonium (reduced form of N) and high temperature. Notably, the suspended macroalgae had multiple eff ects on diff erent groups of phytoplankton communities, such as inhibiting diatoms and promoting dinof lagellates 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 ineff ective 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 diff erent 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 effi ciently ref lect 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 identif ied 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 diff erent taxonomic levels during the three investigations, the composition of the NPEP communities exhibited a seasonal succession from diatom-dominated to dinof lagellate-dominated(Figs.3–4). In August and September, the abundance of dinof lagellates 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 dinof lagellate abundance in the NPEP community and micro-sized phytoplankton community in the Qinhuangdao coastal waters. Dinophyceae(dinof lagellates) 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 f lagellate-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 dinof lagellates 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 diff erent 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 identif ied 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 f irst axis was signif icantly positively correlated with temperature, ammonium,DIP, and the suspended macroalgae and signif icantly negatively correlated with DO and nitrate. The second axis was signif icantly positively correlated with temperature and suspended macroalgae and signif icantly 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 dinof lagellates 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 ref lected the diff erent 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 °C;however, dinof lagellates were dominant groups in August and September when the surface seawater temperature was more than 22 °C, exhibiting a seasonal succession from diatom-dominated to dinof lagellate-dominated (Figs.3–4; Supplementary Tables S6 & S8). Diatoms are considered taxa with low temperature optima; however, dinof lagellates 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 °C 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 stratif ication and shallows mixed-layer, reducing the vertical exchange rate of nutrient cycling (Wang et al., 2021)and compressing the niches for no-f lagellate species(Hallegraeff , 2010; Xiao et al., 2018; Gobler, 2020).

In addition to seawater temperature, the succession of phytoplankton communities from diatomdominated to dinof lagellate-dominant is aff ected 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 signif icant 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. Dinof lagellates,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 signif icant positive correlation with multiple factors, including nitrate, DON, DOP, and DO and signif icant 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 Redf ield N/P ratio (16:1) (Redf ield et al., 1963), implying that the Qinhuangdao coastal waters are DIN-rich and DIPlimited, 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,SkeletonemaandThalassiosira,the dominant diatom species, had diff erent 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 dinof lagellates (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 f lagellates 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 ofA.anophageff erensbrown 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).

Brief ly, the temperature and inorganic N/P ratio have diff erent impacts on diatoms and dinof lagellates,the two dominant NPEP groups in this area. The high temperature and inorganic N/P ratio seem to favor the dinof lagellates (Glibert, 2020; Zhang et al., 2022), resulting in the increasing dominance of dinof lagellates 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 nanoand 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 ofUlva,Bryopsis, andGracilariahave 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, f isheries, and the tourism-based economy but also exerts a negative eff ect 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 diff erent 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 f loating 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-achanges (Xing et al., 2015). However, the preference for the nutrient absorption and utilization of diff erent groups of microalgae change the nutrient composition in the seawater. An example is thatU.proliferarapidly 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 diff erent groups of macroalgae, the decomposition of the suspended and f loating 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 signif icantly increased in the settling region of massive f loating macroalgae in the YS. Moreover,Wang et al. (2012) found that the decomposing green algae promoted the proliferation of theHeterosigmaakashiwo(raphidophyte) but inhibited the growth of the diatomSkeletonemacostatum. Brief ly, the outbreak and decomposition of macroalgae change the composition and content of nutrients and then aff ect 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 signif icant positive correlations withGyrodiniumandKarlodinium(dinof lagellates) andBathycoccus,Ostreococcus,andMicromonas(chlorophytes) but signif icant negative correlations withSkeletonema(diatom),Pelagodinium(dinof lagellate),Pyramimona(chlorophyte), andChrysochromulina(haptophyte)(Fig.6; Supplementary Table S10). Several previous studies generally considered that there are allelopathic interactions between macroalgae and phytoplankton.The macroalgae have adverse eff ects on diatoms;however, their eff ects on the dinof lagellates are controversial: some experiments found that macroalgae inhibit the growth of dinof lagellates,and the others found that macroalgae promote the development of dinof lagellates (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,Skeletonemawas the dominant genus with a relative abundance of more than 20% in the NPEP assemblages in May,GyrodiniumandKarlodiniumbecame the dominant genera with high relative abundance (7.8% and 37.5%) in August, andGonyaulax(dinof lagellate)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 dinof lagellates.

Brief ly, 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 identif ied and annotated, the NPEP species diversity was massive,considerably more signif icant than previously viewed in this area. The composition of the NPEP communities exhibited a seasonal succession from diatom-dominated to dinof lagellate-dominated from May to September, andSkeletonema,Karlodinium,andGonyaulaxwere 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 diatomdominated to dinof lagellate-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 f ield samples and Oceanographic Data Center, IOCAS, for molecular data analysis.