Phylogenetic and Phenogenetic Diversity of Synechococcus along a Yellow Sea Section Reveal Its Environmental Dependent Distribution and Co-Occurrence Microbial Pattern

: Synechococcus is a dominant genus of the coastal phytoplankton with an effective contribution to primary productivity. Here, the phylogenetic and phenogenetic composition of Synechococcus in the coastal Yellow Sea was addressed by sequencing marker gene methods. Meanwhile, its co-occurrence pattern with bacterial and eukaryotic microbes was further investigated based on the construction of networks. The result revealed that Synechococcus abundance ranged from 9.8 × 10 2 cells mL − 1 to 1.6 × 10 5 cells mL − 1 , which was signiﬁcantly correlated to sampling depth and nutrient contents of nitrite, ammonia, and dissolved silicon. A total of eight Synechococcus phylogenetic lineages were detected, of which clade III was dominant in most of the samples. Meanwhile, clade I increased along the water column and even reached a maximum value of 76.13% at 20 m of station B. Phenogenetically, Synechococcus PT3 was always the predominant pigment type across the whole study zone. Only salinity was signiﬁcantly correlated to the phenogenetic constitution. The networks revealed that Synechococcus co-occurred with 159 prokaryotes, as well as 102 eukaryotes including such possible grazers as Gymnodinium clades and Alveolata. Potential function prediction further showed that microbes co-occurring with Synechococcus were associated with diverse element cycles, but the exact mechanism needed further experimentation to verify. This research promotes exploring regularity in the genomic composition and niche position of Synechococcus in the coastal ecosystem and is signiﬁcant to further discuss its potential participation in materials circulation and bottom-up effects in microbial food webs. constraints; and (iii) study the co-occurrence pattern of Synechococcus with multiple eukaryotic and prokaryotic taxa to explore its ecological functions associated with microbiome Synechococcus in the coastal ecosystem.


Introduction
Marine Synechococcus is a main group of photosynthetic picophytoplankton with the feature of self-assembly orange phycoerythrin for light capturing. Taxonomically, it belongs to Cyanobacteria and is demonstrated with rich phylogenetic diversity [1]. Through molecular detections, Synechococcus can be classified into three phylogenetic subclusters including subcluster 5.1 (S5.1), subcluster 5.2 (S5.2), and subcluster 5.3 (S5.3) [2], which at least consist of 20 clades, 2 clades, and 6 clades, respectively [3]. Furthermore, the different constitution of the main light-harvesting antennae for photosystem II, phycobilisomes (PBSs), leads to its rich phenogenetic diversity [4]. Three main pigment types (PTs), PT1, PT2, and PT3, can be identified [5]. Due to the independent evolution between phylogenetics and PTs, it is impossible to simultaneously identify these two assemblage characteristics

Sampling
Field sampling was conducted in the South Yellow Sea by taking Research Vessel Kexue III from August 29 to September 6 in 2018 ( Figure 1). Three sample stations on the same latitude of 36 • N were designated according to their representatives at relatively shallow (29 m), average (43 m), and relatively deep (72 m) depths. Fifteen samples were collected from standard depths according to the specifications for the oceanographic survey (National Standard of China, GB/T 12763-2007). A total of 1.40 mL seawater was taken in sterile tubes in quintuplicate, and then 0.20 mL paraformaldehyde (final concentration 0.5%) was added to determine the count of Synechococcus [39]. The tubes were flash-frozen in liquid nitrogen, and then stored at −80 • C until analysis using flow cytometry. One liter seawater was filtered through 48 µm nylon mesh to remove microphytoplankton, and then through 0.22 µm polycarbonate membrane (Millipore Co., Bedford, MA, USA). Acquired membranes were preserved for high-throughput sequencing analysis by storing in liquid nitrogen. J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 3 of 17

Sampling
Field sampling was conducted in the South Yellow Sea by taking Research Vessel Kexue III from August 29 to September 6 in 2018 ( Figure 1). Three sample stations on the same latitude of 36 °N were designated according to their representatives at relatively shallow (29 m), average (43 m), and relatively deep (72 m) depths. Fifteen samples were collected from standard depths according to the specifications for the oceanographic survey (National Standard of China, GB/T 12763-2007). A total of 1.40 mL seawater was taken in sterile tubes in quintuplicate, and then 0.20 mL paraformaldehyde (final concentration 0.5%) was added to determine the count of Synechococcus [39]. The tubes were flash-frozen in liquid nitrogen, and then stored at −80 °C until analysis using flow cytometry. One liter seawater was filtered through 48 μm nylon mesh to remove microphytoplankton, and then through 0.22 μm polycarbonate membrane (Millipore Co., Bedford, MA, USA). Acquired membranes were preserved for high-throughput sequencing analysis by storing in liquid nitrogen. Physical parameters and chlorophyll a (chl. a) content were detected on boat through a CTD sensor (Sea-Bird Electronics Inc., Bellevue, WA, USA). The nutrient concentrations were determined by an AA3 segmented flow analyzer (Seal Analytical GmbH, Germany) [40].

Cell Abundance of Synechococcus
Synechococcus abundance was detected using a BD FACSAria™ flow cytometer. Yellow-green fluorescent beads (2.0 μm, Polysciences, Warrington, PA, USA) were added to each sample before the detection. The forward and right-angle light scattering and four fluorescence signals were collected and analyzed using WinMDI 2.9. Red fluorescence excited by a 488 nm laser indicates the presence of chl. a. Red fluorescence excited by 635 nm indicates the presence of phycocyanin (PC), while orange fluorescence indicates the presence of phycoerythrin (PE).

