1 INTRODUCTION

Phaeocystis Lagerheim (Haptophyta) is one of the most extensively studied genera of marine phytoplankton worldwide and is recognized both as a nuisance and ecologically important microalga (e.g., Schoemann et al., 2005; Medlin and Zingone, 2007, and references therein). It plays an essential role in the marine ecosystem, particularly the carbon and sulfur biogeochemical cycles by producing substantial amounts of dimethylsulphoniopropionate (DMSP) and volatile catabolite dimethylsulphide (DMS) which can impact the cloud albedo in the atmosphere, hence climate regulation (Stefels et al., 2007; Verity et al., 2007 and references therein). Members of the genus have a complex polymorphic life cycle alternating between solitary free-living cells and gelatinous colonies where hundreds of cells are embedded in a polysaccharide matrix (Rousseau et al., 2007, 2013; Zingone et al., 2011). Colonial species can form nearly monospecific extensive blooms (contributing > 90% of total phytoplankton abundance), which are associated with dramatic ecosystem changes and negative impacts on fisheries and fish farming from the polar to tropical oceans (Hansen et al., 2003; Peng et al., 2005; Schoemann et al., 2005; Smith et al., 2014; Xu et al., 2019). In addition to the well-documented colony-forming species, the non-colonial species of Phaeocystis species have been suggested as potentially important components of open ocean nanoplankton communities (Zingone et al., 1999; Lin et al., 2014; Gran-Stadniczeñko et al., 2017).

Up to now, ten Phaeocystis species have been described around the world since the establishment of the genus Phaeocystis in 1893 (Medlin and Zingone, 2007; Andersen et al., 2015; Shen and Qi, 2021). Among them, three colonial species, P. globosa Scherffel, P. pouchetii (Hariot) Lagerheim, and P. antarctica Karsten are the most frequently recorded species and are known to form dense blooms in major climatic regions, i.e., the Arctic, the Antarctic, and temperate/tropical regions (e.g., Schoemann et al., 2005). Compared with the well-documented colony-forming species, less is known about the occurrence and significance of solitary species, including P. scrobiculata Moestrup, P. cordata Zingone et Chrétiennot-Dinet, P. jahnii Zingone, which are supposed to be distributed worldwide (Moestrup, 1979; Zingone et al., 1999). The newly discovered unicellular species P. rex Andersen, Bailey, Decelle & Probert was observed from the Arabian Sea for the first time in 2015 (Andersen et al., 2015). Other 3 species, including P. amoeboidea Büttner (1910), P. sphaeroides Büttner (1911), and P. brucei Mangin (1922) have never been found again after their first description (Medlin and Zingone, 2007). Besides the recognized species, several undescribed Phaeocystis species were reported from the South China Sea, and the western Mediterranean Sea based on morphological characteristics (Doan-Nhu and Larsen, 2010; Novarino, 2014). On the other hand, many uncultured clones belonging to this group were detected worldwide based on DNA metabarcoding (de Vargas et al., 2015; Wu et al., 2015; Shih et al., 2019). Most importantly, two unknown Phaeocystis species have been found in the Indian and Pacific Oceans to form endosymbiotic associations with the Acantharia (Decelle et al., 2012), and these symbiont sequences form the basal clades in the phylogeny of Phaeocystis, indicating a wider diversity of the genus need to be characterized.

The first P. globosa bloom in China was recorded in 1997 on the southeastern coast of China (Chen et al., 2002). Since then, P. globosa bloom has spread rapidly and occurred repeatedly along the South China Sea (Shen et al., 2018 and references therein). In 2004, P. globosa bloom broke out in Bohai Bay for the first time (Qu et al., 2008) and bloomed frequently after that in the Bohai Sea and the northern Yellow Sea (Tu et al., 2011; Li et al., 2022). Until 2018, up to 97 Phaeocystis bloom events have been recorded, covering an area of nearly 20 000 km², but only P. globosa was reported and described (Shen et al., 2018; Shen and Qi, 2021; Wang et al., 2021), indicating that the species diversity and spatial distribution of Phaeocystis are poorly understood in China. Moreover, studies have demonstrated that P. globosa strains isolated from coastal China display remarkable genetic diversity, which might include cryptic species or different ecotypes (Song et al., 2021; Zhang et al., 2021; Li et al., 2022). Moreover, bloom of unicellular species has been reported in the East China Sea, where two new record species P. cordata and P. jahnii were found at high density (Lin et al., 2014), indicating the solitary species may have been undersampled, primarily due to its small size and offshore location. Therefore, to clarify the species diversity and geographical distribution of Phaeocystis in coastal China, high-throughput sequencing was used to investigate the species composition and abundance of Phaeocystis in the Bohai Sea, the Yellow Sea, and the East China Sea in the present study.

