1 INTRODUCTION

Phytoplankton are widely distributed in the marine and estuarine environments and play key roles in supporting the productivity of the ecosystem. These microscopic organisms are one of the initial biological components from which energy is transferred to higher organisms through the food chain (Raymont, 1983). The structure and succession of phytoplankton communities are very sensitive to environmental variation and are regulated by a combination of physical, chemical, and biological factors (Raymont, 1983). Thereby, variation in the populations can act as important indicators and evaluation indices for environmental assessment of the aquatic ecosystems.

Phytoplankton are traditionally described based on morphological characteristics, making it difficult to distinguish the cryptic species complexes, species with multiple life stages, or species with distinct sexual dimorphism (Utermöhl, 1958; de Vargas et al., 2015). As a result, pico- and nano-sized phytoplankton have been commonly overlooked in earlier studies due to their small sizes and inconspicuous morphology, although they are very important contributors to the primary production (López-García et al., 2001). Recently, molecular approaches have been widely applied to study the diversity and abundance of phytoplankton in various environments, and treated as vital approaches to the taxonomic revision of phytoplankton (de Vargas et al., 2015; Lima-Mendez et al., 2015). The 18S rDNA gene, containing both conserved and hypervariable regions, is universally present in living eukaryotes and has become the most widely used biomarker to detect and classify known species present in marine eukaryotic communities (Ki, 2012; Morard et al., 2018). In particular, the 18S rDNA metabarcoding can detect rare species (López-García et al., 2001), cryptic species (Chen et al., 2019b), and pico-phytoplankton (Moon-van der Staay et al., 2001), and has been successfully applied to explore the diversity and community structure of marine eukaryotic organisms (de Vargas et al., 2015; Lima-Mendez et al., 2015).

Jiaozhou Bay, located in the southwest part of the Yellow Sea, is a semi-enclosed bay with an area of approximately 367 km² and average depth of 7 m. Since the 1990s, coastal bay environment has changed greatly due to the increased nutrient loadings along with the rapidly expanding economic activities and human population, which leads to subsequent variations in the species composition and community structure of marine phytoplankton (Shen, 2001; Liu et al., 2011). Long-term investigations on phytoplankton communities in Jiaozhou Bay had accumulated valuable historical data (Guo et al., 2019). While most studies focusing on the species composition and biodiversity of marine phytoplankton were carried out using traditional microscopic methods, it is challenging to identify small cryptic species. In recent years, many new and revised species, genera, and higher taxonomic groups have been proposed with assistance of molecular evidences (Mahon et al., 2011; Liu et al., 2020). The physico-chemical parameters vary dramatically among seasons in the Jiaozhou Bay, and the detailed seasonal variation in the phytoplankton community and the determining factors remain unclear. Herein, a comprehensive field investigation coupled with microscopic observation and high-throughput sequencing were conducted to reveal the spatial-temporal distribution and composition of marine phytoplankton in the bay, and the correlation between phytoplankton community and major environmental parameters were further assessed to highlight the effects of climate and environmental changes.

2 MATERIAL AND METHOD 2.1 Study area and cruises

Jiaozhou Bay is a landlocked bay, connecting to the Yellow Sea through a narrow opening. Four cruises covering three zones (inner, mouth, and outer bay) were conducted in March (spring), June (summer), September (autumn), and December (winter) 2019 (Fig. 1). Basic physical parameters (temperature and salinity) of surface water were measured in situ with a conductivity-temperaturedepth (CTD) sensor, and the data were provided by the Jiaozhou Bay National Marine Ecosystem Research Station, Institute of Oceanology, Chinese Academy of Sciences.

2.2 Sampling and processing of environmental water samples