DNA Extraction
Total nucleic acids were extracted from the membrane bead-beaten with 500 mg of 0.10 mm silica/zirconia glass beads (BioSpec Products Inc., Bartlesville OK, USA), 0.50 mL 120 mM (pH 8.0) K2HPO4 butter, and 0.50 mL phenol-chloroform-isoamyl alcohol (25:24:1) (pH 8.0) by Mini-BeadBeater (BioSpec, USA). The supernatants were collected to extract DNA using the Fast DNA spin kit (MP BIO, Santa Ana, CA, USA) following the instruction of the kit. Agarose gel electrophoresis and the Nano Drop 2000c spectrophotometer Physical parameters and chlorophyll a (chl. a) content were detected on boat through a CTD sensor (Sea-Bird Electronics Inc., Bellevue, WA, USA). The nutrient concentrations were determined by an AA3 segmented flow analyzer (Seal Analytical GmbH, Germany) [40].

Cell Abundance of Synechococcus
Synechococcus abundance was detected using a BD FACSAria™ flow cytometer. Yellowgreen fluorescent beads (2.0 µm, Polysciences, Warrington, PA, USA) were added to each sample before the detection. The forward and right-angle light scattering and four fluorescence signals were collected and analyzed using WinMDI 2.9. Red fluorescence excited by a 488 nm laser indicates the presence of chl. a. Red fluorescence excited by 635 nm indicates the presence of phycocyanin (PC), while orange fluorescence indicates the presence of phycoerythrin (PE).

DNA Extraction
Total nucleic acids were extracted from the membrane bead-beaten with 500 mg of 0.10 mm silica/zirconia glass beads (BioSpec Products Inc., Bartlesville, OK, USA), 0.50 mL 120 mM (pH 8.0) K 2 HPO 4 butter, and 0.50 mL phenol-chloroform-isoamyl alcohol (25:24:1) (pH 8.0) by Mini-BeadBeater (BioSpec, USA). The supernatants were collected to extract DNA using the Fast DNA spin kit (MP BIO, Santa Ana, CA, USA) following the instruction of the kit. Agarose gel electrophoresis and the Nano Drop 2000c spectrophotometer (Thermo Fisher, Waltham, MA, USA) were used to assess the quality and concentration of the extracted DNA. Then, extracted DNA was kept at −20 • C until further analysis.

PCR Reactions
Polymerase chain reaction (PCR) was performed in the thermocycler GeneAmp ® PCR system 9700 (Perkin Elmer, Foster City, CA, USA) in triplicate [41]. The rpoC1 gene of Synechococcus was amplified by the primers rpoC1-39F (5 -GGNATYGTYTGYGAGCGYTG-3 ) and rpoC1-462R (5 -CGYAGRCGCTTGRTCAGCTT-3 ) with a unique 6 nt barcode attached at the end of the forward primer [42]. The PCR reactions for rpoC1 gene were conducted with a 2 × Rapid Taq Master Mix (Vazyme, Nanjing, China) in a 20 µL reaction volume containing 20 ng template DNA. The PCR process started with a 3-min pre-denaturation at 95 • C, followed with 37 cycles of denaturation at 95 • C for 30 s, annealing at 55 • C for 30 s, and elongation at 72 • C for 45 s, and ended with a final extension at 72 • C for 10 min. The cpeBA operon of Synechococcus was amplified in the same way but by the primers peBF (5 -GACCTACATCGCWCTGGGYG-3 ) and peAR (5 -CCMACAACCARGCAGTAGTT-3 ) with annealing at 53 • C [16]. The V4 regions of the 16S RNA gene and 18S RNA gene were amplified with the same 20 µL reaction volume by the universal prokaryotic primers 515F/806R and the specific eukaryotic planktonic primers 547F (5 -CCAGCASCYGCGGTAATTCC-3 )/952R (5 -ACTTTCGTTCTTGATYRA-3 ), respectively [43,44]. The PCR process for 18S started with a 2-min pre-denaturation at 94 • C, followed by 10 cycles of denaturation at 94 • C for 30 s, annealing at 65 • C for 30 s, 18 cycles of denaturation at 94 • C for 30 s, annealing at 55 • C for 30 s and elongation at 72 • C for 30 s, and ended with a final extension at 72 • C for 10 min [44].

Sequencing
For the rpoC1, 16S, and 18S genes, high-throughput sequencing for the PCR products after being gel purified was conducted at Majorbio Co., Ltd. (Shanghai, China). Sequencing libraries were generated using NEB Next ® Ultra™DNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) following the manufacturer's recommendations, and index codes were added. The library quality was assessed on the Qubit@ 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. At last, the library was sequenced on an Illumina MiSeq platform and 300 bp paired-end reads were generated. For the cpeBA gene, PCR products were purified by agarose gel electrophoresis with a TaKaRa MiniBEST DNA Fragment Purification Kit (TaKaLa, Kusatsu, Japan) and then cloned into a pMD-19T vector with a Cloning Kit (TaKaLa, Dalian, China). The recombinant products were subsequently transformed into Escherichia coli DH5α competent cells. After incubation, recombinant clones were picked up and sequenced by the RuiBo Sequencing Company (Qingdao, China) for constructing clone libraries.