2 MATERIAL AND METHOD 2.1 Sample collection

Surface water samples were collected from 82 sites in the Bohai Sea (BS), the Yellow Sea (YS), and the East China Sea (ECS) in April 2021 (Fig. 1). The samples of the BS and the YS, and the ECS were collected onboard R/V LanHai 101, and Xiangyanghong 18 using SBE 32 Water Sampler with a SeaBird CTD system (SeaBird, USA), respectively. After collection, 3-L water samples were pre-filtered by a 200-μm mesh sieve to remove large particles and then collected using a vacuum pump onto a 2-μm pore size polycarbonate membrane filter (Millipore, USA). A total of 246 samples were obtained and stored under -80 ℃ before DNA extraction.

At each site, in-situ water temperature and salinity were measured by SeaBird CTD during the field sampling. Water samples for dissolved nutrients (nitrate, phosphate, ammonium, nitrite, and silicate) and chlorophyll-a (Chl-a) determination were filtered through 0.45-μm pore size cellulose filters, then both the water and filters were preserved at -20 ℃ in darkness before further laboratory processing. Nutrient concentrations were determined by SKALAR sa3000/5000 chemistry unit (SKALAR, the Netherlands) for the BS and the YS samples, and SEAL QuAAtro39 nutrient analyzer (SEAL Analytical, Germany) for the ECS samples, respectively. Chlorophyll-a content was measured by TD10-AU fluorometer (Turner Designs, USA) according to the method of Parsons et al. (1984).

2.2 DNA extraction and amplification

Planktonic community genomic DNA was extracted using the FastDNA^®Spin Kit for Soil (MP Biomedicals, USA) according to the manufacturer's instruction. The DNA extract was checked on 1% agarose gel, and its concentration and purity were determined with a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, USA).

The hypervariable region V4 of the 18S rDNA was amplified with primer pairs V4-F (5′-GCGG TAATTCCAGCTCCAATA-3′) and V4-R (5′-GATC CCCHWACTTTCGTTCTTGA-3′) (Song et al., 2016) by 9700 PCR thermocycler (ABI, USA). The PCR amplification was performed as follows: initial denaturation at 94 ℃ for 3 min, followed by 33 cycles of denaturing at 94 ℃ for 30 s, annealing at 58 ℃ for 45 s and extension at 72 ℃ for 45 s, and final extension at 72 ℃ for 5 min, and end at 4 ℃. PCR reactions were performed in triplicate and the PCR products were extracted and purified using the DNA Gel Extraction Kit (Axygen Biosciences, USA), and then quantified using Quantus^TM Fluorometer (Promega, USA).

2.3 Sequencing and bioinformatics analysis

Purified amplicons were pooled in equimolar and paired-end sequenced on an Illumina NovaSeq PE250 platform (Illumina, USA) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (China). The sequences reported in this paper were deposited into the National Center for Biotechnology Information (NCBI) GenBank database (Accession Nos: ON129525–ON129547). The raw 18S rDNA sequence reads were deposited to NCBI under the project number PRJNA876366.

The raw 18S rDNA sequence reads were demultiplexed, quality-filtered by fastp version 0.20.0 (Chen et al., 2018), and merged by FLASH version 1.2.7 (Magoč and Salzberg, 2011). Only overlapping sequences longer than 10 bp were assembled and the maximum mismatch ratio of the overlap region is 0.2. Reads that could not be assembled were discarded. Samples were distinguished according to the barcode and primers, and the sequence direction was adjusted.

Then the high-quality sequences were de-noised using DADA2 (Callahan et al., 2016) plugin in the Qiime2 pipeline with recommended parameters (maxEE=2, truncQ=0). After that, the amplicon sequence variants (ASVs) were obtained and then annotated for species classification through the Protist Ribosomal Reference (PR²) database, setting the comparison threshold at 70%. For the ASVs that could not be annotated by the PR² database, a further taxonomic assignment was used by BLAST against the NCBI nucleotide database. Each ASV was annotated to the species level with the highest percentage identity (PID) at a cutoff level of 99%. The ASVs of Phaeocystis should meet the following conditions in sequence alignment with known Phaeocystis: Query cover=100%, E value=0, PID > 99%. ASVs information and sequence abundance of Phaeocystis are shown in Supplementary Table S1.

