Diversity and seasonal variation of marine phytoplankton in Jiaozhou Bay, China revealed by morphological observation and metabarcoding*

Tiantian CHEN , Yingxin ZHANG , Shuqun SONG , Yun LIU ,Xiaoxia SUN , Caiwen LI ,4 ,

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

2 Laboratory of Marine Ecology and Environmental Science, Qingdao National Laboratory for Marine Science and Technology,Qingdao 266200, China

3 Center for Ocean Mega-Science, Chinese Academy of Sciences, Qingdao 266071, China

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

5 Jiaozhou Bay National Marine Ecosystem Research Station, Chinese Academy of Sciences, Qingdao 266071, China

Abstract Phytoplankton are central components of marine environments, and are major players in the production and respiration budgeting. However, their diversity and distribution patterns are still poorly understood due largely to their small sizes and inconspicuous morphology that have been determined via the application of traditional morphology methods over the past two decades. To better understand the composition and diversity of phytoplankton in Jiaozhou Bay, China, seasonal sampling was carried out in 2019 and samples were analyzed with morphological observations and high-throughput sequencing,from which obvious seasonal variations in phytoplankton composition and proportional abundances were uncovered. Metabarcoding revealed far more diversity and species richness of phytoplankton than morphological observations, especially with respect to dinof lagellates. Diatoms were the most dominant phytoplankton group throughout the year, of which Thalassionema and Skeletonema were co-dominant in the bay. Parasitic dinof lagellates (e.g. Amoebophrya), which is often overlooked in the morphological observations, were in dominance and high diversity in the metabarcoding dataset, thus more attention should be paid to exploring the potential role of parasitic dinof lagellates. Temperature, chlorophyll a, and nutrient levels were the main inf luential factors on the distribution of phytoplankton. This study provided a comprehensive morphological and molecular description of phytoplankton and clearly demonstrated the importance of molecular technology in exploring phytoplankton communities. More-widespread use of molecular technology will facilitate deeper understanding of the ecological importance of the diff erent species.

Keyword: phytoplankton; high-throughput sequencing; diversity; morphological observation; Jiaozhou Bay

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 diffi cult 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 km2and 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 f ield 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 eff ects of climate and environmental changes.

Fig.1 Location of the sampling stations (black dots) in the Jiaozhou Bay

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 f ixed stations. Surface water was f irstly f iltered through a 200-μm mesh sieve to remove most of the mesozooplankton and large particles. The chlorophylla(chla) was collected by f iltering 250 mL of the prescreened water samples through the Whatman GF/F f ilters (0.7-μm pore size, 25 mm in diameter, Whatman Inc.) with a vacuum pump(＜0.04 MPa). Then the f ilters were wrapped against light in aluminum foils, and stored at -20 °C for further processing. Chl-aconcentrations were measured in the laboratory using a f luorometric method after 24 h in 90% acetone in the dark at 4 °C(Parsons et al., 1984) with a Turner Designs 10-005R f luorometer. Water samples in triplicates were f iltered through 0.45-μm pore-size cellulose f ilters for analysis of inorganic nutrients (nitrate/NO3ˉ, nitrite/NO2ˉ, ammonium/NH4+, phosphate/PO43ˉ, and silicate/SiO44ˉ), which were analyzed by a continuous f low analyzer (SANplus system, SKALAR Inc.) by colorimetric methods (Parsons et al., 1984).

Aliquots of 250 mL were preserved with buff ered formaldehyde (2% f inal concentration). Species identif ication and cell enumeration were performed with the Utermöhl’s sedimentation assay (1958).Brief ly, 10-25 mL of the subsamples were settled in Hydro-bios chambers for 24 h. Samples were examined at 200× or 400× magnif ication 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 pref iltered 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 f ilters 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 PowerSoilTMDNA Isolation Kit (MO BIO Laboratories, USA) according to the manufacturer’s protocol. The V4 region of the 18S rRNA gene was amplif ied 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 conf irmed through visualization on a 2% agarose gel electrophoresis, and then purif ied with GeneJETTMGel 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 f irst trimmed to remove the adapters and primers using the Quantitative Insights Into Microbial Ecology (QIIME v. 1.9.1) and f iltered according to previously described methods (Caporaso et al., 2010). Brief ly, all low-quality or unassembled sequences, including chimeric sequences, ambiguous bases, short sequences (＜200 bp), nonspecif ic amplif ication, 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 affi liated withAmoebophryawere f iltered for further analysis. TheAmoebophryasequences 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

At-test was performed to determine the diff erence 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 signif icance level ofP＜0.05.

Fig.2 Variations of major environmental parameters in the surface layer of the Jiaozhou Bay

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

