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.

Figure 1.

Schematic illustration of surface current around Japan (left) and location of survey site (right, red circle) in Sagami Bay, Japan.

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).

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²⁺, 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.

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).

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, but with slight differences. To avoid overestimation of the OTUs, these artificial OTUs were merged into a single OTU with the greatest similarity score. The phylum and class/order (prokaryotes) level taxonomic grouping were also confirmed by the feature-classifier command in Qiime2 (https://docs.qiime2.org/2020.8/) (Bokulich et al. 2018) using the Silva database (https://docs.qiime2.org/2020.8/data-resources/common/silva-138-99-nb-classifier.qza) (Yilmaz et al. 2014) as the taxonomic database.

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).

Vertical profiles of temperature, salinity, density, chlorophyll fluorescence, dissolved oxygen and nutrients in Sagami Bay during the observation period are shown in Fig. 2. Water temperature in the sampling layer ranged from 26.9 °C at the surface in August 2019 to 13.9 °C at 100 m depth layer in April 2019. Salinity varied from 34.6 to 32.8 and the highest value was observed at 100 m depth layer in February 2019 and the lowest value was observed at the surface in July 2019. Density was highest at 100 m depth layer in April 2019 (25.8 σ_t) and lowest at the surface in July (21.2 σ_t). The mixed layer depths (MLDs) were the deepest in February when vertical mixing was strong and < 20 m in May to November when the seasonal thermocline was well developed. Chlorophyll a varied from 2.24 μgl^-1 to 0.11 μgl^-1 with highest at 10 m depth layer in April 2019 and the lowest at 100 m depth layer in October 2018. Nutrient concentration ranged from 14.53 μM to 0.04 μM for nitrate + nitrite and 1.012 μM to 0.054 μM for phosphate with the maximum value observed at 100 m depth layer in April 2019 and the minimum value at 10 m depth in May 2019 for nitrate + nitrite and in July for phosphate. The maximum value of silicate was observed at the surface in July 2019 (22.60 μM) and the minimum value at the surface in September 2019 (0.98 μM).

Figure 2.

Vertical profiles of temperature, salinity, density, chlorophyll fluorescence, dissolved oxygen (DO) and nutrients in Sagami Bay during observation period. Crosses show water sampling layer. The Mixed Layer Depths (MLDs), calculated

In the study area, 1,660 OTUs (985 OTUs from the surface) of eukaryotic plankton were detected. Supergroup Alveolata was the largest of 658 OTUs (377 OTUs from the surface), followed by Stramenopiles of 368 OTUs (259 OTUs), Opisthokonta of 262 OTUs (131 OTUs), Rhizaria of 160 OTUs (77 OTUs), Hacrobia of 90 OTUs (66 OTUs) and Viridiplantae of 77 OTUs (50 OTUs). The number of OTUs detected in the study area were ca. three times larger in the bacterial community (4477 OTUs and 2838 OTUs from the surface). Phylum Proteobacteria was the largest of 2483 OTUs (1573 OTUs from the surface), followed by Bacteroidota of 1067 (786 OTUs), Verrucomicrobiota of 112 (66 OTUs), Actinobacteriota of 87 (70 OTUs), Bdellovibrionota of 78 OTUs (31 OTUs) and Cyanobacteria of 71 (56 OTUs).

Diversity analysis was conducted for eukaryotic community amongst depths (Fig. 3). The mean numbers of OTUs were lower in the surface (175) and the bottom depth layer (100 m; 169) and highest in the intermediate depth layer (30 m; 211). Pielou’s Evenness was the lowest in the surface (J′ = 0.64) and higher in the intermediate depth layers (10–50 m; J′ = 0.69–0.71). The Shannon-Wiener Diversity Index was also lower in the surface (H′ = 3.3) and the bottom layer (H′ = 3.4) and higher in the intermediate depth layers (10–50 m; H′ = ca. 3.7). 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).

Figure 3.

Species richness, Shannon-Wiener Diversity (H′) and Pielou’s Evenness (J′) of eukaryotic and bacterial communities in Sagami Bay comparing amongst sampling depth layers.

Figure 4.