Sequence Data Processing
The raw data were quality-controlled [45], chimeras-filtered [46], and redundancyeliminated [47] to downstream analyses using QIIME2. Sequences of rpoC1 and cpeBA genes were aligned using Clustal W [48] and trimmed by BioEdit [49]. Sequences were grouped into OTUs (operational taxonomic units) at 95% similarity. Representative sequences of OTU were first identified using BLASTn against the nt database, and those not belonging to Cyanobacteria were picked out. The 15 most abundant OTUs sequenced based on the rpoC1 gene (covering 96.3% of total reads) and cpeBA operon (covering 94.2% of total reads) were aligned with the reference sequences (Tables S1 and S2). Maximum likelihood phylogenetic trees were constructed by Mega 7.014 with the model GTR+G+I and 200 bootstraps [9]. For 16S and 18S sequences, the OTU clustering was performed considering a similarity of 97% to obtain the representative OTU sequences [46]. The Greengenes v.135/16S database and the Silva v.132/18S eukaryota database were used to obtain the taxonomic information corresponding to each OTU. The function of 16S sequencing data was annotated using the script collapse_table.py based on the Functional Annotation of Prokaryotic Taxa (FAPROTAX) database in Python v.2.7 [41,50].

Construction of Co-Occurrence Network
Rare OTUs of 16S and 18S genes with a total relative abundance < 0.01% in all samples were filtered out prior to the construction of microbial community networks [29]. Spearman correlation matrixes among OTUs were calculated based on the relative abundance [51,52]. The connections in the network with the absolute value of Spearman coefficients > 0.6 and p values < 0.05 were considered significant [53]. Network properties were calculated using the psych package [54], and the networks were drawn using Gephi v0.9.2 [39,55]. Global genetic and functional correlation networks encompassing all genera and functions were constructed from the correlation matrix. Metabolic cycling networks were constructed using similar methods based on the relative abundance of metabolic functions annotated by the FAPROTAX database. To determine whether the proportion and type of metabolic cycles involving Synechococcus vary along the vertical water column, the samples were selected to represent the upper (A1, A2, B1, B2, and C1) and bottom (A5, B5, C4, and C5) layers according to the depths and photosynthetically active radiation (PAR) [56].

Statistical Analysis
Alpha diversity indexes and beta diversity distance were calculated using the core_diversity_analyses.py script. The environmental data were normalized and scaled to statistical analysis. The vegan package v.2.5-6 in R v3.6.1 was used to calculate the variance inflation factor (VIF), and then we picked nonredundant environmental variables to conduct redundancy analysis (RDA) with OTU tables [57]. Statistical parameters of co-occurrence networks, such as modularity and density, were calculated using the builtin program of Gephi v0.9.2 [58]. Correlationships among environmental variables and between environmental variables and Synechococcus cell abundance were examined by Spearman's correlations using the corrplot package v.0.84 in the R programming language v3.6.1. Differences in Synechococcus cell abundance among three stations were tested using analysis of variance (ANOVA) with a post hoc test (Tukey test) in the software Statistical Product and Service Solutions (SPSS v.17.0). A p-value less than 0.05 denoted statistical significance. The software TBtools (Toolbox for Biologists) v1.055 [59], and the vegan package v.2.5-6 in R v3.6.1 was utilized to draw a heatmap.

Environmental Variables
The environmental variables of each sample are shown in Table 1, and their Spearman's correlations are shown in Table S3. The water depth increases from the shallow coastal station A to C. Positive correlations between longitude and oxygen (ρ = 0.89) and with salinity (ρ = 0.59) were found. Along the vertical water column with the increase in sampling depth, salinity (ρ = 0.84), the concentration of dissolved total phosphorus (DTP; ρ = 0.85), dissolved silicon (DSi; ρ = 0.78), and phosphate (PO 4 3− ; ρ = 0.73) gradually increased, whereas temperature (ρ = −0.91) and concentration of ammonia (NH 4 + ; ρ = −0.75) continually decreased. In particular, the presence of YSCWM (Yellow Sea Cold Water Mass) caused a sharp drop in temperature in the lower water layers. Correlations between nutrient concentrations were also revealed (Table S3).

Phylogenetic and Phenogenetic Diversity and Structure of Synechococcus Assemblage
The rpoC1 gene generated 216,077 quality sequences from 14 samples, with an average of 15,434 sequences per sample (Table S4). The good coverage of all rpoC1 sequencing samples was >99.9%. The genomic diversities of Synechococcus assemblages were estimated by the NP-Shannon index and the Chao1 index (predicted OTUs) ( Table S4). The Synechococcus counts showed no significant difference among the three stations. However, they correlated with the vertically sampled water layer significantly (ρ = 0.62; Figure 2). The concentrations of nitrite (NO 2 − ; ρ = 0.58) and DSi (ρ = 0.57) were positively correlative to the cell abundance of Synechococcus, while NH 4 + (ρ = −0.66) was negatively correlative to the Synechococcus counts.