Four cruises were carried out in the Jiaozhou Bay, China in 2019, and each cruise was completed in two or three days (March 11, 13, and 14; June 11–13; September 9, 10, and 12; and December 12–13). Surface water samples (< 3-m depth) were collected with 12-L Niskin bottles (Wildlife Supply Company, MI, USA) from 12 fixed stations. Surface water was firstly filtered through a 200-μm mesh sieve to remove most of the mesozooplankton and large particles. The chlorophyll a (chl a) was collected by filtering 250 mL of the prescreened water samples through the Whatman GF/F filters (0.7-μm pore size, 25 mm in diameter, Whatman Inc.) with a vacuum pump (< 0.04 MPa). Then the filters were wrapped against light in aluminum foils, and stored at -20 ℃ for further processing. Chl-a concentrations were measured in the laboratory using a fluorometric method after 24 h in 90% acetone in the dark at 4 ℃ (Parsons et al., 1984) with a Turner Designs 10-005R fluorometer. Water samples in triplicates were filtered through 0.45-μm pore-size cellulose filters for analysis of inorganic nutrients (nitrate/NO₃^-, nitrite/NO₂^-, ammonium/NH₄⁺, phosphate/PO₄^3-, and silicate/SiO₄^4-), which were analyzed by a continuous flow analyzer (SANplus system, SKALAR Inc.) by colorimetric methods (Parsons et al., 1984).

Aliquots of 250 mL were preserved with buffered formaldehyde (2% final concentration). Species identification and cell enumeration were performed with the Utermöhl's sedimentation assay (1958). Briefly, 10–25 mL of the subsamples were settled in Hydro-bios chambers for 24 h. Samples were examined at 200× or 400× magnification using an inverted light microscope (IX71, Olympus, Japan). The lower size limit of resolution for analysis was ~5 μm. The phytoplankton community was characterized by identifying taxa to genus or species level.

2.3 Sampling, processing, and sequencing of environmental DNA samples

Two liters of prescreened waters were prefiltered through 3-μm pore-size Nuclepore membranes (47-mm diameter, Whatman, Piscataway, NJ, USA) using a vacuum pump at low pressure (< 0.03 MPa). The filters were then transferred into a 2.0-mL DNAfree cryogenic vials (Axygen, USA) and frozen immediately in liquid nitrogen until being further processed.

The genomic DNA of water samples was performed using the PowerSoil^TM DNA Isolation Kit (MO BIO Laboratories, USA) according to the manufacturer's protocol. The V4 region of the 18S rRNA gene was amplified using the primers 528F (5′-GCGGTAATTCCAGCTCCAA-3′) and 706R (5′-AATCCRAGAATTTCACCTCT-3′) (Elwood et al., 1985). PCR reactions were performed as described in Cheung et al. (2010). PCR products were confirmed through visualization on a 2% agarose gel electrophoresis, and then purified with GeneJET^TM Gel Extraction Kit (Qiagen, North Rhine-Westphalia, Germany). 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 was deposited in the NCBI Sequence Read Archive (SRA) under accession number PRJNA658427.

2.4 Bioinformatic of metabarcoding amplicons

Raw sequences were first trimmed to remove the adapters and primers using the Quantitative Insights Into Microbial Ecology (QIIME v. 1.9.1) and filtered according to previously described methods (Caporaso et al., 2010). Briefly, all low-quality or unassembled sequences, including chimeric sequences, ambiguous bases, short sequences (< 200 bp), nonspecific amplification, homologous regions, were eliminated for further analysis (Cheung et al., 2010). Remaining reads were then strictly de-replicated and aligned to the Silva reference alignment. The resulting sequences were clustered into operational taxonomic unit (OTU) based on 97% similarity (Edgar, 2010). The representative sequence of each OTU (most abundant) was assigned using BLAST (Altschul et al., 1990). All OTU assigned to metazoans and fungi were removed from the data set.

Sequences affiliated with Amoebophrya were filtered for further analysis. The Amoebophrya sequences obtained in the present study were previously submitted to the National Center for Biotechnology (NCBI) GenBank under accession nos. MT269074–MT269265. Multiple alignments were performed by Clustal X ver. 1.8 (Van de Peer et al., 2000), then ambiguous nucleotides and gaps were removed before phylogenetic analysis. Maximumlikelihood (ML) analyses were conducted using the MEGA 6.06 with General Time Reversible (GTR) model. One thousand bootstrap replicates were performed to test the robustness of each node in the phylogenetic tree.

2.5 Statistical analysis

A t-test was performed to determine the difference in physical parameters between the three zones in Jiaozhou Bay. All statistical tests were performed with the SPSS 22.0 software (SPSS, Chicago, IL, USA), with the significance level of P < 0.05.

The dominant species of phytoplankton in morphological datasets were determined by the McNaughton index (Y), according to the following formula:

where n_i is the cell abundance of i species, f_i is the frequency of occurrence for i species, and N is the total number for all species.

The diversity of phytoplankton in morphological datasets was determined by the Shannon-Wiener diversity index (H′) (Shannon and Weaver, 1949), according to the following formula:

where P_i represented the proportion of i species densities in the total dinoflagellate density in the sample.

The morphological and metabarcoding datasets of phytoplankton abundance were transformed into relative proportions prior to conducting multivariate analysis. Distance-based redundancy analysis (dbRDA) was used to correlate the phytoplankton community with the environmental variables. The contribution of environmental variables was quantified by the marginal and sequential tests. dbRDA was performed in the statistics program R using the "vegan" package. Spearman correlation analysis between Amoebophrya and dinoflagellate genera was performed using the R package of "corrplot" and "ggcorrplot", respectively. Prior to these analyses, relative abundances of each genus and the environmental parameters were log (x+1) transformed and normalized.

3 RESULT 3.1 Environmental characteristics

The surface water temperature and salinity were relatively consistent throughout Jiaozhou Bay, and exhibited no significant difference among the three designated areas (the inner bay, mouth bay, and outer bay) in four seasons (P > 0.05). The variation of temperature followed the classical seasonal dynamics of northern China seas, characterized by a minimum of 6.58 ℃ (March) and a maximum of 26.56 ℃ (September) (Fig. 2a). Salinity ranged from 31.02 to 32.07 in Jiaozhou Bay (Fig. 2b). Chlorophyll a concentrations showed an annual cycle characterized by the highest value in June and lowest value in December at all sampling sites, ranging from 0.92 to 5.80 μg/L. In general, chl a was relatively higher in the inner bay, followed by the mouth area and the outer bay (Fig. 2c), that was in accordance with the nutrient concentrations of the seawaters from the inner bay to the outer bay areas (Fig. 2d–f). Highest values of DIN (nitrate, nitrite, and ammonium), dissolved reactive phosphate (DRP), and silicate were found in September in the inner bay, while there was no significant difference among the three designed areas in four seasons (P > 0.05, Fig. 2d–f).

3.2 Microscopic examination of marine phytoplankton

A total of 112 phytoplankton species under 60 genera were identified in the Jiaozhou Bay during the sampling period. Most of the ecotypes of phytoplankton were temperate coastal species, with total cell abundance ranging from 41.44×10³ to 310.48×10³ cells/L in March, 105.88×10³ to 864.12×10³ cells/L in June, 42.20×10³ to 1 001.24×10³ cells/L in September and 23.03×10³ to 55.42×10³ cells/L in December (Fig. 3a). The phytoplankton community was mainly composed of diatoms and dinoflagellates, dominated by Skeletonema sp., Thalassiosira sp., Skeletonema costatum, Chaetoceros curvisetus, Nitzschia sp., Thalassiosira rotula, and Bacillaria paxillifera for the diatoms and Prorocentrum minimum for the dinoflagellate (Fig. 3a). The distribution of phytoplankton abundances showed dominance of diatoms in four seasons, with total cell abundance ranging from 33.17×10³ to 252.07×10³ cells/L in March, 38.75×10³ to 795.91×10³ cells/L in June, 23.05×10³ to 979.08×10³ cells/L in September and 19.07×10³ to 54.55×10³ cells/L in December, and the abundances were always high in inner bay and mouth area of the bay (Fig. 3b). The cell abundance of dinoflagellate was low in the bay (Fig. 3c), and the distribution of the dinoflagellate was quite different from that of the total phytoplankton and diatoms.

Four groups of phytoplankton were identified at the class level, namely Bacillariophyceae, Mediophyceae, Dinophyceae, and Dictyochophyceae. Mediophyceae always dominated the phytoplankton communities in March and September, and higher abundances of Dinophyceae was observed in June, especially in the inner and outer bay seawaters. Bacillariophyceae was more abundant in December, and was accompanied by Mediophyceae and Dinophyceae (Fig. 4a). For the 35 dominant genera, a detailed abundance pattern across all samples was illustrated by a heatmap (Fig. 4b). The dominant phytoplankton genera in four seasons were quite different. Chaetoceros was abundant throughout the year. In March, Prorocentrum and Gymnodinium had a higher relative abundance. Thalassionema was abundant in March and December, especially in the inner and mouth areas. Higher abundances of Skeletonema and Cerataulina were observed in March and September.