The Maximum Likelihood (ML) phylogenetic tree was constructed using the MEGA X (Kumar et al., 2018) software package. The phylogenetic tree was inferred by the maximum likelihood (ML) method, using a Kimura 2-parameter model of DNA substitution. The bootstrap consensus tree was inferred from 1 000 replicates. The sampling sites and distribution of species were drawn with Surfer 16 (Golden Software LLC, USA) and GraphPad Prism 8 (GraphPad Software, USA). Spearman correlation coefficient was used to explore the relationship between Phaeocystis abundance and environmental factors. The spearman corresponding significance was tested and drawn using the corrplot package in R. Community network analysis among the Phaeocystis and other plankton species was carried out using the Spearman correlation coefficient in the R package psych. Principal component analysis (PCA) of the environmental factors was performed using software Canoco 5 (Microcomputer Power, USA) and the data was lg(n+1) transformed prior to performing PCA.

3 RESULT 3.1 Environmental condition

In April, the surface water temperature varied from 6.4 ℃ to 24.7 ℃, with an average value of 13.2±4.2 ℃. It was generally high in the ECS, with the highest value (24.7 ℃) recorded at the southernmost station S55, and the lowest value (6.4 ℃) occurred in the BS at the northernmost station B09. The salinity ranged from 29.1 to 34.5, with an average value of 32.4±1.3 across the whole studied area. In general, the salinity of the nearshore is lower than that of the open sea (Fig. 2).

The concentration of nitrate, phosphate, and silicate showed a similar spatial distribution pattern, with the highest content occurring near the Changjiang (Yangze) River estuary and the ECS coast (Fig. 3a, d, & e). The nitrate concentration ranged from 0.02 to 15.12 μmol/L, with a mean of 2.40± 3.06 μmol/L. It was higher than 11.56 μmol/L around the Changjiang River estuary and higher than 6.63 μmol/L along the ECS coast. In addition, phosphate concentration ranged from 0.01 to 0.90 μmol/L (mean 0.21±0.17 μmol/L) and silicate concentration ranged from 0.52 to 21.20 μmol/L (mean 5.51±5.17 μmol/L). The nitrite concentrations were between 0.01 and 0.85 μmol/L (average 0.20± 0.20 μmol/L), with the high values generally found near the Changjiang River estuary and the ECS coast (Fig. 3b). The ammonium levels varied from 0.12 to 14.44 μmol/L, with an average of 1.64±2.15 μmol/L (Fig. 3c). Two maximum of ammonium concentrations were found at stations H05 and H15 of the southern Yellow Sea. In addition, chlorophyll-a content showed remarkable variability, with concentrations varying between 0.51 and 14.47 μg/L (mean 2.47±1.98 μg/L) (Fig. 3f).

3.2 Species diversity and distribution of Phaeocystis

In this study, a total of 26 381 ASVs were identified from the environment samples in April 2021, and 18 ASVs were annotated to four recognized Phaeocystis species (Figs. 4–5), including P. globosa, P. cordata, P. pouchetii, and P. jahnii using the DADA2 analysis based on partial 18S rDNA sequences. Among them, the warm-water species P. globosa distributed widely throughout the Chinese seas, especially in the BS (Figs. 4–5). The highest P. globosa abundance was found at Station H15 in the southern YS (Fig. 4a). P. cordata is also a widespread species and occurred at most sampling stations, with the higher abundances at Stations H19 and H37 in the southern YS (Fig. 4b). Cold-water species P. pouchetii distributed only in the BS and the northern YS, with a substantially elevated abundance occurring at Station B09 in the BS, which was much higher than that of other stations (Fig. 4c). The solitary species P. jahnii was only recorded in the southeastern ECS with relatively low abundance (Figs. 4d & 5).

Overall, Phaeocystis was found at 58 out of 82 sampling stations. From a biodiversity perspective, the number of Phaeocystis species was relatively low in the southern Yellow Sea, which was dominated by P. globosa in terms of abundance (Fig. 5). In addition, several ASVs from the southern Yellow Sea were annotated yet uncharacterized Phaeocystis species (Fig. 4e), which can be clustered into one distinguished clade (Fig. 6). Clade 1 is a sister taxon to the P. globosa, but is significantly different from the known P. globosa species, suggesting the existence of possible cryptic species of P. globosa.

3.3 Correlation between environmental factors and Phaeocystis