Phylogenetic and Phenogenetic Diversity and Structure of Synechococcus Assemblage
The rpoC1 gene generated 216,077 quality sequences from 14 samples, with an average of 15,434 sequences per sample (Table S4). The good coverage of all rpoC1 sequencing samples was >99.9%. The genomic diversities of Synechococcus assemblages were estimated by the NP-Shannon index and the Chao1 index (predicted OTUs) ( Table S4). The high phylogenetic diversity usually appeared in the middle or bottom layers along the water column of stations B and C, while the surface layer of station A had higher diversity indexes. In the phylogenetic tree, a total of eight Synechococcus lineages were identified from 14 samples based on the rpoC1 gene (Figure 3a). Three subclusters of Synechococcus, S5.1, S5.2, and S5.3, all appeared in the studied area. Among them, S5.1 accounted for the highest proportion, with an average of 90.92%, followed by S5.3 (2.09%), while S5.2 as a rare subcluster accounted for only 0.54%. Both S5.1 and S5.2 were more abundant in station C, while S5.3 had a higher relative abundance in station B. Synechococcus clade III of S5.1 was the dominant clade, which possesses a high proportion of up to 60.82%, on average. In particular, it occupied more than 80% of the total clades in station A, which was about twice that in stations B and C. The proportion of clade I increased with the decrease in water depth and reached a maximum value of 76.13% at B4. Clade IV in station A was significantly less than that of station C and the lower water areas of station B. A higher relative abundance of clades II and IX was detected in station A, while that of clade VI was in station C. The proportion of Synechococcus S5.2 was not obviously different among samples, while that of S5.3 was relatively high in the surface area. The hierarchical clustering based on Bray-Curtis demonstrated that the phylogenetic composition of Synechococcus in samples from lower water layers of stations B and C was different from others ( Figure 3a).  Figure 3a). The cpeBA gene generated 497 quality sequences from 13 samples, with an average of 38 sequences per sample, through the method of clone library (Table S4). The good coverage of the cpeBA gene was >73.0%, except in A2 (only 1 OTU). In the tree based on the cpeBA operon, five well-separated clusters were formed (Figure 3b). PT2 and 3a Synechococcus could be classified by the sequencing of the cpeBA sequence, while Synechococcus PT3cB, 3d, and 3eA could not be distinguished from each other. PT3a sequences from clade I and V formed a cluster (hereafter PT3a-I/V) and were separated from the cluster formed by those from clade II and WPC1 (hereafter PT3a-II/WPC1). The sample A2 only The cpeBA gene generated 497 quality sequences from 13 samples, with an average of 38 sequences per sample, through the method of clone library (Table S4). The good coverage of the cpeBA gene was >73.0%, except in A2 (only 1 OTU). In the tree based on the cpeBA operon, five well-separated clusters were formed (Figure 3b). PT2 and 3a Synechococcus could be classified by the sequencing of the cpeBA sequence, while Synechococcus PT3cB, 3d, and 3eA could not be distinguished from each other. PT3a sequences from clade I and V formed a cluster (hereafter PT3a-I/V) and were separated from the cluster formed by those from clade II and WPC1 (hereafter PT3a-II/WPC1). The sample A2 only contained 1 OTU, so it was removed from statistical analyses. Synechococcus PT3, including PT3a-I/V, PT3a-II/WPC1, and PT3cB/3d/3eA, was dominant in the studied area with a mean proportion of 73.04% in all samples, while that of PT2, including PT2A and PT2B, decreased to 13.42%. Among PT3, PT3cB/3d/3eA had the highest relative abundance in the studied section, with a peak of 52.94% in C2. PT3a-I/V had a higher proportion in station C, while PT3a-II/WPC1 was more abundant in station A. For PT2, the relative abundance of PT2A was relatively higher than PT2B. In samples C2 and C3, both PT2A and PT2B were not detected. Furthermore, PT2B was not found in B4 and most water layers of station C. The clustering analysis showed that stations A and C each formed a separated cluster, while samples from station B alternatively clustered in them (Figure 3b).
The first axis and the second axis of the RDA analysis based on the rpoC1 gene explained 61.33% and 10.15% of the variance in Synechococcus assemblages, respectively ( Figure 4). Five environmental variables, NO 2 − (ρ = 0.60), oxygen (ρ = 0.58), sampling depth (ρ = 0.56), DTN (ρ = 0.54), and NO 3 − (ρ = 0.49), significantly correlated to the community. Clade I and IV relative abundances were associated with nitrite/nitrate-N, while clade III showed correlations with oxygen. Clade VI was associated with sampling depth. The first axis and the second axis of the RDA analysis based on the cpeBA gene totally explained 61.65% of the variance of Synechococcus assemblage. Only salinity (ρ = 0.64) was significantly correlated to the phenotypic constitution. contained 1 OTU, so it was removed from statistical analyses. Synechococcus PT3, including PT3a-I/V, PT3a-II/WPC1, and PT3cB/3d/3eA, was dominant in the studied area with a mean proportion of 73.04% in all samples, while that of PT2, including PT2A and PT2B, decreased to 13.42%. Among PT3, PT3cB/3d/3eA had the highest relative abundance in the studied section, with a peak of 52.94% in C2. PT3a-I/V had a higher proportion in station C, while PT3a-II/WPC1 was more abundant in station A. For PT2, the relative abundance of PT2A was relatively higher than PT2B. In samples C2 and C3, both PT2A and PT2B were not detected. Furthermore, PT2B was not found in B4 and most water layers of station C. The clustering analysis showed that stations A and C each formed a separated cluster, while samples from station B alternatively clustered in them ( Figure  3b).
The first axis and the second axis of the RDA analysis based on the rpoC1 gene explained 61.33% and 10.15% of the variance in Synechococcus assemblages, respectively ( Figure 4). Five environmental variables, NO2 − (ρ = 0.60), oxygen (ρ = 0.58), sampling depth (ρ = 0.56), DTN (ρ = 0.54), and NO3 − (ρ = 0.49), significantly correlated to the community. Clade I and IV relative abundances were associated with nitrite/nitrate-N, while clade III showed correlations with oxygen. Clade VI was associated with sampling depth. The first axis and the second axis of the RDA analysis based on the cpeBA gene totally explained 61.65% of the variance of Synechococcus assemblage. Only salinity (ρ = 0.64) was significantly correlated to the phenotypic constitution. The Spearman correlation analysis was performed to examine the links between Synechococcus clades and PTs ( Figure S1). It showed that S5.1 clade I and IV as well as S5.2 positively correlated with PT3a-I/V (low phycourobilin: phycoerythrobilin), while clade III negatively correlated with it. No other significant correlation between Synechococcus clades and PTs was detected. The result supported the independent evolution between phylogenetic and phenogenetic diversity of Synechococcus.