3.3 Metabarcoding of phytoplankton community

A total of 6 173 414 raw rDNA sequences were obtained from high-throughput sequencing of the 48 environmental samples collected in Jiaozhou Bay and 5 929 352 sequences were retained after quality filtering and were further clustered into 4 689 OTUs (290 to 852 per sample), of which 2 405 OTUs (51.29%) were common in all sampling sites (Fig. 5a). The community was dominated by diatoms (Mediophyceae and Mamiellophyceae) and dinoflagellate (Dinophyceae) (Fig. 5b). Mediophyceae contributed substantially to the diatom sequences in March, June, and September (comprising over 51.02%, 36.68%, and 52.10% of the sequences, respectively), and co-dominated with Dinophyceae in September and December (51.31% and 38.16% of the sequences, respectively).

The communities of phytoplankton in the mouth area and outer bay area shared 3 122 OTUs, which were mainly dominated by Dinophyceae and Mediophyceae, contributing to approximately 21.00%–61.10% and 12.31%–48.24% of total phytoplankton sequences, respectively (Fig. 5b), while the contributions of sequences from other classes varied between seasons. Particularly in March, Dinophyceae and Mediophyceae were accompanied by Syndiniophyceae, Mamiellophyceae, Cryphophyceae, Dictyochophyceae, and Bacillariophyceae. In June and December, there was a higher contribution of sequences from Mamiellophyceae and Cryptophyceae besides Dinophyceae and Mediophyceae. In September, higher abundances of Chlorophyceae and Syndiniophyceae were observed in the inner and outer bay waters, respectively. The active phytoplankton communities in the inner bay area shared 2 737 and 2 615 OTUs with those in the mouth and outer bay area, respectively (Fig. 5a). The sequences of Mamiellophyceae were more abundant in the inner bay in the four seasons (22.31%, 13.33%, 23.04%, and 12.18%). The sequence contribution of Dinophyceae and Syndiniophyceae was increased in June and September, while there was higher sequence contribution of Mediophyceae in March and December, and a higher contribution of sequences from Chlorophyceae in September. It is noteworthy that members of the Chlorodendrophyceae were detected only inside of the Jiaozhou Bay in September (Fig. 5b).

The seasonal variation in the relative abundance of 35 dominant genera was further illustrated by a heatmap (Fig. 5c). The most prevalent genera throughout the sampling period were Thalassiosira, Skeletonema, Ostreococcus, Amoebophrya, and Gyrodiniellum. Thalassiosira had a higher relative abundance on a per sample basis, especially the inner bay, but with a lower abundance in December when it constituted less than 10% of total phytoplankton sequences. Skeletonema dominated the phytoplankton sequences in the inner and mouth area of the bay in March and December, while lower abundance was observed in September (less than 1%). Sequences representing Ostreococcus was also more abundant in the inner and mouth area of the bay than in the outer bay area. Whereas the proportion of Gyrodiniellum sequences was higher in the mouth area and the outer bay area in March, they were relatively lower inside the bay in September. Specifically, there were high proportions of the parasitic dinoflagellate Amoebophrya identified in the bay throughout the year, which was neglected in the morphological observations and further described in the following section.

3.4 Comparison of the phytoplankton community identified with the two methods

A significantly higher number of phytoplankton OTUs (360–718) was found compared to the number of morphospecies (13–49) (Fig. 6a). The community composition determined by metabarcoding exhibited higher Shannon-Wiener index values than that obtained through morphological identification (Fig. 6b). The phytoplankton communities characterized by morphological and metabarcoding methods were compared at the generic level due to morphological and molecular complexity. Metabarcoding detected four times more phytoplankton OTUs than morphospecies in our 48 samples considering that each OTU represented one species at 97% sequence similarity. Forty-seven genera were detected by the two methods, while 13 and 193 genera were detected solely in the morphological observations and metabarcoding, respectively. The genera identified by morphological methods in Dictyochophyceae and Mediophyceae were all observed in the metabarcoding dataset, while 13 genera from other two phytoplankton classes were not fully detected by metabarcoding.