Principal component analysis of environmental factors basically divides the sampling stations into three groups (BS, YS, and ECS), which are consistent with their geographical areas (Fig. 7). The results indicate that 56.0% and 21.6% of the variation could be explained by PC1 and PC2, respectively. Samples of the ECS are closely related to high temperature, salinity, and high silicate, while samples from the BS showed opposite characteristics with low temperature, low salinity, but high ammonium and phosphate. No obvious features were observed for the YS's samples, with half of them characterized by low nutrients (Fig. 7). To further explore the impact of environmental factors on Phaeocystis abundance and distribution, the correlation between them was analyzed. The P. globosa abundance was positive significantly related to phosphate and ammonium concentrations (P < 0.05, n=82), but negatively related to the water temperature, salinity, and silicate concentration (P < 0.05, n=82), respectively. The correlation between the P. pouchetii and environmental factors was similar to those of P. globosa, but the significant levels were much higher, in particular the effect of water temperature and salinity (P < 0.001, n=82). On the contrary, P. jahnii was positively related to the water temperature and salinity (P < 0.05, n=82), and negatively related to phosphate and ammonium concentrations (P < 0.01, n=82).

In addition, community network analysis among all the eukaryotes species indicated the complicated relationships among the plankton, including competition and predation. In general, different Phaeocystis species showed negative response to ciliates, copepods, and some dinoflagellates and diatom species (e.g., Noctiluca sp. and Chaetoceros decipiens), but usually positive related to phytoplankton, in particular several clades of Syndiniales dinoflagellates (Syndiniales Groups Ⅰ and Ⅱ; Supplementary Table S2).

4 DISCUSSION 4.1 Phaeocystis species and their distribution

In the present study, 4 described Phaeocystis species, including P. globosa, P. pouchetii, P. cordata, and P. jahnii, were detected in the Chinese seas in April 2021. Except for P. globosa, the other 3 species are new records to several seas, implying the species number of the genus is severely underestimated on Chinese coasts (Shen and Qi, 2021). Although the South China Sea (SCS) was not investigated this time, previous studies have demonstrated that P. globosa distributed widely in the SCS (Hai et al., 2010; Shen et al., 2018). Most interestingly, two cold-water species P. antarctica and P. pouchetii were found in a cold eddy in the tropical SCS with high abundance (Wu et al., 2015). Moreover, P. scrobiculata was also observed on the Vietnam coast along with the SCS based on morphological features (Doan-Nhu and Larsen, 2010), revealing a much higher diversity of the Phaeocystis than previously recognized in the SCS, which is now recognized as a hotspot of Phaeocystis distribution and biodiversity (Shen and Qi, 2021).

There is no doubt that P. globosa is the most widespread and abundant species in the world's oceans, including off the China coasts (Schoemann et al., 2005, and references therein). From 1997 to 2022, P. globosa has caused over 100 bloom events in China, including 13 events in the Bohai Sea, 5 blooms in the East China Sea, and most remaining (over 80 events) in the northern South China Sea (e.g., Wang et al., 2021; Li et al., 2022). In late November 2021, a large-scale P. globosa bloom broke out in the southern Yellow Sea for the first time (Li et al., 2022), and just before that in April, the highest abundance of P. globosa was detected near this region at Station H15 (Fig. 4a), suggesting a great potential risk of the bloom outbreak. In addition, the ammonium concentration at this station and during the P. globosa bloom in the southern Yellow Sea in late November were both high, implying an important effect of the ammonium on the growth of P. globosa (this study; Wang et al., 2011; Li et al., 2022). Besides ammonium, phosphate also showed a positive significant influence on the P. globosa abundance in this study.

In fact, numerous studies have demonstrated that nutrients play an essential role in regulating the outbreak and maintenance of P. globosa blooms (Schoemann et al., 2005; Lancelot et al., 2007; Xu et al., 2019), but the results were depending and other effects such as temperature and induced allelopathy should be taken into account. Moreover, although no Phaeocystis blooms occurred in the Bohai Sea in recent years, P. globosa is still widely distributed in the Bohai Sea, indicating the high potential risk of P. globosa bloom reoccurrence in this area.

Phaeocystis pouchetii is another important bloom-forming species that has caused dense blooms in the high latitude northern seas, including the New England continental shelf, Massachusetts Bay, and the Arctic Ocean (Saiz et al., 2013; Borkman et al., 2016; Smith et al., 2021). It can survive at water temperatures between -1.7 and 9 ℃, but is better adapted to temperature below 5 ℃ (Schoemann et al., 2005; Simo-Matchim et al., 2017). In this study, extremely high abundance of P. pouchetii was found in the northern Yellow Sea and the Bohai Sea, where the water temperature was the lowest (6.4–9.6 ℃). Moreover, a strong negative correlation between P. pouchetii and temperature was detected, indicating that low temperature is critical to the distribution and abundance of P. pouchetii in the Chinese seas. On the other hand, based on the high abundance of P. pouchetii, it is reasonable to speculate that previously reported Phaeocystis blooms from the Bohai Sea were possibly initiated by P. pouchetii and/or P. globosa since only limited morphological and molecular evidence is available for most of them (Qu et al., 2008; Tu et al., 2011). Hence, identification of the causative species needs to be strengthened, especially in the Bohai Sea.