Microbial Co-Occurrence with Synechococcus
The 16S rRNA gene generated 840,321 quality sequences from 15 samples, with an average of 56,021 sequences per sample. The 18S rRNA gene generated 850,827 quality sequences from 15 samples, with an average of 56,722 sequences per sample. Good coverages of all sequencing samples were >99.9%. The average value of the NP-Shannon in prokaryotic communities was 3.70, while the Chao1 index was 741 (Table S5). For eukaryotic communities, the average value of the NP-Shannon index was 3.65, and the Chao1 index was 466. The composition of the prokaryotic and eukaryotic communities was shown in Figure S2. The relative abundance of Synechococcus in the total prokaryotic community was high, maximally reaching 47.94% in A3 ( Figure S3). Six OTUs were detected The Spearman correlation analysis was performed to examine the links between Synechococcus clades and PTs ( Figure S1). It showed that S5.1 clade I and IV as well as S5.2 positively correlated with PT3a-I/V (low phycourobilin: phycoerythrobilin), while clade III negatively correlated with it. No other significant correlation between Synechococcus clades and PTs was detected. The result supported the independent evolution between phylogenetic and phenogenetic diversity of Synechococcus.

Microbial Co-Occurrence with Synechococcus
The 16S rRNA gene generated 840,321 quality sequences from 15 samples, with an average of 56,021 sequences per sample. The 18S rRNA gene generated 850,827 quality sequences from 15 samples, with an average of 56,722 sequences per sample. Good coverages of all sequencing samples were >99.9%. The average value of the NP-Shannon in prokaryotic communities was 3.70, while the Chao1 index was 741 (Table S5). For eukaryotic communities, the average value of the NP-Shannon index was 3.65, and the Chao1 index was 466. The composition of the prokaryotic and eukaryotic communities was shown in Figure S2. The relative abundance of Synechococcus in the total prokaryotic community was high, maximally reaching 47.94% in A3 ( Figure S3). Six OTUs were detected belonging to the genus Synechococcus, which dominated the phylum Cyanobacteria with 46.72-99.22% occupation ( Figure S4). No OTU was categorized in Prochlorococcus.
belonging to the genus Synechococcus, which dominated the phylum Cyanobacteria with 46.72-99.22% occupation ( Figure S4). No OTU was categorized in Prochlorococcus.
In the co-occurrence network, genetic correlation linked 869 nodes via 37,072 connections, corresponding to 30,707 positive and 6365 negative correlations ( Figure S5). Synechococcus was connected to 159 prokaryotic OTU nodes, including 5 archaea and 154 bacteria ( Figure 5). Proteobacteria was the most abundant prokaryotic group (76 OTUs), followed by Bacteroidetes (27), Actinobacteria (11), and Planctomycetes (11). A total of 102 of the 18S OTUs were associated with Synechococcus, which can be classified into eight known kingdoms (Alveolata, Chloroplastida, Cryptophyceae, Fungi, Haptophyta, Metazoa Animalia, Rhizaria, and Stramenopiles; Figure 5; Table S6). Among these kingdoms, Alveolata was represented the most (46), which included three phyla of Protalveolata (19), Dinoflagellata (17), and Ciliophora (10). The second-most abundant kingdom was Stramenopiles (21), which contained eight known phyla associated with Synechococcus. Figure 5. Microbial community co-occurrence network associated with Synechococcus. (a) The network of microbial community relationships encompassing all prokaryotic and eukaryotic OTUs was constructed from the correlation matrix using a force-directed layout algorithm; red nodes represent six OTUs belong to Synechococcus; gray nodes represent OTUs with the absolute value of Spearman coefficients < 0.6 and p > 0.05 with Synechococcus; green nodes represent prokaryotic OTUs with the absolute value of Spearman coefficients > 0.6 and p < 0.05 with Synechococcus; blue nodes represent eukaryotic OTUs with the absolute value of Spearman coefficients > 0.6 and p < 0.05 with Synechococcus. (b) Heatmap showing the number of predominant kingdoms/phyla associated with Synechococcus using color legend; two domains were represented with different colors (prokaryota: green; eukaryota: blue).