Sample-based rarefaction curves of each sampling depth layer of eukaryotic and bacterial community.

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 depths have not reached asymptote, but were higher in the surface and lower in the intermediate depth layers (30 m and 50 m), contrary to the eukaryotic community (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.

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.

Figure 5.

Non-metric multidimensional scaling (NMDS) ordination of eukaryotic (upper figures) and bacterial communities (lower figures) in Sagami Bay. Coloured area showed plot ranges of each depth layer (left figures) and each month (right figures). Numbers in parenthesis after month represent year of sampling; “18” as 2018, “19” as 2019 and “20” as 2020.

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.

Figure 6.

Distance-based redundancy analysis (DBRDA) ordination of eukaryotic and bacterial community. The distance-based redundancy analysis was constrained by the best-fit explanatory variables from a distance-based multivariate linear model (DISTLM).

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.

Bacterial communities were classified into four clusters at 30% of the Bray-Curtis Similarity in cluster analysis (Fig. 7, Suppl. material 1: Fig. S1). Moreover, clusters of 3 and 4 were further classified into 3 (3-1, 3-2 and 3-3) and 4 groups (4-1, 4-2, 4-3 and 4-4) at 39% and 41% of the Bray-Curtis Similarity, respectively (Fig. 7, Suppl. material 1: Fig. S1). In cluster 1, limited in deep samples, SAR11_clade was most dominant (14 ± 8%), followed by Flavobacteriales (8.7 ± 2.4%), Sphingomonadales (6.9 ± 6.2%) and Oceanospirillales (6.2 ± 3.2%) (Fig. 8, Suppl. material 1: Fig. S2). SAR324_clade and SAR406_clade were also relatively abundant (5.0 ± 1.4% and 4.9 ± 1.6%, respectively). Cluster 2 consisted of surface samples from July (190709) to November (191115). In this cluster, Favobacteriales (32.7%), Rhodobacterales (11.6%), SAR11_clade (10.8%), Synechococcales (10.6%) and Cellvibrionales (11.7%) were dominant. In cluster 3-1, Synechococcales dominated the communities (32.0%), followed by SAR11_clade (21.3%) and Flabovacteriales (9.4%). In cluster 3-2, SAR11_clade (17.9%), Flavobacteriales (17.3%) and Alteromonadales (10.8%) were dominant. Cluster 3-3 consisted of samples from November (191125) to February (200220) and was characterised by a dominance of SAR11_clade (26.1%) and Flavobacteriales (22.3%). In cluster 4-1, Flavobacteriales (14.6%), SAR11_clade (14.6%) and Parvibaculales (13.9%) were dominant, followed by Alteromonadales (8.2%) and Thiomicrospirales (7.0%). In cluster 4-2, Flavobacteriales (39.6%) and Rhodobacterales (27.4%) were highly dominant. Cluster 4-3 consisted of the sample at 100 m depth in May and Thiomicrospirales (19.3%) was dominant, following Flavobacteriales (25.9%). In cluster 4-4, Flavobacteriales (39.9%) was the most dominant taxa, followed by Rhodobacterales (14.3%), SAR11_clade (10.0%) and Parvibaculales (6.4%).

Figure 7.

Spatiotemporal distributions of bacterial community clusters.

Figure 8.

Taxonomic relative contributions of top 30 abundant Order (on average) to the bacterioplankton faction for each cluster.

SAR11 subclade Ia was dominant in SAR11 clade (53.1 ± 16.1%, Suppl. material 1: Fig. S3) and the contribution of this subclade to the bacterial community accounted for 7.6 ± 5.2% (up to 21.7%, Fig. 9). This subclade was ubiquitously found with relatively high contributions, except for the surface samples from spring to summer, when Flavobacteriales and Rhodobacterales were abundant (Suppl. material 1: Fig. S4). Subclade IV was the second dominant taxa in this clade (23.6 ± 18.6%, Suppl. material 1: Fig. S3) and showed relatively high contributions to the bacterial community in cluster 1 and 3-1, 3-2 and 3-3 (autumn to winter, 8.2 ± 4.6%, up to 18.5%, Fig. 9, Suppl. material 1: Fig. S4), while the relative contribution of subclade IV was small in the other clusters (1.3 ± 1.6%). Subclade II was the third dominant in this clade (15.5 ± 10.0%, Fig. 9). Except for clusters 4-2, 4-3 and 4-4 (0.89 ± 0.26%), subclade II showed relatively stable contributions to the bacterial community in the other clusters (1.8 ± 1.1%). The occurrence of subclade II was positively correlated to that of subclade Ia (n = 74, r = 0.79, p < 0.001, figure not shown).

Figure 9.

Relative abundance of SAR11 subclades in bacterial community in each cluster.

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).