3.5 The high abundance and diversity of Amoebophrya in Jiaozhou Bay

A total of 17 724 sequence reads, assigned to 192 OTUs, were identified as affiliated with the genera of Amoebophrya. Among those OTUs, 33 OTUs were assigned to 9 Amoebophrya strains, including Amoebophrya spp. infecting Alexandrium affine, Gonyaulax polygramma, Prorocentrum minimum, Akashiwo sanguinea, Scrippsiella sp., Karlodinium micrum, Prorocentrum micans, Cochlodinium polykrikoides, and Dinophysis sp. (Fig. 7a). The assigned Amoebophrya species possessed unique patterns of distribution in the sampling area, while the distribution of unidentified Amoebophrya species was relatively uniform (Fig. 7b). The Amoebophrya sp. ex C. polykrikoides dominated the Amoebophryidae community, representing approximately 14.38% of Amoebophrya sequences. While C. polykrikoides had a lower relative abundance on a per sample basis, but with a higher abundance in December when it constituted approximately 10% of the total phytoplankton sequences. Amoebophrya sp. ex P. micans was the second dominant species (4.84%), followed by Amoebophrya sp. ex Scrippsiella sp. (2.49%) and G. polygramma (1.03%). Other species constituted less than 1% of Amoebophrya sequences, respectively; however, there were approximately 76.04% of Amoebophrya sequences (159 OTUs) not assigned to particular putative hosts, demonstrating the current limited understanding on their high abundance and diversity in marine environments.

Five distinct clades of members of the Amoebophryidae organisms were categorized when the genetic distances of Amoebophrya sequences were evaluated with ML analyses (Fig. 8). The three most common clades (clade I, clade II, and clade III) accounted for 93.75% of the total number of OTUs. Clade I contained the largest number of Amoebophrya OTUs, including sequences from Amoebophrya spp. infecting Dinophysis sp., A. sanguinea, G. polygramma, A. affine, P. micans, and C. polykrikoides. Clade II contained Amoebophrya strains infecting Dinophysis sp., P. minimum, and K. micrum. Clade III included Amoebophrya sp. infecting Dinophysis sp., C. polykrikoides, and Scrippsiella sp. No particular putative hosts was assigned to clade IV. Clade V included Amoebophrya stains infecting C. polykrikoides. Further, the phylogeny of Amoebophrya strains did not appear to be associated with the phylogeny of their host dinoflagellates, further manifested by the fact that Amoebophrya was a species complex and implied the high genetic divergence of the parasitoid. For example, the Ameobophrya strains infecting P. micans and C. polykrikoides were more strongly related the Amoebophrya sp. ex A. affine and Dinophysis sp., respectively.

3.6 Associations between phytoplankton communities and environmental variables

As inferred from the dbRDA ordination plot, phytoplankton communities inferred from morphological observation and metabarcoding methods were clearly separated by seasons (Fig. 9a–b). The two RDA dimensions collectively explained 23.53% and 33.58% of the variation in the phytoplankton composition according to the morphological and metabarcoding methods, respectively. Temperature seemed to be the most influential factor affecting the composition of phytoplankton, and Gymnodinium and Amoebophrya were more positively related to temperature in the morphological and metabarcoding data, respectively. Nutrients (nitrate/NO₃^-, nitrite/NO₂^-, and ammonium/ NH₄⁺) were also determing factors. For example, Prorocentrum was positively related to NH₄⁺, while Bacillaria was negatively related to NH₄⁺ according to the morphological method (Fig. 9a). Neoceratium was positively related to NH₄⁺, while Thalassiosira was negatively related to NH₄⁺; Noctiluca was positively correlated to NO₃^-, and Ostreococcus was negatively related to NO₂^- in the metabarcoding data (Fig. 9b).

Further variation partitioning analysis showed that physical factor (temperature) and chemical factors (NO₃^-, NO₂^-, NH₄⁺) accounted for 8.36% and 17.88% of phytoplankton community variations in the morphological data, and 27.86% and 39.16% in the OTU data, respectively (inserted picture, Fig. 9a–b). In addition, the interaction between physical and chemical parameters could explain 6.29% and 8.67% of the variation from the two methods, respectively. Whereas, 73.76% and 32.98% of the phytoplankton community variation from the two methods could not be explained by these components, indicating that other biotic and/or abiotic parameters might have impacts in this process.