Distribution Features of Synechococcus Cell Abundance
In the studied Yellow Sea area, Synechococcus counts maximally reach 1.6 × 10 5 cells mL −1 , with an average of 2.4 × 10 4 cells mL −1 (Figure 2). Along the seawater column, the counts of Synechococcus increased first and then decreased with the maximum cell abundance appearing in the middle layer as A3 (10 m, 3.3 × 10 4 cells mL −1 ), B4 (20 m, 3.2 × 10 4 cells mL −1 ), and C4 (30 m, 1.6 × 10 5 cells mL −1 ). The count of Synechococcus in the studied area generally accorded with that in other coastal seas of the world: using the parametric model, the global quantitative distribution of Synechococcus reached 10 4 cells mL −1 except for the Arctic Ocean and the Southern Ocean [60]. The similar distribution characteristics Figure 5. Microbial community co-occurrence network associated with Synechococcus. (a) The network of microbial community relationships encompassing all prokaryotic and eukaryotic OTUs was constructed from the correlation matrix using a force-directed layout algorithm; red nodes represent six OTUs belong to Synechococcus; gray nodes represent OTUs with the absolute value of Spearman coefficients < 0.6 and p > 0.05 with Synechococcus; green nodes represent prokaryotic OTUs with the absolute value of Spearman coefficients > 0.6 and p < 0.05 with Synechococcus; blue nodes represent eukaryotic OTUs with the absolute value of Spearman coefficients > 0.6 and p < 0.05 with Synechococcus. (b) Heatmap showing the number of predominant kingdoms/phyla associated with Synechococcus using color legend; two domains were represented with different colors (prokaryota: green; eukaryota: blue).

Distribution Features of Synechococcus Cell Abundance
In the studied Yellow Sea area, Synechococcus counts maximally reach 1.6 × 10 5 cells mL −1 , with an average of 2.4 × 10 4 cells mL −1 (Figure 2). Along the seawater column, the counts of Synechococcus increased first and then decreased with the maximum cell abundance appearing in the middle layer as A3 (10 m, 3.3 × 10 4 cells mL −1 ), B4 (20 m, 3.2 × 10 4 cells mL −1 ), and C4 (30 m, 1.6 × 10 5 cells mL −1 ). The count of Synechococcus in the studied area generally accorded with that in other coastal seas of the world: using the parametric model, the global quantitative distribution of Synechococcus reached 10 4 cells mL −1 except for the Arctic Ocean and the Southern Ocean [60]. The similar distribution characteristics of Synechococcus abundance that ranged from 10 2 to 10 4 cells mL −1 in the coastal Yellow Sea waters were also observed [61]. The YSCWM-a bulk of cold water below the seasonal thermocline with a lower temperature (5-12 • C) and relatively constant salinity (31.5-32.5 PSU)-is an important phenomenon in the Yellow Sea, which has an impact on the summer strategy and vertical distribution of zooplankton and phytoplankton [62,63]. The marked decrease in Synechococcus counts in lower layers, particularly in station C, may attribute to the existence of YSCWM. Sequencing results also demonstrated that Synechococcus was dominant in the prokaryotic community with the absence of Prochlorococcus. Consequently, as the only highly abundant picocyanobacteria, Synechococcus is an important prokaryotic contributor to primary production in the studied coastal Yellow Sea.
The distribution of Synechococcus was significantly correlated to the sampling depth and contents of nitrite, ammonia, and dissolved silicon in the studied area. Consistent with the result of previous studies, the distribution of Synechococcus was affected by the sampling depth (Table 1) [64,65]. It could be closely related to the vertical attenuation in PAR, which is known to significantly affect the physiological activity of Synechococcus, including the amount of chl. a, total carotenoids, and D1 protein in the photosystem II reaction center [65,66]. The correlation between Synechococcus and nitrogen source (NO 2 − and NH 4 + ) was also worthy of investigation. Marine Synechococcus isolates have been reported to utilize NH 4 + , NO 2 − , NO 3 − , urea, and amino acids as nitrogen sources via systems facilitating the active uptake of these compounds [67]. This effect of the nitrogen source was considered to be related to the light-harvesting pigment protein, phycoerythrin, of Synechococcus [68]. Nitrogen shortage was another factor affecting the vertical distribution of Synechococcus, as the low light levels in the deep euphotic zone cannot provide sufficient energy for assimilating NO 2 − and NH 4 + for their N-rich phycobilisomes [69]. However, planktons prefer to incorporate NH 4 + rather than nitrate in the coastal ecosystem [70]. The enrichment of NH 4 + may attract relatively larger planktons such as diatoms, which may compete for nutrients and niches with Synechococcus. Some of these larger planktons may even prey on Synechococcus, which may explain the negative correlation between NH 4 + concentration and Synechococcus counts. Our results also identified a correlation between Synechococcus and another material, DSi, which has been shown to favor its biomass [23,71,72].