whereniis the cell abundance ofispecies,fiis the frequency of occurrence forispecies, andNis 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:

wherePirepresented the proportion ofispecies densities in the total dinof lagellate 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 quantif ied by the marginal and sequential tests.dbRDA was performed in the statistics program R using the “vegan” package. Spearman correlation analysis betweenAmoebophryaand dinof lagellate 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 signif icant diff erence 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 °C (March) and a maximum of 26.56 °C(September) (Fig.2a). Salinity ranged from 31.02 to 32.07 in Jiaozhou Bay (Fig.2b). Chlorophyllaconcentrations 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, chlawas 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 signif icant diff erence among the three designed areas in four seasons (P＞0.05, Fig.2d-f).

Fig.3 Variation of phytoplankton abundance in diff erent seasons in the Jiaozhou Bay

3.2 Microscopic examination of marine phytoplankton

Four groups of phytoplankton were identif ied 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 diff erent.Chaetoceroswas abundant throughout the year. In March,ProrocentrumandGymnodiniumhad a higher relative abundance.Thalassionemawas abundant in March and December, especially in the inner and mouth areas. Higher abundances ofSkeletonemaandCerataulinawere observed in March and September.

3.3 Metabarcoding of phytoplankton community

Fig.4 The composition of phytoplankton communities at class level (a) and genus level (b) in the Jiaozhou Bay determined by morphological method

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 f iltering 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 dinof lagellate (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 wereThalassiosira,Skeletonema,Ostreococcus,Amoebophrya, andGyrodiniellum.Thalassiosirahad 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.Skeletonemadominated 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 representingOstreococcuswas also more abundant in the inner and mouth area of the bay than in the outer bay area. Whereas the proportion ofGyrodiniellumsequences was higher in the mouth area and the outer bay area in March, they were relatively lower inside the bay in September.Specif ically, there were high proportions of the parasitic dinof lagellateAmoebophryaidentif ied in the bay throughout the year, which was neglected in the morphological observations and further described in the following section.

Fig.5 Venn diagrams showing number of unique and shared OTUs for the three zones of the Jiaozhou Bay (a), and the composition of phytoplankton communities at class level (b) and genus level (c) determined by metabarcoding approach

Fig.6 Phytoplankton diversity compared between morphological and metabarcoding methods

3.4 Comparison of the phytoplankton community identif ied with the two methods

A signif icantly 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 identif ication (Fig.6b).

Fig.7 Abundance (a) and occurrence (b) of the Amoebophrya taxa detected in Jiaozhou Bay

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 identif ied 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 identif ied as affi liated with the genera ofAmoebophrya. Among those OTUs, 33 OTUs were assigned to 9Amoebophryastrains, includingAmoebophryaspp. infectingAlexandriumaffi ne,Gonyaulaxpolygramma,Prorocentrumminimum,Akashiwosanguinea,Scrippsiellasp.,Karlodinium micrum,Prorocentrummicans,Cochlodinium polykrikoides, andDinophysissp. (Fig.7a). The assignedAmoebophryaspecies possessed unique patterns of distribution in the sampling area, while the distribution of unidentif iedAmoebophryaspecies was relatively uniform (Fig.7b). TheAmoebophryasp. exC.polykrikoidesdominated the Amoebophryidae community, representing approximately 14.38% ofAmoebophryasequences.WhileC.polykrikoideshad 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.Amoebophryasp. exP.micanswas the second dominant species (4.84%), followed byAmoebophryasp. exScrippsiellasp. (2.49%) andG.polygramma(1.03%). Other species constituted less than 1% ofAmoebophryasequences, respectively;however, there were approximately 76.04% ofAmoebophryasequences (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 ofAmoebophryasequences 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 ofAmoebophryaOTUs, including sequences fromAmoebophryaspp.infectingDinophysissp.,A.sanguinea,G.polygramma,A.affi ne,P.micans, andC.polykrikoides. Clade II containedAmoebophryastrains infectingDinophysissp.,P.minimum, andK.micrum. Clade III includedAmoebophryasp. infectingDinophysissp.,C.polykrikoides, andScrippsiellasp.No particular putative hosts was assigned to clade IV.Clade V includedAmoebophryastains infectingC.polykrikoides. Further, the phylogeny ofAmoebophryastrains did not appear to be associated with the phylogeny of their host dinof lagellates,further manifested by the fact thatAmoebophryawas a species complex and implied the high genetic divergence of the parasitoid. For example, theAmeobophryastrains infectingP.micansandC.polykrikoideswere more strongly related theAmoebophryasp. exA.affi neandDinophysissp.,respectively.

Fig.8 Phylogenetic tree based on maximum-likelihood (ML) analysis of Amoebophrya sequences obtained from environmental water samples of the Jiaozhou Bay via high-throughput sequencing

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 inf luential factor aff ecting the composition of phytoplankton, andGymnodiniumandAmoebophryawere more positively related to temperature in the morphological and metabarcoding data, respectively.Nutrients (nitrate/NO3ˉ, nitrite/NO2ˉ, and ammonium/NH4+) were also determing factors. For example,Prorocentrumwas positively related to NH4+, whileBacillariawas negatively related to NH4+according to the morphological method (Fig.9a).Neoceratiumwas positively related to NH4+, whileThalassiosirawas negatively related to NH4+;Noctilucawas positively correlated to NO3ˉ, andOstreococcuswas negatively related to NO2ˉ in the metabarcoding data(Fig.9b).

Fig.9 Distance-based redundancy analysis (dbRDA) ordination plot of phytoplankton groups in morphological (a) and metabarcoding datasets (b) versus environmental variables, and the Spearman’s correlation analysis between Amoebophrya and top 30 dinof lagellate genera (c)

Further variation partitioning analysis showed that physical factor (temperature) and chemical factors(NO3ˉ, NO2ˉ, NH4+) 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.

Fig.10 The Spearman’s correlation analysis between Amoebophrya clades and dinof lagellate hosts

Amoebophryais an obligate endoparasitoid infecting wide ranges of marine dinof lagellates, 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 ofAmoebophryawere negatively related to the relative abundances of most genera of dinof lagellates,particularlyHaplozoon,Prorocentrum,Lepidodinium,andDinophysis(Fig.9c). Moreover,Amoebophryacorrelated negatively with their dinof lagellate hosts(Fig.10), implying that the parasitoid may be a factor in shaping the local dinof lagellate assemblage, and the ecological roles of this group should not be neglected.

“I think it is gone,” he said. But it was not gone; it was one of those bits of the looking-glass—that magic mirror, of which we have spoken—the ugly glass which made everything great and good appear small and ugly, while all that was wicked and bad became more visible, and every little fault could be plainly seen. Poor little Kay had also received a small grain in his heart, which very quickly turned to a lump of ice. He felt no more pain, but the glass was there still. “Why do you cry?” said he at last; “it makes you look ugly. There is nothing the matter with me now. Oh, see!” he cried suddenly, “that rose is worm-eaten, and this one is quite crooked15. After all they are ugly roses, just like the box in which they stand,” and then he kicked the boxes with his foot, and pulled off the two roses.

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 dinof lagellates (e.g.AmoebophryaandDuboscquella).Morphological observation provided cell abundance of dominant groups and reliable assessment on the abundance and community structure of phytoplankton,which can be a benef icial supplement to the molecular sequencing methods.

The two methods diff ered 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 identif ied. 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 diff erences 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 aff ected 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 f iltering 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 quantif ication of the rare species and pico-phytoplankton, particularly the dinof lagellates in the class of Syndiniophyceae. Most species in the group were parasites or symbionts,which were diffi cult 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 dinof lagellates were signif icantly 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 dinof lagellates (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 affi liated withAmoebophrya,Duboscquella, and other parasitic dinof lagellates (e.g.Hematodinium).Amoebophryawas 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). TheTaraOceansstudy further verif ied thatAmoebophryawas one of the dominant groups in oceans worldwide (de Vargas et al., 2015; Lima-Mendez et al., 2015). The host ranges ofAmoebophryaare extremely diverse,including dinof lagellates, radiolarians, ciliates,cercozoans, and even f ish 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,Amoebophryashowed positive correlations with temperature and nutrient conditions,and presented the highest relative abundance in the outer area of Jiaozhou Bay. Temperature and nutrient conditions inf luenced the success ofAmoebophryaby altering parasitoid reproductive output and infectivity of progeny (Yih and Coats, 2000). Negative correlations were also found betweenAmoebophryaand multiple groups of dinof lagellate, especially signif icant forProrocentrumandDinophysis(Fig.9c).In theTaraOceansstudy,Amoebophryawas also enriched in negative associations with Dinophyceae(Lima-Mendez et al., 2015). The fact thatAmoebophryalives 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 dinof lagellate 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 (NO3ˉ, NO2ˉ, and NH4+) were found to aff ect 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 diff erent 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.ThalassiosiraandSkeletonemawere the dominant species and preferred lower temperatures. Generally, temperature and nutrients are considered to be the signif icant 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 signif icant correlation of phytoplankton abundance with those nutrients.

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 inf luenced by temperature, chla, 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 identif ication. Morphological observation provided reliable assessment on the abundance and community structure of phytoplankton,which can be a benef icial 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/VChuangxinfor logistic support during cruises.Comments and suggestions from two anonymous reviewers are gratefully appreciated.