Amoebophrya is an obligate endoparasitoid infecting wide ranges of marine dinoflagellates, it was also the fourth dominant genus in the metabarcoding dataset, but was commonly overlooked in the morphological observations because of its parasitic characteristics. Additionally, the relative abundances of Amoebophrya were negatively related to the relative abundances of most genera of dinoflagellates, particularly Haplozoon, Prorocentrum, Lepidodinium, and Dinophysis (Fig. 9c). Moreover, Amoebophrya correlated negatively with their dinoflagellate hosts (Fig. 10), implying that the parasitoid may be a factor in shaping the local dinoflagellate assemblage, and the ecological roles of this group should not be neglected.

4 DISCUSSION

Phytoplankton communities in the Jiaozhou Bay have been well studied in light microscopy for many years (Shen, 2001; Guo et al., 2019). However, rare and cryptic species could not be easily recognized in light microscopy, and their molecular studies are still rare. In the present study, the community structure and diversity of marine phytoplankton in the Jiaozhou Bay was investigated based on microscopy and metabarcoding methods. Metabarcoding revealed much more diversity and species richness of phytoplankton in comparison to morphological studies, especially with respect to the parasite dinoflagellates (e.g. Amoebophrya and Duboscquella). Morphological observation provided cell abundance of dominant groups and reliable assessment on the abundance and community structure of phytoplankton, which can be a beneficial supplement to the molecular sequencing methods.

The two methods differed greatly in the number of phytoplankton groups detected and the estimated abundance of phytoplankton groups. Twenty two classes of phytoplankton were revealed by metabarcoding, which was much higher than the parallel morphological observation in which only four classes were identified. In the present study, 47 genera were observed by the two methods, 13 and 193 genera were detected solely in the morphological observations and metabarcoding, respectively. These differences may be partially due to the fact that indepth sequencing based on high-throughput sequencing retrieved a higher number of sequences compared with counting according to morphological observations (Stoeck et al., 2009). Secondly, the investigated water volumes of microsopy (10–25 mL) were much less than that of metabarcoding (2 L), which may also result in the less number of phytoplankton groups in microsopy. Thirdly, the metabarcoding results may be affected by the DNA from non-living cells. The stable DNA can remain preserved outside the cell for several years, which may be collected during the filtering process and contributed to the phytoplankton taxon richness and abundance in molecular sequencing (Vlassov et al., 2007; Massana et al., 2015). Fourthly, metabarcoding may reveal concealed phytoplankton diversity when morphologically identical taxa exhibit distinct genetic variation (Huo et al., 2020). Additionally, microscopy is not suitable for diagnosis or quantification of the rare species and pico-phytoplankton, particularly the dinoflagellates in the class of Syndiniophyceae. Most species in the group were parasites or symbionts, which were difficult to identify microscopically because of their small size or lack of morphological characteristics in microscopic observations (Coats, 1999; Jephcott et al., 2016). It is noteworthy that the taxon richness and relative abundance of dinoflagellates were significantly higher in the metabarcoding data than that in the morphological data. The abundance of diatoms were relatively lower in metabarcoding compared to that of the microscopic observation, even higher relative abundance of diatom sequences was observed in sites where the diatom cells were relatively abundant. The bias could be introduced in the metabarcoding processes by DNA extraction (Martin-Laurent et al., 2001), which might result in the absence of diatoms in our metabarcoding process. Besides, extremely high copy numbers of the rDNA gene have been found in dinoflagellates (Gong et al., 2013), and thus might be over-represented in sequencing (Liu et al., 2017).

Members of the Syndiniophyceae are distributed widely with high genetic diversity and abundance in the world oceans (de Vargas et al., 2015; Lima-Mendez et al., 2015). They formed a bulk abundance of phytoplankton as well as various species composition in Jiaozhou Bay in the present study. Environmental sequences belonging to the Syndiniophyceae are mainly affiliated with Amoebophrya, Duboscquella, and other parasitic dinoflagellates (e.g. Hematodinium). Amoebophrya was the fourth dominant genus in the present study, which was observed throughout the survey area and was the most diverse and highly represented in the coastal waters of China (Li et al., 2014; Liu et al., 2017; Chen et al., 2019a, 2020). The Tara Oceans study further verified that Amoebophrya was one of the dominant groups in oceans worldwide (de Vargas et al., 2015; Lima-Mendez et al., 2015). The host ranges of Amoebophrya are extremely diverse, including dinoflagellates, radiolarians, ciliates, cercozoans, and even fish eggs, and one infection can produce hundreds of dinospores (~400) in 2–3 d, which contributes to the relative higher abundance of this group in natural waters (Coats, 1999; Chambouvet et al., 2011).