Structure and Niche Occupation of Synechococcus Assemblages
Despite the ubiquitous and cosmopolitan distribution of the genus Synechococcus, there might be unseen barriers among the world's oceans, and thus diverse genotypes of Synechococcus evolve. Some genotypes widely distribute, and others might restrict their habitat to narrow areas of marine environments, resulting in the current world distribution pattern of Synechococcus. Consistent with the results of Xia, Synechococcus S5.1 was the most abundant and widespread subcluster in the coastal Yellow Sea [8]. S5.2 and S5.3 were also widespread but occurred at low abundance. Generally, I, II, and IV are the most abundant in the Western Pacific Ocean, and the global ocean in general, in which clades I and IV appear to be the dominant type in temperate, coastal environments, while clade II widely dominated in subtropical/tropical and off-shore, oligotrophic environments [73,74]. Clade III also occurred in warm waters, but with a lower abundance and a narrower distribution [8]. However, in our study, clade III dominated in most samples, while clade I, II, and IV added up to one third of clade III on average. Another study by Xia showed that clade III is highly abundant in samples when the salinity is intermediate to high [9]. Tai and Palenik deemed the presence of clade III in this coastal environment to be more reflective of physical oceanographic patterns, such as stratification and the Inshore Countercurrent [75]. Synechococcus clade I was generally typical in cold, mesotrophic, and eutrophic waters [74]. Our research also exhibited the increase in clade I along with the vertical water column, especially at the bottom area located in the YSCWM zone.
The co-occurrence of pigment types in the studied coastal Yellow Sea was observed, and PT3 was always the dominant PT across the whole study area. This result is consistent with the opinion that one PT usually predominates in a sea area, although different PTs can co-occur [76]. The clone library method may not detect the rare taxa, while it is available to respond to the dominant species [16]. Attempts were made to detect the cpcBA operon, which encodes the phycocyanobilin (PCB). However, these attempts were not successful in any samples using the primer set of SyncpcB-Fw and SyncpcA-Rev [77]. Our gene marker cannot identify PC-only Synechococcus without the cpcBA operon, so cannot study the niche occupation of PT1 in the studied area. It was accepted that PT1 Synechococcus was found in turbid estuarine waters [78], while it was almost absent from the oceanic waters [10]. PT2 thrives in turbid coastal or continental shelf waters [79], whereas PT3 would mainly occur in mesotrophic or oligotrophic oceanic waters [9]. Consistently, it was ensured that from station A to C, the proportion of Synechococcus containing phycourobilin (PUB) increased gradually, while that containing no PUB decreased and was even absent in some water layers, demonstrating that underwater light spectral properties have a strong selective pressure on Synechococcus populations.

Potential Ecological Functions of Synechococcus in the Coastal Ecosystem
Synechococcus may play roles in multiple metabolic cycles in ecosystems via interactions with various microbial assemblages. From network results, Synechococcus was directly correlated with 159 prokaryotes and 102 eukaryotes ( Figure 5). Correlations between Synechococcus and other prokaryotes could be explained by a similar adaptability to the external environment, reciprocity on metabolite supply, competition for nutrients, and the antergic effect of pathogenic bacteria. Microbiota roles in biochemical cycles can be revealed from genomics-based analyses [50]. From genomic data, Synechococcus is found to be possibly involved with some metabolic activities, such as photoautotrophy [80], nitrogen fixation [81], and extracellular degradation of hydrocarbons ( Figure S6) [82]. The connection of Synechococcus with other metabolic activities, such as manganese oxidation and dark hydrogen oxidation, may result in the exchange of metabolic products (oxygen, carbon dioxide, and ammonia). Based on these multiple metabolic activities, Synechococcus appears to participate in carbon, nitrogen, hydrogen, and manganese cycling, outlining its feedback to the ecosystem. In addition, predation-based interactions can contribute to the transfer of materials and energy assimilated by Synechococcus to higher trophic levels.
In coastal areas, the environmental factors and biogeochemical cycles in different water layers vary, caused by the attenuation of light, exchanges between the atmosphere and ocean, and an inflow of freshwater rivers [56,83,84]. Using the function prediction of FAPROTAX, Synechococcus was found to associate with microbiota that are involved in different metabolic pathways and element cycles in upper and bottom water layers (Figures 6 and S6). In addition to photosynthesis and nitrogen fixation, Synechococcus was associated with the degradation of certain hydrocarbons carried by fresh rivers in upper layers. Comparatively, in the bottom layer, Synechococcus was related to more diverse microbial metabolic cycles, including carbon, nitrogen, hydrogen, and manganese cycles. Among these, 77.85% of the metabolic activities were aerobic chemoheterotrophy, which was on account of the supply of oxygen, and the decomposition of its biomass. In addition, previous studies have suggested that the prasinophyte Tetraselmis sp. (belonging to the Alveolata kingdom) and Gymnodinium-like dinoflagellates (belonging to the Chloroplastida kingdom) were Synechococcus grazers under certain conditions, and play important roles in determining the fluctuation of the Synechococcus population [85]. Here, significant correlations between 102 eukaryotes and Synechococcus were revealed, which included five Gymnodinium clades. No significant correlation was observed between Synechococcus and Tetraselmis sp., while the association with 46 extra Alveolata was observed. These eukaryotes may be potential grazers of Synechococcus, which could transfer biomass assimilated by Synechococcus to upper trophic levels of the food chain for maintaining the stability and coordination of the ecosystem (bottom-up effects on the food web; Figure 6).

