Annual dynamics of eukaryotic and bacterial communities revealed by 18S and 16S rRNA metabarcoding in the coastal ecosystem of Sagami Bay, Japan

Sagami Bay, Japan is influenced by both the warm Kuroshio Current and the cold Oyashio Current and rich nutrients are supplied from multiple river sources and the deep-sea, forming a dynamic ecosystem. The aim of the present study was to investigate eukaryotic and bacterial communities in the coastal waters of Sagami Bay, using 16S rRNA and 18S rRNA sequencing and to assess the seasonal and vertical dynamics in relation to physicochemical and biological conditions. Eukaryotic and bacterial communities showed synchronous seasonal and vertical changes along with environmental variability. Diversity of plankton community suspended in the surface was lower than those at the subsurface layers in both the eukaryotes and bacteria communities; however, community diversity showed different characteristics in the subsurface where the eukaryotic community at the deeper layer (100 m) was as low as the surface and highest in intermediate depth layers (10–50 m), while that of bacterial community was highest in the deeper layer (100 m). The annual variability of the coastal microbial communities was driven, not only by the seasonal changes of abiotic and biotic factors and short-term rapid changes by river water inflow and phytoplankton blooms, but also largely influenced by deep-seawater upwellings due to the unique seafloor topography.


Introduction
Marine pelagic ecosystems are interconnected by two types of food chains: grazing food chains and microbial food chains (microbial loop). Classical grazing food chains begin with phytoplankton photosynthesis (producers), which become food for zooplankton, such as herbivorous copepods and further preyed upon by larger organisms, such as fishes (Fenchel 1988;Azam 1998). The microbial loop is a relatively new food chain concept beginning with bacteria introducing dissolved organic material (DOM) to the classical food chain via heterotrophic protist and microzooplankton (Pomeroy 1974;Azam et al. 1983;Pomeroy et al. 2007). In the marine ecosystem, physical factors, such as currents, temperature, density and light, chemical factors, such as nutrients and oxygen and biological factors related directly or indirectly with marine organisms, regulate the system. Several studies have shown that diversity, abundance and the community structure of marine plankton correspond with changes in physical and chemical factors and interaction between marine plankton (Sin et al. 1999;Hays et al. 2005;Sogawa et. al. 2013; and bacteria influence one another, having significant impact on biogeochemical cycles and form the basis of marine ecosystem community and diversity Little et al. 2008;Amin et al. 2015). The dynamics of coastal marine microbial communities are also affected by seasonal change in abiotic and biotic factors and short-term effects of river water and phytoplankton blooms (Lucas et al. 2015). Thus, a comprehensive understanding of plankton communities in marine ecosystems and revealing changing mechanisms are essential for sustainable use and effective management of marine resources (Ministry of Environment of Japan 2011; Sogawa et al. 2018). However, it is difficult to grasp the diverse plankton community by morphological identification using microscopy as many of the marine plankton are small and amorphous (Herndl and Peduzzi 1988;Bochdansky and Herndl 1992;Hirai et al. 2017) and require advanced techniques and experience for identification (Abad et al. 2016(Abad et al. , 2017. The advances of molecular biological techniques, such as DNA metabarcoding using high-throughput sequencing, have become an efficient and powerful tool for monitoring and assessing marine communities (Valentini et al. 2016;Ruppert et al. 2019;Suter et al. 2020).
Sagami Bay, located on the Pacific Ocean side of middle Honshu Island, Japan, is known for its variety of aquatic landform, containing rivers (terrestrial runoff), coastal zones and deep seas (deeper than 1,000 m). Abundant nutrients are supplied from multiple river sources and deep ocean upwelling, forming rich harvest grounds containing a wide variety of migratory and local fishery catch. The Bay is influenced by both warm Kuroshio and cold Oyashio currents: Kuroshio, the western boundary current of the North Pacific subtropical gyre, transports warm and saline subtropical waters as well as a diverse assortment of organisms and Oyashio, the western boundary current of the western subarctic gyre, branches and flows from the opposite direction of Kuroshio in the mesopelagic layer and flows into Sagami Bay. Consequently, the Bay would be an appropriate location for monitoring diverse plankton communities and survey of the area could provide indispensable information towards understanding the coastal aquatic ecosystem. Previous studies reported temporal/seasonal variation of bacterio-, phyto-and zooplankton in Sagami Bay and their ecology in relation to environmental factors. (Takasuka et al. 2003;Ara et al. 2006;Miyaguchi et al. 2006;Shimode et al. 2006;Hashihama et al. 2008;Baek et al. 2009;Baki et al. 2009;Tsuchiya et al. 2013Tsuchiya et al. , 2016Sugai et al. 2016Sugai et al. , 2018Sugai et al. , 2020. A few studies analysed diversity of plankton using 18S and 16S rRNA gene (Kita-Tsukamoto et al. 2006;Kok et al. 2012a, b). However, no studies have conducted comprehensive seasonal analysis on plankton communities in Sagami Bay that address environmental change and interactions between organisms. In this study, a continuous survey over a year was conducted at a longterm monitoring site of Sagami Bay and 18S-RNA and 16S-rRNA genes were analysed to reveal the annual vertical dynamics of eukaryotic and bacterial communities in the temperate coastal pelagic ecosystem.

Field sampling and DNA extraction
Seasonal observations were carried out at the shelf station (120 m depths, latitude 35°09.75'N, 139°10.00'E longitude), located in the western part of Sagami Bay (Fig.  1). Continuous surveys were conducted 23 times over a one-year period from October 2018 to February 2020 (Table 1) on board R/V "Tachibana" of the Manazuru Marine Center for Environmental Research and Education (MMCER), Graduate School of Environment and Information Sciences, Yokohama National University, Yokohama National University.
Vertical profiles of temperature, salinity, density, chlorophyll a fluorescence and dissolved oxygen were measured continuously using a RINKO-Profiler (ASTD102, JFE Advantech Co., Ltd). Mixed layer depths (MLDs) were calculated from density with the following equation: MLD: Δσt(z) = σt(z) -σt(5m) > 0.125 kg/m3 (Monterey and Levitus 1997). The euphotic zone (1% light level; EZ1%) was measured by Profiling Reflectance Radiometer (PRR) / Compact Optical Profiling System (C-OPS). Seawater samples were collected at several depth layers (maximum depth 100 m) in Niskin bottles and directly from the surface layer in sampling bottles. Seawater samples of 1litre for 16S and 18S rRNA analysis were filtered through Sterivex filter (0.22 μm) and DNA was extracted using 5% Chelex buffer (Nagai et al. 2012;Tanabe et al. 2016). Sampling bottles and filtering equipment were washed with a 10% commercial bleach solution and rinsed with distilled water before usage. Seawater samples for discrete vertical profiles of nutrient (nitrate + nitrite, silicate, phosphate) concentrations were collected at the same depths and measured using a continuous flow injection analyser (QuAAtro, BL-Tech).

Paired-end library preparation, metabarcoding sequencing
For detection of eukaryotic and bacterial communities, respectively primer pairs of V7 -9 region in 18S-rRNA gene (18S-V7F: TGGAGYGATHTGTCTGGTTDATTCCG and 18S-V9R: TCACCTACGGAWACCTTGTTACG; Sildever et al. 2019) and V3-4 region in 16S-rRNA gene (16S-341F: CCTACGGGNGGCWGCAG and 16S-805R: GACTACHVGGGTATCTAATCC) (V3-V4 region, Sinclair et al. 2015) were used as universal primers for metabarcoding analysis. Two-step PCR for the construction of paired-end libraries and HTS on Illumina Miseq 300 PE platform (Illumina, San Diego, CA, USA) followed the protocol described by Dzhembekova et al. (2017). The first PCR was performed using a thermal cycler (PC-808; ASTEC, Fukuoka, Japan) in a 25 μl reaction mixture containing 1.0 μl template DNA (< 1 ng), 0.2 mM of each dNTP, 1× PCR buffer, 1.5 mM Mg 2+ , 1.0 U KOD-Plus-ver.2 (Toyobo, Osaka, Japan) and 1.0 μM of each primer. PCR conditions were first denaturation at 94 °C for 3 min, followed by 30 -32 cycles (18S-rRNA gene) or 28 cycles (16S-rRNA gene) at 94 °C for 15 s, 56 °C for 30 s and 68 °C for 40 s. PCR amplification was verified by 1.5% agarose gel electrophoresis. PCR products were purified using Agencourt AMPure XP (Beckman Coulter, Brea, CA, USA) and eluted in 25 μl TE buffer according to the manufacturer's protocol. The second PCR was performed using the same conditions for the first PCR, except the reaction mixture volume was 50 μl with the addition of 2.0 μl diluted PCR product from the first round as a template. The PCR cycling conditions were as follows: initial denaturation at 94 °C for 3 min, followed by 10-12 cycles at 94 °C for 15 s, 56 °C for 30 s and 68 °C for 40 s. PCR amplification was again verified by agarose gel electrophoresis and the PCR products were purified using an Agencourt AMPure XP (Beckman Coulter, Brea, CA, USA). The amplified PCR products were quantified, mixed equally and stored at -30 °C until used for sequencing.

Treatment processes and operational taxonomic unit picking
Nucleotide sequences were demultiplexed depending on the 5′-multiplex identifier tag and primer sequences according to the default format in MiSeq. Sequences longer than 300 bp were truncated to 300 bp by trimming the 30 tails. The trimmed sequences shorter than 250 bp were filtered out. Demultiplexing and trimming were performed using Trimmomatic version 0.35. The remaining sequences were merged into paired reads using Usearch version 8.0.1517. In addition, singletons were removed. Sequences were then aligned using Clustal Omega version 1.2.0. Multiple sequences were aligned with each other and only sequences that were contained in more than 75% of the read positions were extracted. Filtering and part of the multiple alignment process were performed as described in the Miseq standard operating procedure using Mothur (Schloss et al. 2011). Erroneous and chimeric sequences were detected and removed using Mothur (Edgar et al. 2011). Demultiplexed and filtered, but untrimmed sequence data are deposited into the DDBJ Sequence Read Archive (access no. DRA013426). The sequences were clustered to select representative sequences at a 0.99 level of sequence identity by CD-HIT-EST version 4.6.8 (Li and Godzik 2006).

Taxonomic identification of the OTUs
The sequence database used to identify OTUs was downloaded from GenBank on 20 August 2020. The taxonomic identification of each OTU was performed by the BLAST search (Cheung et al. 2010) with NCBI BLAST+ 2.10.1+ (Camacho et al. 2009). The taxonomic information was obtained from the BLAST hit with the top bitscores for each query sequence. The OTUs of the same top-hit were merged and the OTUs were also merged when an OTU matched several data with the same bitscore and top-hit similarity. Several data, therefore, were sometimes merged at an OTU. The removal of sequences containing errors was imperfect after the successive processes of MPS data treatment; sequences containing different types of errors derived from the original ones remained in the following analytical steps. Therefore, these sequences were detected as artificial OTUs with the same top BLAST hit name,

Statistical analysis
All multivariate analyses of eukaryotic community structure and diversity were performed using PRIMER version 7 with PERMANOVA+ add-on software (Anderson et al. 2008;Clarke and Gorley 2015). For multivariate analyses, Bray-Curtis similarities (Bray and Curtis 1957) were calculated amongst samples, based on the log-transformed abundance data of sequence reads. The OTU hit as Chloroplast was excluded from the 16S analysis. A non-metric multidimensional scaling (NMDS) ordination was used to visualise the differences in the community amongst samples and analysis of similarity (ANOSIM) tested differences amongst water depths and months (Clarke 1993). A distance-based linear model (DISTLM) and distance-based redundancy analysis (DBRDA) were conducted to explain changes in eukaryotic and bacterial community with respect to environmental variables. Environmental data were normalised and square root transformed prior to the analysis. Multicollinearity was tested prior to DISTLM and DBRDA and variables with correlations stronger than 0.95 were removed. Plankton diversity was measured according to the total number of OTUs together with the Shannon-Wiener Diversity Index (Shannon and Weaver 1949) and the Evenness (Pielou 1966).

Physical and chemical oceanographic environment
Vertical profiles of temperature, salinity, density, chlorophyll fluorescence, dissolved oxygen and nutrients in Sagami Bay during the observation period are shown in Fig. 2 Diversity analysis was conducted for eukaryotic community amongst depths (Fig. 3). The mean numbers of OTUs were lower in the surface (175)  The number of OTUs were higher in intermediate depths layers (10-50 m) and lower in the surface and bottom depth layer, although rarefaction curves of all depths did not reach asymptote (Fig. 4).
Diversity analysis was also conducted for bacterial community amongst depths (Fig. 3). Characteristics of the mean numbers of OTUs were contrary to the eukaryotic community that were higher in the surface (321) and the bottom depth layer (100 m; 318) and lower in the intermediate depth layer (30 m and 50 m; 282 and 289, respectively). Pielou's Evenness was the lowest in the surface (J′ = 0.61), the same as the eukaryotic community; however, they gradually increased with increasing depth and reached the highest in the bottom depth layer (100 m; J′ = 74). The Shannon-Wiener Diversity Index was also lowest in the surface (H′ = 3.5) and highest in the bottom depth layers (100 m; H′ = 4.2). Rarefaction curves of all  (Fig. 4). Diversity of plankton community suspended in the surface was lower than those at the subsurface layers in both the eukaryotes and bacteria communities; however, the diversity of the eukaryotic community at the bottom layer (100 m) was as low as the surface, while that of the bacterial community was highest, contrary to the surface.

Annual dynamics of eukaryotic and bacterial community
The ANOSIM test revealed significant differences amongst months for both eukaryotic and bacterial community (global test: R = 0.395, P = 0.001 for 18S and R = 0.701, P = 0.001 for 16S). On the contrary, there were no significant differences amongst depth layers for both the eukaryotic and bacterial communities (ANOSIM global test: R = 0.009, P = 0.389 for 18S and R = 0.005, P = 0.470 for 16S). NMDS ordination represented samples collected in similar seasons (months) and depth layers are plotted closer, particularly during winter (Fig. 5). In winter, plankton communities were more vertically homogeneous as samples are plotted closer compared to spring to autumn. Distance-based multivariate linear model (DISTLM) analysis revealed significant relationships (p < 0.01) between eight out of nine variables and eukaryotic and bacterial community structures, respectively (Table 2). One variable, phosphate, was excluded prior to analysis as it was strongly correlated (|r| > 0.95) with nitrate + nitrite. MLDs accounted for the largest proportion of fitted variance for eukaryotic community, followed by nitrate + nitrite, density, chlorophyll a, EZ1% and temperature (11.5-4.2%) ( Table 2). Axis 1 of the distance-based redundancy analysis (DBRDA) accounted for 46.3% of the fitted model (11.9% of the total variation) and strongly correlated with MLDs (Fig. 6). Axis 2 of DBRDA accounted for 28.6% of the fitted model (7.3% of the total variation) and strongly correlated with nitrate + nitrite.
In the top 10 OTU of eukaryotic communities, Amobophyra (Alveolata, Dinophyceae) and Astrosphaera (Rhizaria, Polycystinea), which were abundant at 100 m depths, were positively related to nitrate + nitrite and negatively related to chlorophyll a (Fig. 6). On the contrary, Paracalanus (Opisthokonta, Metazoa), which were abundant on the surface (0 m), was positively correlated to chlorophyll a. Karlodinium veneficum, Heterocapsa rotundata (Alveolata, Dinophyceae) and Leptocylindrus convexus (Stramenopiles, Bacillariophyta), which were abundant in spring-autumn, were negatively correlated to MLDs. On the other hand, Ostreococcus, Bathycoccus prasinos and Micromonas pusilla (Viridiplantae, Chlorophyta), which were abundant in winter, were positively correlated to MLDs. MLDs also accounted for the largest proportion of fitted variance for bacterial community, followed by density, EZ1%, nitrate + nitrite and chlorophyll a (11.6-8.7%) ( Table 2). Axis 1 of the DBRDA accounted for 33.6% of the fitted model (14.3% of the total variation) and strongly correlated with MLDs, EZ1%, density and chlorophyll a (Fig. 6). Axis 2 of DBRDA accounted for 24.7% of the fitted model (10.5% of the total variation) and strongly correlated with density and MLDs. These two axes were strongly related to water mass structure and seasonality.
In the top 10 OTU of bacterial communities, Flavobacteriales, Cellvibrionales and Rhodobacterales, which were abundant in spring and summer, were positively related to chlorophyll a and negatively related to MLDs (Fig. 6). Synechococcus (CC9902) was located in the second quadrant in summer and autumn and positively correlated to temperature and negatively correlated to density. SAR11 subclade IV was abundant in autumn, where chlorophyll a was relatively low and negatively related to Flavobacteriales and Rhodobacterales. SAR11 subclade Ia was negatively correlated to Flavobacteriales and Rhodobacterales, as well as SAR11 subclade IV.

Cluster analysis of eukaryotic (phytoplankton) community
The eukaryotic community was classified into seven clusters, which represent clear seasonal and vertical patterns (Fig. 10, Suppl. material 1: Fig. S1). Cluster 7 was further classified into four clusters (7-1, 7-2, 7-3 and 7-4). Cluster 1, 3 and 4 composed of samples from 100 m depth of October, Spring (April-May) and February, respectively. Cluster 2 composed mostly from surface samples in late spring to autumn (May-September). Cluster 5 composed of samples from all vertical depths in November and the subsurface layer of different months. Cluster 6 composed of samples in all vertical depths of late spring to summer (May -July). Cluster 7 composed of samples in all vertical depths of winter to early spring (December-April). The taxonomic groups of phytoplankton (diatom, dinoflagellates, haptophytes and green algae) were represented by heatmap at class level (Fig. 11). Noctilucales (dinoflagellates) showed a sudden increase in the surface in May and July. The abundance of reads in Peridinales (dinoflagellates) increased on the surface in May and deeper than 30 m in November, whereas Gymnodiniales (dinoflagellates) increased in the shallower depths (less than 10 m) in July and November to February. Syndiniales (dinoflagellates) showed large abundance in sequence reads in the deeper depths (deeper than 30 m with the largest at 100 m depth) of spring to autumn. Amongst diatoms, Coscinodiscophyceae showed large abundance in sequence reads from spring to autumn in the shallower depths near the surface (≥ 10 m) and especially large in May from the surface to 100 m depth. Mediophyceae (diatoms) showed a sudden increase on the surface in April. Amongst green algae, Mamiellophyceae showed large abundance in reads during winter in all depths and in April with peak in 30 m depth. Amongst haptophytes, Prymnesiales showed the largest abundance in reads in April with a peak in the surface (≥ 30 m) and during winter at all depths.

Eukaryotic community
Astrosphaera hexagonalis/ Arachnosphaera myriacantha (radiolarian, Rhizaria) dominated in cluster 1 (October 2018 at 100 m depth) of eukaryotic community (75.3%). These species form a monophyletic group by molecular phylogenetic analysis using 18S rRNA gene (Yuasa et al. 2009). They were recorded in the warm surface water in the Pacific  Ocean including the water area around the Kuroshio Current and surface sediment (Suzuki and Sugiyama 2001;Yuasa et al. 2009;Boltovskoy et al. 2010). The surface layer of Sagami Bay is strongly influenced by the oligotrophic Kuroshio (Kawabe and Yoneno 1987;Iwata and Matsuyama 1989;Hashihama et al. 2008;Baek et al. 2009). Furthermore, the same species was also most abundant in sequence reads amongst Rhizaria in the waters around the Kuroshio current south of Japan, with its distribution peak at 500-2000 m depths (Sogawa et al., submitted manuscript, 2022). These species were located in the direction of high nitrate + nitrite concentration, based on the DBRDA results, suggesting their coupled distribution to high nutrient deep-sea water.
The genus Amoebophrya is known to parasitise on a wide range of other dinoflagellates (Gunderson et al. 2002;Park et al. 2004). It also parasitises red tideforming species (Taylor 1968) and is detrimental to harmful bloom algae (Miller et al. 2012). Parasitism was reported to occur during the bloom events in the surface or sub-surface layers (Nishitani et al. 2021). However, in our study, the parasitic syndinian (order Sindiniales) Amoebophrya sp. (dinoflagellates) dominated in cluster 4 (February at 100 m depth) of the eukaryotic community. No significant difference was observed between 100 m depths in October 2018 and February 2020 in the cluster analysis of only dinoflagellates. Amoebophrya sp. showed positive correlation with three dinoflagellates species; another parasitic, but non-toxic dinoflagellate Pentapharsodinium tyrrhenicum that forms resting cysts (r = 0.52), Euduboscquella sp. (r = 0.55) and unarmoured dinoflagellates Takayama cf. pulchellum (synonym Gymnodinium pulchellum) that causes a red tide bloom which may lead to fish kills (r = 0.48). Many species of dinoflagellates form resting cysts and have different life cycles that consist of a cyst and motile stage. The genus Euduboscquella parasitises Tintinnids and produces cysts (Dolan et al. 2013;Shang et al. 2019;Liu et al. 2020). The temporal dominance of Amoebophrya sp. could be considered caused by turbulence from bottom sediment from sudden upwelling events.
In the present study area, early summer is the rainy season and direct typhoon passes occasionally occur during the summer. High precipitation and abundant nutrients supplied by terrestrial run-off cause an increase in chlorophyll a concentration a few days later (Baek et al. 2009;Tsuchiya et al. 2015). Subsequently, abundance of harmful algal bloom (HAB) species of dinoflagellates (Tripos furca, previously described as Ceratium furca, Tripos fusus, previously described as Ceratium fusus and Noctiluca scintillans) and diatoms (Rhizosolenia delicatula, Hemiaulus sinensis and Navicula spp.) increase and exhaust nitrogen and phosphorus supply, whereas silicate concentrations are not limited (Miyaguchi et al. 2006;Baek et al. 2008Baek et al. , 2009. It can be inferred that the same succession occurred in our study period, based on nutrient and salinity results. Dinoflagellate species are common bloom-forming species in coastal waters; however, our study showed increased abundances for N. scintillans, but not for T. furca and T. fusus. Furthermore, no recent blooms of the two Tripos species occurred at the monitoring site (Shimode's periodical observation, data not shown). On the other hand, small mixotrophic dinoflagellates Heterocapsa rotundata (3-10 μm) and Karlodinium veneficum dominated. The former species, widely distributed in tropical and temperate

Jul_0m
Jul_0m Dinoflagellates-Phytodiniales coastal areas, is known to form large blooms in temperate estuaries (Millette et al. 2017). The latter species is ichthyotoxic and is distributed worldwide causing massive fish kills in coastal areas (Müller et al. 2019). The Kuroshio Current formed a large meander, which began in August 2017, continuing for a long period and was also observed during the survey periods (Japan Meteorological Agency, 20 May 2021). The surface layer of Sagami Bay is strongly influenced by the oligotrophic Kuroshio as it changes inflow to the Bay, based on the Kuroshio path (Kawabe and Yoneno 1987;Iwata and Matsuyama 1989;Hashihama et al. 2008;Baek et al. 2009). Therefore, the long-term large meandering might have caused changes in the plankton community and occurrence of bloom-forming species in Sagami Bay.

Bacterial community
The bacterial communities transitioned seasonally in the surface waters: from autumn-winter in 2018 (cluster 3-1) to winter-spring (cluster 4-4), rainy season before summer (cluster 4-2) and summer-autumn (cluster 2) in 2019 and to winter in 2019-2020 (cluster 3-3). The succession could be derived from the change of dominant taxa. Synechococcales and SAR11_clade were dominant in cluster 3-1, Flavobacteriales and Rhodobacterales became dominant in clusters 4-4 and 4-2. In cluster 2, Flabovacteriales, Rhodobacteriales, SAR11_clade, Synechococcales and cellvibrionales were evenly dominant. After that, SAR11_ clade and Flavobacteriales became dominant. Flavobacteriales and Rhodobacterales are well known as copiotrophic taxa and respond quickly to the flux of fresh labile DOM derived from phytoplankton, subsequently outcompeting other community members as observed during phytoplankton blooms (Buchan et al. 2014;Luria et al. 2017).
In this study, samples clustered in 2, 4-2 and 4-4 were located in the direction of high chlorophyll a concentration, based on the DBRDA results, suggesting Flavobacteriales and Rhodobacterales may have responded quickly to the DOM supply released from phytoplankton. SAR11 was negatively correlated to Flavobacteriales and Rhodobacterales in the present study. The inverse correlation between those taxa was also observed in the coastal Southern Ocean (Williams et al. 2013). Flavobacteriia is essential for the breakdown of organic matter and they make labile organic compounds available to Proteobacteria, including SAR11 (Williams et al. 2013). SAR11 cells oxidise many labile, low-molecular-weight organic compounds and their small cells with periplasms and many substrate-biding proteins (SBPs) are likely to be superior competitors for low-molecular-weight DOM in oligotrophic ocean environments (Giovannoni 2017). Indeed, dissolved organic carbon (DOC) concentration increases during spring to summer and decreases during winter in the mixed layer of Sagami Bay (Sugai et al. 2016), suggesting that the DOC dynamics and the succession from Flavobacteriales and Rhodobacterales to SAR11 are consistent. Their different ecological roles could separate their occurrences spatially and seasonally. SAR11 can be divided into several subclades, each of which has a different response to the environment (Giovannoni 2017;Bolanos et al. 2021). In this study, three SAR11 subclades (Ia, II and IV) were dominant. Subclade Ia was found with relatively regular contributions in all clusters, whereas subclade IV sporadically increased in clusters 1, 3-1, 3-2 and 3-3. Subclade Ia (Ia.1 in cold environment and Ia.3 in temperate and tropical environments) was characterised by surface ecotypes (Giovannoni 2017) and dominated in subpolar regions in the Atlantic Ocean (Bolanos et al. 2021) and poly-and euhaline sites in San Francisco Bay (Rasmussen et al. 2021), which are consistent with the present study. Since there was no distinct spatiotemporal variation in subclade Ia, DBRDA analysis did not show a relationship between environmental factors and subclade Ia occurrence. At BATS site, subclade IV was found at lower abundance in the surface layer, predominantly around 80 m depth (Vergin et al. 2013). Besides, subclade IV was varying only slightly and not dominant throughout the seasons in the North Atlantic (Bolanos et al. 2021). Subclade II was relatively low (0.89 ± 0.26%) in clusters 4-2, 4-3 and 4-4 (mainly in spring) compared to that in the other clusters (2.2 ± 0.4%). There are seasonal blooms in subclade IIa on the surface and IIb in the mesopelagic during spring (Carlson et al. 2009;Carlson et al. 2010;Vergin et al. 2013), which disagreed with the results of the present study. The function of SAR11 subclades can differ depending on the ocean area and environments, although there is a presumed relationship between ecotype classification and taxonomy in SAR11 subclades (Giovannoni 2017).
Cyanobacteria were relatively abundant in clusters 2 and 3-1, accounting for up to 61% and DBRDA ordination revealed that cyanobacteria (Synechococcus_ CC9902) were highly correlated with water temperature. The sporadic increase of cyanobacteria from summer to autumn is consistent with previous studies in the same Sagami Bay area (Mitbavkar et al. 2009;Ara et al. 2019). Cyanobacterial abundance and growth rate are dependent principally on the temperature (Chen and Liu 2010;Chen et al. 2014;Wang et al. 2020). In Sagami Bay, Synechococcus biomass (216 ± 43 mg C m -2 on average annually) was much higher than that of Prochlorococcus (11.5 ± 1.7 mg C m -2 on average annually) (Mitbavkar et al. 2009). This is consistent with the present results that the ratio of Synechococcus read numbers to those of Cyanobacteria was 97 ± 7%. These two cyanobacterial genera show different ecological distribution: Synechococcus increases in nutrient-rich waters (generally coastal waters), whereas Prochlorococcus is most abundant in offshore nutrient-depleted areas (Partensky et al. 1999a). The global distribution prediction for Prochlorococcus and Synechococcus indicate that the abundance of Prochlorococcus and Synechococcus will increase overall by the end of the 21 st century under the RCP4.5 scenario; however, the community structure will change significantly in the mid-latitudes (around 35°N or 35°S) and the abundance of Synechococcus may decrease (Flombaum et al. 2013). Flombaum et al. (2013) found that temperature is the primary control of the biogeography of both taxa and nutrient concentration had a limited influence on it. However, Prochlorococcus has been reported to be in a competitive disadvantage in elevated nutrient environments (Partensky et al. 1999b) and its future prediction in coastal waters is still debatable. Therefore, it is essential to understand the dynamics of both Prochlorococcus and Synechococcus in the coastal areas to predict the dynamics of primary production in the future since marine cyanobacteria could be responsible for ~ 25% of ocean net primary productivity (Flombaum et al. 2013).
SAR324_clade increased at deep waters (100 m depth) in summer and autumn, especially in October (5.4%), November (3.3%), December (6.1%) 2018 and July 2019 (5.8%). The results suggested that SAR324 gradually increased in their relative abundance several months after the MLD exceeded 100 m and vertical mixing occurs to the deeper layers. The seasonal variation agreed with a previous study reporting that SAR324 increased in the upper mesopelagic in autumn in the North Atlantic (Bolanos et al. 2021). SAR324 has been associated with a chemolithotrophic lifestyle and their distribution was well documented in oxygen-deficient waters (Zaikova et al. 2010), hydrothermal plumes (Sheik et al. 2015) and aphotic layers (160-3000 m, Wright et al. 1997;López-García et al. 2001;Bolanos et al. 2021). Draft genome analysis has revealed that SAR324 is also involved in new metabolic features, such as aromatic compound degradation and C1 metabolism, such as methanotrophy (Cao et al. 2016). This group might have contributed to carbon cycling in the stratified, deep layers in temperate coastal waters.

Overall plankton community and diversity dynamics
The present study revealed clear seasonal and vertical variations in the eukaryotic and bacterial communities consistent with changes in the environmental parameters. Eukaryotic and bacterial communities showed synchronous seasonal and vertical changes corresponding to the physical change; summer when the thermocline develops, winter when the mixed layer depth (MLD) is deep and transition periods that MLD gradually shallows (spring) and gradually deepens (autumn). In particular, DBRDA results suggested that the succession of eukaryotic and bacterial communities is closely related to the variations of the MLD, followed by nutrients for the eukaryotic community and temperature and chlorophyll a for the bacterial community.
The planktonic community was also vertically uniform in winter when vertical mixing formed a uniform water column and abiotic environment, although there was a slight time lag between the eukaryotic community (November-February) and the bacterial community (December-April). Productivity and biomass of bacterial communities are known to depend on eukaryote-derived resources (Williams et al. 2013;Buchan et al. 2014;Tsuchiya et al. 2015;Sugai et al. 2016;Luria et al. 2017;Tsuchiya et al. 2019); therefore, the bacterial commu-nity would also be formed dependent on the eukaryotic community. Significantly different clusters formed during different seasons at 100 m depths (clusters 1, 3 and 4 for eukaryotic community and clusters 1 and 4-3 for bacterial community). Upwelling of deep seawater supplies rich nutrients to the coastal areas from Sagami Trough (more than 1000 m deep), located in the centre of the Bay (Takahashi et al. 1980;Taira and Teramoto 1985;Kawabe and Yoneo 1987). In our study, deep-waters, rich in nutrients, were supplied to deeper layers (100 m) when vertical mixing weakened. It is interesting to note that the different eukaryotic communities formed during each upwelling event (cluster 1 in October 2018, cluster 3 in April-May 2019 and cluster 4 in February 2020). The annual dynamics of the plankton and bacterial communities are driven not only by the seasonal changes of abiotic and biotic factors and short-term rapid changes by river water inflow and phytoplankton blooms, but also largely influenced by deep-seawater upwelling due to the unique seafloor topography.
Vertical characteristics of community diversity were different between the eukaryotic and bacterial communities; the lowest diversity was observed in the surface water in both communities, whereas the highest diversity was observed in intermediate depth layers (10-50 m) for the eukaryotic community and in the deeper layer (100 m) for the bacterial community. The number of OTUs and the Evenness were both low at the surface and 100 m layer in the eukaryotic community. The higher OTU number in the intermediate depth layers could be explained by seasonally changing MLD; the subsurface chlorophyll maxima (SCM), formed at the optimal depth where light availability and upward nutrient flux are sufficient to sustain phytoplankton growth in the water column, is characterised by a high species diversity of phytoplankton (Furuya and Marumo 1983). It is also well known that ultraviolet radiation is harmful to many organisms and the surface is illuminated most with ultraviolet light . The Evenness could be very low in the surface and 100 m layers due to phytoplankton and dinoflagellate blooms caused by shortterm rapid environmental change by river water inflow and upwelling events as noted above. In the bacterial community, the number of OTUs does not vary much amongst depths, but the Evenness gradually increased with increasing depth and reached the highest in the 100 m layer along with the Diversity Index. This could be explained by low oxygen concentration in the deeper layers (except temporal vertical mixing in winter) and higher concentration near the surface as reported in the other water areas (Stevens et al. 2008;Signori et al. 2014).
network of environment and fisheries information) by the Ministry of Agriculture, Forestry and Fisheries of Japan. We would like to thank all members of Team Manazuru in Yokohama National University and Soka University for observation support and members of Fisheries Research Institute, Japan Fisheries Research and Education Agency for measurement support. Ia, II and IV shown by their relative reads abundances in total bacterial community. Figure S5. Cluster analysis of diatom and relative reads abundance of top 15 OTUs classified to species level. Figure S6. Cluster analysis of dinoflagellates and relative reads abundance of top 15 OTUs classified to species level. Figure S7. Cluster analysis of green algae and relative reads abundance of top 15 OTUs classified to species level. Figure S8. Cluster analysis of haptophytes and relative reads abundance of top 15 OTUs classified to species level. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/mbmg.6.78181.suppl1