In the present study, Amoebophrya showed positive correlations with temperature and nutrient conditions, and presented the highest relative abundance in the outer area of Jiaozhou Bay. Temperature and nutrient conditions influenced the success of Amoebophrya by altering parasitoid reproductive output and infectivity of progeny (Yih and Coats, 2000). Negative correlations were also found between Amoebophrya and multiple groups of dinoflagellate, especially significant for Prorocentrum and Dinophysis (Fig. 9c). In the Tara Oceans study, Amoebophrya was also enriched in negative associations with Dinophyceae (Lima-Mendez et al., 2015). The fact that Amoebophrya lives within their host and infection of the parasitoid eventually leads to the mortality of host cells, implied that the parasitoid may be a factor in shaping the local dinoflagellate assemblage, and the ecological roles of this group should not be neglected (Kim et al., 2004; Montagnes et al., 2008). The distribution, diversity, abundance, and composition of phytoplankton are linked with a variety of environmental parameters such as temperature, salinity, and nutrients in oceanic environments (Perumal et al., 2009; Saifullah et al., 2019). Temperature and nutrients (NO₃^-, NO₂^-, and NH₄⁺) were found to affect the composition and structure of the phytoplankton communities inferred from both morphological observation and metabarcoding methods in the present study, which were well recognized in different studies (Liu et al., 2017; Saifullah et al., 2019). The total biomass of phytoplankton was higher in the inner bay in the present study, which were eutrophic waters in Jiaozhou Bay. Thalassiosira and Skeletonema were the dominant species and preferred lower temperatures. Generally, temperature and nutrients are considered to be the significant components playing vital roles in phytoplankton diversity in the estuarine environment (Saifullah et al., 2019). Nitrite, nitrate, and ammonium seemed to be of particular importance in Jiaozhou bay. A similar observation was also reported in Pichavaram mangrove, India (Perumal et al., 2009), and in the Changjiang (Yangtze) River, China (Liu et al., 2017). Silica and phosphate were also reported as controlling factors of phytoplankton communities (Perumal et al., 2009; Liu et al., 2017; Saifullah et al., 2019), though there was no significant correlation of phytoplankton abundance with those nutrients.

Based on the variation partitioning analysis, the selected physical and chemical factors explained 8.36% and 17.88%, 27.86% and 39.16% of the observed variation in the phytoplankton community according to the morphological and metabarcoding methods, respectively (Fig. 9a–b); however, 73.76% and 32.98% of the phytoplankton community variation from the two methods remained unexplained, suggesting that there are other factors, such as oxygen conditions, flow or conductivity, as well as biological interactions, which may co-affect phytoplankton abundance in the natural environments. Further investigations should be performed to confirm the influence of unmeasured parameters and relationships in phytoplankton communities.

5 CONCLUSION

The community structure and diversity of marine phytoplankton in the Jiaozhou Bay, China was studied based on morphological observations and metabarcoding. Phytoplankton communities inferred from the two methods were both influenced by temperature, chl a, and nutrients. Metabarcoding revealed high diversity and species richness, and was shown to be more valuable for assessing picophytoplankton and novel phytoplankton diversity than morphological identification. Morphological observation provided reliable assessment on the abundance and community structure of phytoplankton, which can be a beneficial supplement to the molecular method. Combination of both traditional microscopic observations and molecular sequencing will provide a comprehensive view of the phytoplankton diversity and its correlation with environmental variables in coastal and oceanic waters.

6 DATA AVAILABILITY STATEMENT

Data are available on request from the corresponding author.

7 ACKNOWLEDGMENT

We thank all the crew and captain of the R/V Chuangxin for logistic support during cruises. Comments and suggestions from two anonymous reviewers are gratefully appreciated.