Conclusions
The distribution of Synechococcus was studied in a typical section of the Yellow Sea. The cell abundance in the studied area was 2.4 × 10 4 cells mL −1 on average, which was significantly correlated to the sampling depth, nitrogen source (NO2 − and NH4 + ), and DSi. The intense decrease in the count in the lower layer probably attributes to the YSCWM. Synechococcus clade III was dominant in most samples. Clade I increased along the water column, and its proportion exceeded that of clade III in some bottom layers. Phenogeneticly, PT3 was always the dominant PT across the whole study area. From station A to C, the proportion of PT3 increased gradually, while PT2 decreased and was even absent in some water layers. RDA showed that multiple environmental variables were significantly correlated to the Synechococcus phylogenetic structure and salinity was significantly correlated to the Synechococcus phenogenetic structure. Cross-correlation co-occurrence networks revealed that Synechococcus significantly correlated with various microbial taxa, including 159 prokaryotes and 102 eukaryotes. Metabolic cycling networks showed that Synechococcus may be related to more diverse microbial metabolic cycles, including the carbon, nitrogen, hydrogen, and manganese cycles in the bottom layer of the photic zone. Furthermore, various in situ grazers of Synechococcus were predicted, including Gymnodinium clades and 46 extra Alveolata.
Due to the unseen barriers among the world's oceans, diverse genotypes of Synechococcus evolved and occupied different areas of the global ocean. To further explore the current world distribution pattern of Synechococcus, additional studies that examine the distribution of Synechococcus species and their evolution are required. Additionally, functional roles and grazer observations based on amplicon sequencing in this study tend to be preliminarily sketched, and mesocosm experiments combined with isotopic tracer and molecular omics methods should be further applied to prove the conclusions in the future.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figures S1 -S6, Tables S1-S6. Figure S1: Spearman rank correlation analysis between Synechococcus clades as determined using rpoC1 and pigment types as indirectly assessed using the cpeBA operon; Figure  S2: The relative abundances of predominant taxa in the prokaryotic community at the phylum level (a) and the eukaryotic community at the kingdom level (b); Figure S3: The relative abundance of taxa in the prokaryotic community at the genus level. (a) Community bar map showing the relative

Conclusions
The distribution of Synechococcus was studied in a typical section of the Yellow Sea. The cell abundance in the studied area was 2.4 × 10 4 cells mL −1 on average, which was significantly correlated to the sampling depth, nitrogen source (NO 2 − and NH 4 + ), and DSi. The intense decrease in the count in the lower layer probably attributes to the YSCWM. Synechococcus clade III was dominant in most samples. Clade I increased along the water column, and its proportion exceeded that of clade III in some bottom layers. Phenogeneticly, PT3 was always the dominant PT across the whole study area. From station A to C, the proportion of PT3 increased gradually, while PT2 decreased and was even absent in some water layers. RDA showed that multiple environmental variables were significantly correlated to the Synechococcus phylogenetic structure and salinity was significantly correlated to the Synechococcus phenogenetic structure. Cross-correlation co-occurrence networks revealed that Synechococcus significantly correlated with various microbial taxa, including 159 prokaryotes and 102 eukaryotes. Metabolic cycling networks showed that Synechococcus may be related to more diverse microbial metabolic cycles, including the carbon, nitrogen, hydrogen, and manganese cycles in the bottom layer of the photic zone. Furthermore, various in situ grazers of Synechococcus were predicted, including Gymnodinium clades and 46 extra Alveolata.
Due to the unseen barriers among the world's oceans, diverse genotypes of Synechococcus evolved and occupied different areas of the global ocean. To further explore the current world distribution pattern of Synechococcus, additional studies that examine the distribution of Synechococcus species and their evolution are required. Additionally, functional roles and grazer observations based on amplicon sequencing in this study tend to be preliminarily sketched, and mesocosm experiments combined with isotopic tracer and molecular omics methods should be further applied to prove the conclusions in the future.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/jmse9091018/s1, Figures S1 -S6, Tables S1-S6. Figure S1: Spearman rank correlation analysis between Synechococcus clades as determined using rpoC1 and pigment types as indirectly assessed using the cpeBA operon; Figure S2: The relative abundances of predominant taxa in the prokaryotic community at the phylum level (a) and the eukaryotic community at the kingdom level (b); Figure S3: The relative abundance of taxa in the prokaryotic community at the genus level. (a) Community bar map showing the relative abundance of the genus; The unclassified does not represent a single taxon, but the sum of all taxa not classified to genus level; (b) Chord diagram reveals the relationship between samples and Synechococcus abundance. The length of the left semicircle represents the relative abundance of the Synechococcus in the corresponding samples. The right semicircle represents the distribution proportion of Synechococcus in different samples; Figure S4: The relative abundances of OTUs belonging to phyla Cyanobacteria; Figure S5: Network of microbial community encompassing all prokaryotic and eukaryotic OTUs constructed from the correlation matrix. OTUs with Spearman coefficients > 0.6 and p < 0.05 were connected and graphed using a force-directed layout algorithm; Figure Table S1: Reference sequences of rpoC1 gene; Table S2: Reference sequences of cpeBA operon; Table S3: Spearman correlation coefficients between environmental variables; Table S4: Sequencing information and alpha diversity indices of rpoC1 and cpeBA gene; Table S5: Sequencing information and alpha diversity indices of 16S and 18S gene; Table S6: Correlative eukaryotic taxa of Synechococcus as potential in-situ grazers in the ecosystem.

Data Availability Statement:
The sequence data of this study were deposited in the sequence read archive of NCBI (US National Center for Biotechnology Information; https://www.ncbi.nlm.nih. gov/, accessed on 13 September 2021) under the accession number PRJNA595432.