Two new record species, P. jahnii and P. cordata, were further observed in coastal China (Lin et al., 2014). Based on the phylogenetic analysis, P. jahnii and P. cordata are located in the base branch of the phylogenetic tree and differentiated early in this genus, which should be cosmopolitan species in the world's oceans (Zingone et al., 1999; Lange et al., 2002). The distribution of P. jahnii and P. cordata on the Chinese coasts confirmed the prediction of widespread species once again (Zingone et al., 1999). Furthermore, it is worthy to note that a bloom with a high density of P. jahnii and P. cordata (2.2×104 cells/mL) was observed in the ECS in the spring of 2009 (Lin et al., 2014). In this study, an elevated abundance of P. jahnii was found in the southeastern East China Sea, indicating a high potential risk of solitary cell bloom in this area, where the temporal and spatial distribution of P. jahnii may have been under-sampled. Moreover, P. jahnii was positively related to the water temperature and salinity, while negatively related to phosphate and ammonium concentrations. These features may further indicate that P. jahnii prefers warm open ocean areas; hence, the blooms of solitary cells have been reported on occasion (e.g., Wassmann et al., 2005; Lin et al., 2014).

4.2 Uncharacterized Phaeocystis species

Besides the above 4 described species, 2 uncharacterized Phaeocystis species have also been detected in the southern Yellow Sea. One is a possible cryptic species of P. globosa, and the other is a sister taxon to the species P. rex. In fact, uncharacterized Phaeocystis species are widely distributed in the world's oceans (Schoemann et al, 2005; Medlin and Zingone, 2007; Shen and Qi, 2021). For example, 3 undescribed species and 6 single-cell morphology from the Mediterranean Sea have been reported (Medlin and Zingone, 2007; Novarino, 2014). However, the common problem is that it is uncertain whether these cell types can represent biologically independent Phaeocystis species or the cell morphology of the same species at different life stages (Shen and Qi, 2021). Anyway, the existence of uncharacterized Phaeocystis species indicates that the diversity of the genus remains poorly understudied worldwide. It is partly due to the tiny size of the single-cell species of Phaeocystis, which is easy to be ignored in routine microscopic detection. At present, high-throughput sequencing of rRNA genes combined with metabonomics analysis can be used to investigate species diversity in a wide range of sea areas, especially single-celled algae with complex life cycles (Stoeck et al., 2010; De Vargas et al., 2015). Recently developed analysis methods based on amplicon sequence variants (ASVs) are more sensitive to distinguishing sequence variation. Therefore, the situation of the classification could be improved significantly in the future.

5 CONCLUSION

High-throughput sequencing data from 82 stations expanded our understanding of the species diversity and spatial distribution of Phaeocystis in coastal China. The results showed that at least 4 described Phaeocystis species, including P. globosa, P. pouchetii, P. cordata, P. jahnii, and 2 uncharacterized species, coexisted on the Chinese coast. P. globosa and P. cordata are eurythermal and euryhaline species with a wide distribution. The cold-water species P. pouchetii has also been found off the coast of China, particularly in the Bohai Sea, and P. jahnii prefers warm open ocean areas. Furthermore, environmental factors play an important role in affecting the abundance and geographic distribution of Phaeocystis, especially the temperature. More new species of Phaeocystis are yet to be discovered, and the metabarcoding analysis is an effective method for detecting species and profiling the biodiversity of the Phaeocystis.

6 DATA AVAILABILITY STATEMENT

The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

7 ACKNOWLEDGMENT

Samples of the BS and the YS were collected onboard R/V Lanhai 101 implementing the open research cruise NORC2021-01 supported by NSFC Shiptime Sharing Project (No. 42049901). Sample collections of the ECS were supported by NSFC Open Research Cruise (No. NORC2021-02+NORC2021-301). These cruises were conducted onboard R/V Xiangyanghong 18 by the First Institute of Oceanography, Ministry of Natural Resources, China.

Electronic supplementary material

Supplementary material (Supplementary Tables S1–S2) is available in the online version of this article at https://doi.org/10.1007/s00343-022-2256-1.