Figure 10.

Spatiotemporal distributions of eukaryotic community clusters.

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.

Figure 11.

Heatmap of taxonomic groups that include phytoplankton community (diatom, dinoflagellates, haptophytes and green algae) at class levels identified by 18S rRNA in Sagami Bay. The colour scale of the heatmap indicates relative abundance in reads between the samples, colour red indicates larger abundance and blue indicates smaller abundance.

Additional cluster analysis was conducted for each taxonomic group that include phytoplankton community (diatom 223 OTUs, dinoflagellates 489 OTUs, green algae 77 OTUs and haptophytes 60 OTUs) and relative abundance top 15 OTUs classified to species level (Suppl. material 1: Figs S5­–S8). Amongst diatoms, Leptocylindrus convexus (16.3%) most dominated in relative reads abundance, followed by Minutocellus polymorphusin (14.2%), Skeletonema menzelii (5.6%), Thalassiosira oceanica (4.2%), Brockmanniella brockmannii (3.0%), Stephanopyxis turris (2.8%), Pseudo-nitzschia delicatissima (2.3%), Thalassiosira tenera (1.4%), Cylindrotheca closterium (1.4%) and Thalassiosira allenii (1.3%). Leptocylindrus convexus was abundant in May from the surface to 100 m depths and June–July in the subsurface layer. Thalassiosira allenii dominated from the surface to 30 m depths in June–August. Cylindrotheca closterium and C. rostratus was abundant in July and especially C. rostratus dominated in the surface layer. Stephanopyxis turris dominated in the subsurface layer during winter (October­­–March) and most dominated in October at 100 m depth. Skeletonema menzelii dominated in winter near the surface. Minutocellus polymorphus, Thalassiosira oceanica and Brockmanniella brockmannii were abundant mainly in winter. Pseudo-nitzschia australis dominated locally on the surface in May and the subsurface layer in June.

Amongst dinoflagellates, Karlodinium veneficum (12.6%) most dominated in relative reads abundance, followed by Heterocapsa rotundata (8.0%), Noctiluca scintillans (5.2%), Pentapharsodinium tyrrhenicum (4.4%), Takayama cf. pulchellum (2.3%), Prorocentrum micans (1.8%) and Lepidodinium viride (1.3%). Karlodinium veneficum was relatively abundant throughout the survey period, especially dominant at 50 m depths in March–June. Heterocapsa rotundata dominated at the surface in May­–September and at the subsurface layer (30–100 m) in November. Noctiluca scintillans was locally abundant at the surface in May and July and near the surface in November. Pentapharsodinium tyrrhenicum was relatively abundant in the subsurface layer throughout the survey period, especially at 100 m depth and Parvodinium inconspicuum at 100 m depths in February.

Sequence reads abundance of haptophytes had an overall peak in April at depths shallower than 30 m. Amongst haptophytes, Chrysochromulina scutellum most dominated in relative reads abundance (13.5%), followed by Prymnesium pigrum (13.3%), Chrysochromulina simplex (10.9%), Phaeocystis globose (10.6%), Prymnesium kappa (3.6%), Phaeocystis cordata (3.0%), Chrysochromulina campanulifera (1.7%), Prymnesium palpebrale (1.5%), Haptolina fragaria (1.3%) and Chrysochromulina strobilus (1.3%). Chrysochromulina scutellum occurred throughout the survey period. Prymnesium pigrum was relatively abundant in winter at depths shallower than 50 m. Chrysochromulina simplex also occurred throughout the survey period, but dominated in May–July and P. globosa at the 100 m depths.

Sequence reads abundance of green algae overall peaked in winter. Amongst green algae, Ostreococcus tauri (%) most dominated in relative reads abundance (27.0%), followed by Micromonas pusilla (20.5%), Bathycoccus prasinos (14.6%), Mantoniella squamata (5.1%), Pycnococcus provasolii (3.3%), Pyramimonas disomata (0.4%) and Pyramimonas australis (0.3%). Ostreococcus tauri dominated in December­–June. Micromonas pusilla and Bathycoccus prasinos dominated in October–July and October–March, respectively.

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 tide-forming 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. 2008, 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 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.

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.

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 community 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 (Kuwahara et al. 2015). The Evenness could be very low in the surface and 100 m layers due to phytoplankton and dinoflagellate blooms caused by short-term 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).