Monitoring took place monthly from March 2018 to June 2019 at five fixed stations (130, 230, 330, 700, and ZG02) (Figure 1) in the BCZ (n = 92), located in the Southern North Sea. BCZ time-series sites were sampled performing several expeditions aboard the RV Simon Stevin (Vlaams Instituut voor de Zee). Both in 2018 and 2019 throughout the spring-summer months one extra monthly cruise was undertaken and samples were collected in order to closely follow the evolution of the phytoplanktonic blooms that occur in this period in Belgian waters (Gypens et al., 2007). Seawater samples were collected at 3 m depth using 4 L Niskin bottles connected to a CTD sensor (Sea-bird SBE25). Physicochemical parameters (temperature, salinity, and major nutrients) were measured and analyzed as part of the national water quality monitoring programs LifeWatch (Mortelmans et al., 2019a) years with the methodology found in Mortelmans et al. (2019b). The daily global solar radiation data was collected at the Oostende station of the Royal Meteorological Institute of Belgium, transformed to Photosynthetic Active Radiation (PAR, μmol m^–2 s^–1) making use of the empirical relation found in Rousseau et al. (2000) and biweekly averaged.

For all above-mentioned expeditions the identical sampling methodology was followed. 100 mL of water from Niskin bottles was stored and frozen in vials for the analysis of each of the nutrients such as: NO₃^–, NO₂^–, NH₄⁺ [Dissolved Inorganic Nitrogen (DIN) = NO₃^– + NO₂^– + NH₄⁺], dissolved inorganic phosphate (DIP = PO₄^3–) measured as soluble reactive phosphorus, and dissolved inorganic silica (DISi). For chlorophyll-a (Chl-a) analysis 100–500 mL of water were filtered using glass-fiber filters (Whatman GF/F) and immediately frozen at −20°C. Chl-a was quantified fluorometrically in the laboratory using acetone 90% for extraction (Strickland and Parsons, 1972).

Additional 50 mL of water from St. 330 were taken and preserved with 40% formaldehyde (final concentration 2%) for heterotrophic bacteria counts. Bacteria were DAPI stained following Porter and Feig (1980) method and enumerated under epifluorescence microscopy. Biovolumes were calculated by treating cocci and rods as cylinders and spheres, and finally converted to C-biomass using the biovolume-dependent conversion factor established by Simon and Azam (1989).

Considering that mixoplankton is largely widespread in both nano and microorganisms, the field sampling technique used in this study targeted the marine protist plankton sizing >0.22 μm. The DNA samples for the study of the protistan community were collected vacuum filtering 500–800 mL of water (from Niskin) through 0.22 μm polycarbonate filters (47 mm) and storing the samples immediately at −20°C. Total DNA was extracted from filters using NucleoSpin Soil extraction Kit (Macherey-Nagel, Düren, Germany) following manufacturer’s protocol. For a maximum efficiency of the extraction from the filters a sample lysis step was added using 10 mL cryotubes (using a high velocity bead beater for 10 min). Up to three filters were pooled and used for DNA extraction when there was not found sufficient biomass on a single filter. Standard polymerase chain reactions (PCR) were performed to amplify the universal eukaryote SSU 18S rRNA gene. Primers TAReuk454FWD1 (5′-CCAGCASCYGCGGTAATTCC-3′), TAReukREV3 (5′-ACTTTCGTTCTTGATYRA-3′) were used to target the V4 region of the 18S rRNA gene (Stoeck et al., 2010). PCR reactions performed had a total volume of 25 μL, containing 2.5 μL of microbial DNA (5 ng/μL), 5 μL of both amplicon forward and reverse primers (1 μM) and 12.5 μL of high-fidelity polymerase (Kapa Biosystems). Plates were sealed and the following PCR-program was run in a thermal cycler: initial denaturation at 95°C for 3 min, followed by 25 cycles of 95°C for 30 s, annealing at 55°C for 30 s; extension at 72°C for 30 s final extension at 72°C for 5 min. All PCR products (480 bp, ∼383 bp + 97 bases of primers) were verified on a 1.5% agarose gel. The following library preparation of 18S ribosomal RNA gene amplicons was performed: PCR clean-up 1, index PCR, PCR clean-up 2, library quantification, normalization and pooling following the 16S Metagenomic Sequencing Library Preparation guide (Illumina, 2013). Library denaturing and sample loading to the Illumina MiSeq system was performed to perform a 300 bp paired-end sequencing using V2 chemistry.

The reads were denoised and merged with DADA2, v1.16. (Callahan et al., 2016) and annotation reads were subsequently classified with assignTaxonomy, the DADA2 implementation of the naive Bayesian classifier method. We used the Protist Ribosomal Reference Database PR2 (Guillou et al., 2012) for the taxonomic assignment of sequences. Preparation of the amplicon sequence variant (ASV) tables was done in R v4.0 (R Core Team, 2018). The number of final sequenced reads per sample ranged between 23,299 and 52,256.

As this study focuses on protists, taxa not belonging to eukaryotic unicellular (protists) plankton were removed (including fungi, cnidarians, and metazoans), as well as singletons and all low abundance ASVs (the ones not seen more than three times in at least 10% of the samples). Final curated results produced ∼5 million reads that were grouped into 686 ASVs representing the protistan communities in the BCZ.

Trophic strategies were manually annotated and curated based on the current accepted forms of protistan plankton nourishment (Flynn et al., 2019). Species-specific knowledge obtained from published literature were used (Armeli Minicante et al., 2019; Faure et al., 2019; Leles et al., 2019; Schneider et al., 2020a) to classify the species into the next trophic groups: proto- phytoplankton or strict autotrophs (AU), proto- zooplankton or strict heterotrophs (HET), constitutive-mixoplankton (CM), non-constitutive mixoplankton (NCM), and Unknown (NA). Similarly, ASVs were grouped in functional taxonomic groups based on the assigned taxonomical information: dinoflagellates (Dinophyta), ciliates (Ciliophora), haptophytes (Haptophyta except Phaeocystis genus), diatoms (Bacillariophyta), cryptophytes (Cryptophyta), heterotrophic nanoflagellates (comprised by mainly Cercozoa, Protalveolata, Choanoflagellida, and Picozoa), green-brown microalgae (mainly Mamiellales, Prasinophyceae, Chrysophyceae), Phaeocystis (referring to Phaeocystis genus), rhizarians, and unclassified.

We used the compositional approach described in Gloor et al. (2017) to analyze our metabarcoding dataset. First, a Bayesian-multiplicative treatments of zeros was performed on the ASV table using cmultRepl() function of the zCompositions package (v1.3.4): this uses sample-wise totals to convert zero counts (which will lead to errors in log-ratios) into near-zero estimates, assuming undersampling rather than absence. Next, a centered log-ratio transformation (clr) was applied to the dataset (Aitchison, 1986) using “microbiome” package (Lathi and Shetty, 2017). All statistical analyses were conducted in R, version 4.0.2 (R Core Team, 2018) and the following packages were used for data visualization: “ComplexHeatmap” (Gu et al., 2016) “ggplot2” (Wickham, 2011).

We examined the distribution of samples along the time-series using a principal component analysis (PCA) on our microbial diversity metadata. Microbial community dissimilarities (β-diversity) between seasons, months, sampling dates (cruises), and stations were tested with PERMANOVA (Anderson, 2001) on Aitchison distance of the clr-transformed ASV tables. Ultimately, we used the adonis() function along with a beta dispersion test. The latter was used to evaluate the homogeneity of dispersion by the application of the betadisper() function in “vegan” package (Oksanen et al., 2019).

Read count proportions of the most abundant species occurring in our dataset (auto-, hetero-, and mixotrophs) were displayed using a heatmap visualization. Proportions were relative to the sum of all observed ASVs corresponding to the each sampling date, which corresponded to five station averaged data. In order to enhance the understanding of mixoplankton temporal diversity within a greater context (i.e., seasonality and temperate seas) samples were classified by seasons respecting the season start days (which showed statistical significant differences, see section “Results”). Species with sample reads proportions reaching more than 5% were considered as significant contributors to the community.

Alpha diversity (α-diversity) was determined by calculating Shannon entropy and species richness for all the samples. Shannon diversity measures reflects true diversity (richness + evenness) and is less susceptible to fluctuations in rarer phylotypes. Trophic group’s specific Shannon and species richness were as well calculated for autotrophs, heterotrophs and both mixoplankton functional types (CM and NCM). The correlations of CM and NCM richness and α-diversity with environmental contextual data was performed constructing Pearson correlation matrices that were calculated and plotted using the corrplot() function of the “corrplot” package (v0.84) and adjusted for multiple testing. We calculated Variance Influence Factor (VIF) to check for multi-collinearity of each of the predictor variables. In order to avoid multi-collinearity between highly correlated variables (i.e., nutrients), we removed the variables that had a high VIF (>5).

The spatial and temporal patterns of environmental data were checked for normal distribution and spatially (by stations) and temporally analyzed by multiple statistical tests (Tukey, Kruskal–Wallis, Pairwise Wilcoxon). To explore trophic groups’ beta diversity (β-diversity) and its relation to environmental and contextual variables (SST, salinity, PAR, DISi, DIP, DIN, Bacteria biomass, and Chl-a) we performed a set of redundancy analyses (RDA) using the clr-transformed ASV tables as response matrices and tables of environmental variables as explanatory matrices. We performed stepwise model to identify most significant environmental variables using ordiR2step() function in “vegan.”

The environmental conditions widely varied across the temporal series in the BCZ (Figures 2A,B), and generally the values recorded in the study period follow the seasonal dynamics already described for the area. PAR values significantly increased in spring and summer (up to 903.0 ± 222.7 μmol m^–2 s^–1) and were lowest in winter months (185.8 ± 100.3 μmol m^–2 s^–1) (Figure 2A). The yearly mean incident PAR values for the years 2018 and 2019 were the two highest ones of the last 35 years with an increase of 13.6% (2018) and 11% (2019) to the calculated median for this period (Oostende station of the Royal Institute of Meteorology, Belgium, not shown). SST dynamics closely followed the light increase effects, reaching the peaks in July 2018 and June 2019 (21.3 ± 0.7°C and 16.3 ± 0.9°C) and the minimum temperature values were found in March 2018 and January 2019 (2.7 ± 0.9°C in 2018 and 5.4 ± 0.6°C in 2019) (Figure 2A). Salinity values in the BCZ differed most between summer and autumn 2018, 33.2 ± 1.5 and 32.2 ± 1.6, respectively (data not shown), but were not significantly different (Kruskal–Wallis, p > 0.05). The only environmental contextual data that showed significant differences among stations was salinity (Kruskal–Wallis, p < 0.05).

Inorganic nutrient concentrations followed the typical temperate sea seasonal dynamics (Figure 2B), with higher concentrations in winter (DIN: 27.2 ± 16.3 μM; DIP: 0.57 ± 0.3 μM; DISi: 14.0 ± 9.5 μM) than in summer (DIN: 2.2 ± 1.7 μM; DIP: 0.28 ± 0.2 μM; DISi: 3.1 ± 1.3 μM), and with the maximum peaks just before the early spring phytoplanktonic bloom around February–March in both years 2018–2019. The highest bacterial biomass peak measured throughout the time-series attained 137 mg C m^–3 in March 2018 (Figure 2C). Chl-a concentrations fluctuated between 0.5 and 32.1 μg L^–1 in the BCZ (Figure 2C) and the station averaged values showed three peak periods with maximum values of 27.3 μg L^–1 in March 2018, 30.6 μg L^–1 in August 2019 and 32.1 μg L^–1 in April 2019. The lowest Chl-a concentrations were found in August and December 2018.

Seasonality was inferred when all the sampled (n = 92) were mapped using a PCA on the Euclidean distance of clr transformed data (Figure 3). Four distinct seasonal clusters were observed, placed in a circular pattern (temporal). Statistical analysis showed significant differences among seasons (PERMANOVA, p < 0.001). The similarities between summer and autumn samples were noted, and community variations within a season between different months. We performed PCAs for all samples, clustered by sampling dates, stations and seasons independently (Supplementary Figure 1). Whereas statistically significant differences among sampling dates and months were found (PERMANOVA, p < 0.01), there were no significant differences among stations (PERMANOVA, p = 0.97). Similarly, there were no significant differences found between interannual variations among samples. Since there were no spatial differences found among communities (between stations), time-series visualization of the proportional number of reads was done on the five stations averaged values (Figures 4, 5).

The temporal succession analysis of plankton groups’ (in relative abundances) (Figure 4A), showed that dinoflagellates were major contributors to the dataset with 22.8–82.2% of the assigned ASV reads. The maximum peak of autotrophic sequence reads coincided exactly with the Phaeocystis (group) proportion increase within the community in April and May 2019 samples, that reached up to 18.6–20.9%, respectively (Figure 4A). In addition, Chl-a concentration peaks during the sampling period coincided in time with the Phaeocystis increments observed during spring (greater in 2019 than in 2018) (Figures 2C, 4A). Besides these two peak periods, the contributions of Phaeocystis to the samples were less than 4% of the reads.

Diatoms ranged from 4.9% (April 2018) up to 26.6% (June 2019). Heterotrophic nanoflagellates and ciliates contributed significantly to the sequence reads, attaining maximum values up to 26.3% (November 2018) and 19.1% (January 2019), respectively. Green-brown microalgae, rhizarians and the ensemble of haptophytes-cryptophytes associated reads accounted less than 7% to the total counts throughout all seasons (Figure 4A). Similarly, unclassified taxa were <10% all along the sampling period, however, they had a significant contribution in May 2018 and 2019 (Figure 4A).

Metabarcoding results revealed that dinoflagellates (Dinophyta) was the taxa comprising most mixotrophic species, accounting for 70.9% of the total mixoplankton reads and a total of 70 identified ASVs out of 169 ASVs (Supplementary Data 1). The second major mixoplankton contributors were ciliates with 13.5% of the mixoplankton reads and 36 assigned ASVs. The remaining taxonomic groups containing mixoplankton contributed all together <20% to the total mixoplankton associated read counts and corresponded to a total of 36 ASVs.

We reported the temporal proportions of the trophic strategies adopted by two major taxa comprising mixoplankton species, dinoflagellates, and ciliates (Figures 4B,C). The average proportion of mixotrophic dinoflagellates and ciliates over the studied period was about 49.5 and 41.0%, respectively. Nevertheless, differences lie in the MFTs associated to each taxonomic group: primarily CM for dinoflagellates and NCM for ciliates. Indeed, dinoflagellates classed as strict autotrophic (Supplementary Data 1) or as NCM did not account for more than 5% of the reads throughout the temporal samples and no ciliate species were classed as CM or strict autotroph (Figures 4B,C). Trophic strategies in both in dinoflagellates and ciliates did not show a clear seasonal pattern, however, heterotrophic dinoflagellates attained >50% of the reads during spring 2018 temporal samples (Figure 4B), with a peak of 81% in April 2018. In the case of ciliates, April 2018 temporal sample showed particularly a high proportion of counts belonging to NCM (90.5%).

The protistan plankton taxonomical diversity and seasonal succession in the BCZ were illustrated and expressed in proportional number of reads (five stations averaged) emphasizing two major MFTs (Figure 5). The Shannon α-diversity indexes calculated for the entire community differed seasonally (Figure 5A), reaching the highest values during summer-winter 2018, with a maximum value of H’ = 4.73 in November 2018 and lowest values in April 2019 H’ = ∼ 0.9. Proportions of the community trophic strategies across the time-series samples (Figure 5B) revealed that highest CM proportions occurred between May–August 2018 (32–34% of the reads), and lowest in April 2018 and May–June 2019 (15–19%). CM proportion average throughout the sampling period attained up to 24.3% amongst the four assigned trophic modes. NCM was the least frequently occurring trophic strategy in our dataset, attaining maximum values of 15.1% in September 2018 and averaging 4.7% throughout the whole sampling period. In contrast, proportions of reads assigned to heterotrophic species were often dominating the community (>40% in 10 temporal samples), averaging 40.8% of the reads and reaching up values up to 72.2% in April 2018. Autotrophic read percentages attained a maximum of 53.4% in May 2019 and averaged 28%.

Out of the 10 most abundant non-mixoplankton ASVs 4 were assigned as autotrophs and 6 as heterotrophs (Figure 5C; blue color scale). The top 15 CM and NCM identified by our molecular dataset are shown below (green and orange color scales, respectively). The low diversity index measurements during spring 2018 matched with the important single-species contributions to the community. The heterotrophic dinoflagellate Gyrodinium spirale (Figure 5C) ranged from 29.3 to 62.7% in proportional number of ASV reads from April to May 2018. A slight presence of the characteristic bloom-forming species of the BCZ -Phaeocystis globosa- was noted (up to 5.8%). Heterocapsa pygmaea was the most abundant CM species (15–16%) together with Alexandrium minutum (7.9–12.7%). CM dinoflagellates Tripos concilians and Tripos fusus proportions reached non-negligible values (∼8–9%) in two temporal samples during spring 2018 (March and May). The only NCM species that had a significant (>5%) contribution to the community during this season was a ciliate belonging to Strombidiidae family (ASV 33: Strombidiida_B_XX_sp.).

In spring 2019 Phaeocystis globosa exerted the dominant role accounting up to 18–20% from April to May, together with high proportional abundance in May 2019 of autotrophic Chloroparvula pacifica (20.9%) and a slight presence of the diatom Thalassiosira pseudonana (∼7–9%). Summer and autumn in 2018 were characterized by a strong presence (avg. = 11.9%) of CM Alexandrium minutum, as well other “Non-abundant CM,” that contributed substantially to the community (∼7–12%). The molecular dataset revealed the highest NCM presence of the time-series in September 2018, led by Acanthometron sp. and the dinoflagellate Dinophysis acuminata. Several ciliates belonging to the Strombidium genus (i.e., Strombidium caudatum 1.9%) contributed modestly to the observed NCM peak.

Heterocapsa pygmaea presented the highest CM proportions during winter 2018–2019, sharing similar proportional read abundances with Alexandrium minutum (max. = 8%). However, most abundant non-mixoplanktonic organisms presented higher proportional number of reads such as the heterotrophic dinoflagellates belonging to Gyrodinium genus (max. = 13.9%).

Non-constitutive mixoplankton, overall, had a relatively low contribution to the total number of reads throughout all the studied period compared to CM (Figure 5C), with a total maximum contribution to the entire protistan community of 7% in September 2018. The most important constituents were marine oligotrich ciliates belonging to the genus Strombidium, ranging from 2 to 7% in summer 2018 and 2.6% and attaining up to 5.3% in April 2018.

Mixoplankton observed species richness (Obs.) and Shannon α-diversity (H’) throughout the time-series was significantly lower (Kruskal–Wallis, p < 0.01) than those observed for autotrophs and heterotrophs (see Supplementary Table 1 and Supplementary Figure 2). CM richness and H’ attained maximum values during winter (Obs. = 48.4 ± 4.7; H’ = 3.0 ± 0.2), and the lowest values were found in spring 2018 (Obs. = 2.1 ± 1.8; H’ = 1.2 ± 0.5). NCM observed diversity was significantly lower (Pairwise Wilcoxon, p < 0.01) than for CM (maximum Obs. = 15.4 ± 1.5; H’ = 2.2 ± 0.3) but Shannon diversity values did not differ. The temporal dynamics of the diversity values (Obs. and H’) followed the same seasonal pattern for both MFTs. Relative to other functional groups, heterotrophs presented the highest α-diversity and species richness, followed by autotrophs.

Mixoplankton in functional types specific Shannon α-diversity and ASV richness Pearson correlation matrices (n = 92) with environmental contextual data for the BCZ time-series are shown in Supplementary Figures 2, 3, respectively (autotrophs and heterotrophs were analyzed as well). DISi and DIN presented high VIF (>5, results not shown) and were removed from the analysis to avoid multi-collinearity. CM diversity and richness showed stronger positive correlation (r = 0.33–0.55, p < 0.01) with major inorganic nutrients (represented by DIP) than NCM (r = 0.33–0.45, p < 0.01). In contrast, PAR showed stronger negative correlation with CM Shannon diversity and observed richness (r = −0.70 to 0.55, p < 0.01) than NCM (r = −0.44 to 0.45, p < 0.01). Whereas both CM richness and H’ showed a significant negative correlation with bacterial biomass, NCM Shannon diversity showed a weak correlation. Salinity and Chl-a values were not strongly correlated with mixoplankton diversity or species richness (both CM and NCM), however, we found a weak negative correlation (r = −0.27, p < 0.01) with observed heterotroph species.

In the beta diversity (β-diversity) patterns performed for CM, NCM autotrophs and heterotrophs across the time-series (Figure 6), seasons (spring-winter and autumn-summer) were well separated along the first axis of our redundancy analysis (RDA) for each trophic group (between 13.2 and 14.8% of variance constrained). There was a strong orientation of temperature toward all summer samples in all the four trophic strategies displayed, and dissolved inorganic nutrient concentrations (DIP, DISi, DIN) were found to be associated with autumn and winter samples. High Chl-a concentrations were correlated with spring samples. Autotrophs β-diversity was found to be more impacted by Chl-a and inorganic nutrients than the remaining trophic groups (Figure 6C). Salinity resulted to be the environmental factor with low impact at the community shifts. It was observed that the samples for each of the trophic strategies were differently ordinated in our RDA plots. Shortest distances were observed between spring and summer sites for NCM (Figure 6B). Across all trophic groups, the autumn samples showed the largest spread among the first RDA axis and partially overlapped with summer and winter samples. Secondary axis inferred the separation of winter-autumn against spring-summer seasons. In our multivariate analysis, β-diversity for the different trophic groups was controlled by multiple environmental variables that shaped microbial communities over the time-series.

For the presented work, we gathered new information regarding the spatial and temporal variability of the most relevant MFTs and marine protist plankton (>0.22 μm) in the North Sea based on the analysis of the V4 region of the 18S rRNA gene. Species trophic modes were assigned to depict the taxonomic composition and diversity of mixoplankton with high resolution, focusing on two MFTs.

While microscopy on mixoplankton research is still often deployed (Flynn et al., 2019), here, we used a Next Generation Sequencing approach by means of MiSeq Illumina sequencing, that provides the potential for rapid and more comprehensive automated examination of environmental samples (e.g., Mäki et al., 2017) to investigate the MFTs dynamics in the marine environment. Mostly dependent on taxonomic expertise for a species level resolution (Kase et al., 2020), optical microscopy methods are time-consuming (Stern et al., 2018), and certain small-sized taxa occurrence like dinoflagellates, which are generally mixotrophs, are underestimated (Stern et al., 2018; Kase et al., 2020). Furthermore, due to mixoplankton resemblances with several strict auto- or heterotrophic protistan plankton, misclassification has been relatively common due difficulties to assess phagotrophic and/or autotrophic capabilities (Anderson et al., 2017; Beisner et al., 2019). These reasons have led to a disrupted understanding of the nature of most marine planktonic primary producers, that in fact, have the ability to potentially express mixotrophy (Leles et al., 2019).

DNA analysis opens the possibility to partially replace or rather complement microscopy but just like optical methods, it has its limitations. In this work, the quantification of planktonic protists was visualized based on sequence amplicon data relative abundances, being aware of the difficulties of data interpretation of approaches based only on proportions (Kase et al., 2020), which has largely been discussed in recent literature (Weisse et al., 2016; Vasselon et al., 2017; Santoferrara, 2019). Besides the inherent technical issues occurring due to the method: sample preservation (Mäki et al., 2017), DNA extraction (van der Loos and Nijland, 2021), primer choice and specificity (Elbrecht and Leese, 2015) and polymerase chain reaction (PCR), protists have a strong variation in cell density, cell biomass and 18S rRNA gene copy number. This natural variability strongly varies between different taxonomic groups and species (Zhu et al., 2005) and has been already shown for dinoflagellates (Galluzzi et al., 2004), diatoms (Connolly et al., 2008; Godhe et al., 2008) or ciliates (Gong et al., 2013). For these reasons, sequencing results based only on the 18S rRNA gene proportions might often be misinterpreted since they differ with relative abundances revealed by microscopic counts (Santi et al., 2021).

In fact, microbiome datasets are compositional (Gloor et al., 2017) and therefore, standard methods of analysis are not applicable and validity of all downstream analyses is critical (Weiss et al., 2017). The total number of reads obtained per sample does not reflect the absolute number of microorganisms present, however, visualization of relative abundances of ASVs throughout a time-series is highly intuitive for the most abundant taxa present in metabarcoding datasets. Moreover, this allows a comprehensive understanding of the taxonomic profile (i.e., occurring species).

Although we cannot fully rely on a quantitative picture of microbial communities based on amplicon data (Lin et al., 2018), the DNA-barcoding approach was able to provide advantages over traditional methods to study MFTs, as, event if not quantitative, it gives a taxonomically enhanced picture of the whole protistan community, overrepresenting taxa on which mixoplankton mostly occurs such as dinoflagellates (Gran-Stadniczeñko et al., 2019). Besides, metabarcoding only analyzed molecules and not per se the organism itself, this means that the taxonomic assignment of a given ASV depends strongly on the reference database used (i.e., PR2). Nevertheless, it has recently been shown that the environmental drivers shaping protist communities are the same for molecular (metabarcoding) and quantitative microscopy methods (Piwosz et al., 2020), and therefore are reliable in the context of ecological interpretations.

In the BCZ, the previous knowledge focusing on protistan plankton relies on microscopy (Rousseau et al., 2000, 2006; Muylaert et al., 2006) or flow cytometry (Martínez et al., 2020), some of which include quantitative phytoplankton counts since 1968 to 2010 (Nohe et al., 2018). Traditionally, trend analyses and dynamics of autotrophic plankton biomass have been often based on pigment (Chl-a) concentration. Nevertheless, an important fraction of Chl-a containing protists are well known to have the ability to perform phagoheterotrophy (Flynn et al., 2019). Long-term marine protist plankton surveys have permitted to characterize in detail the recurrent annual phytoplankton events in BCZ: (1) late February–March diatom spring bloom, (2) directly followed by a huge bloom of the HAB Phaeocystis globosa in April–May, and (3) a late summer diatom bloom (Rousseau et al., 2002; Nohe et al., 2020). Regarding the mixoplankton diversity dynamics, we found that CM (mainly dinoflagellates) were highly present throughout the year, and often dominating the amplicon counts of our dataset. Similar results have already been revealed (for neighboring regions to the BCZ) in terms of absolute cell counts under light microscopy (Leles et al., 2019; Schneider et al., 2020b).

The metabarcoding approach did not reflect the above-mentioned diatom blooms in proportional number of reads (they were identified through microscopy survey, data not shown). Phaeocystis globosa blooms, however, were much more accurately represented (mainly in 2019). In contrast, our approach was able to emphasize mixoplankton, that were present and moderately abundant in all samples, suggesting that they coexist (Edwards, 2019) with strict autotrophic and heterotrophic protists under various environmental conditions in the BCZ. Similarly to previous studies (Taylor and Cunliffe, 2014; Massana et al., 2015; Kase et al., 2020), almost half of our sequence assemblages belonged to dinoflagellates, of which most were classified as CM and heterotrophs, and a smaller number as strict autotrophs and NCM. Compared to the main traditional long-term observations in BCZ based on microscopy (Nohe et al., 2018; Martínez et al., 2020) our results showed a higher taxonomic resolution resulting mainly in an increased identification small sized dinoflagellates (<20 μm).

Throughout the studied period CM-dinoflagellates genera Heterocapsa, Alexandrium, Tripos, and Karlodinium stood out within the mixoplankton read counts. Nohe et al. (2020) noted an increase of various dinoflagellate genera such as Alexandrium or Gyrodinium (heterotrophic) since the 1970s. However, it is possible that historical surveys that describe the area were underestimating dinoflagellate occurrence in the environment, due to resolution limitations of the identification methods used (Kase et al., 2020).

Non-constitutive mixoplankton were very scarce compared to other trophic groups and mostly emerged in summer-autumn under low-nutrients high-light irradiance conditions. Two hypothesis may explain this observation. First, there is lower predation by higher trophic levels during this period of the year and second, specific prey types for NCM (which are often restricted to a particular prey type) might be thriving during the period prior to the NCM abundance increase. NCM were almost exclusively represented by ciliates, and within ciliates, it was found that ∼ 40% on average were classed as mixotrophs. Our results are in accordance with former observations of Dolan (1992), that showed that chloroplast-retaining oligotrichs rarely dominate the ciliate community. In fact, the low biomass of NCM throughout the year in temperate seas was demonstrated by Leles et al. (2020), and previous studies on ciliates in oligotrophic waters showed similar proportions of mixotrophy (Dolan and Pérez, 2000; Stoecker et al., 2009; Romano et al., 2021).

Even though we were limited by the sampling strategy performed for heterotrophic bacterial biomass (data available only for one station of the BCZ), we explored its dynamics along with mixoplankton diversity. In fact, our dataset revealed a high bacterivorous mixoplankton genera (i.e., Alexandrium, Heterocapsa, Karlodinium). Heterocapsa, for instance, (i.e., Heterocapsa triquetra) is known to perform strong bacterivory (up to 12.3 ± 2.9 bacteria Heterocapsa^–1 h^–1) under light or nutrients limitation (Millette et al., 2017) as well as Alexandrium (Lim et al., 2019) and both can form large blooms (Galluzzi et al., 2004; Seong et al., 2006; Millette et al., 2015). Indeed, the importance of mixoplankton bacterivory (Zubkov and Tarran, 2008; Unrein et al., 2014) has been discussed to alleviate nutrient (Stoecker et al., 2017) or light stress (Fischer et al., 2017). Protists can strongly impact bacterial diversity (Suzuki, 1999) but the understanding of how mixoplankton impacts the bacterial biomass and nutrient cycling in the microbial loop remains unclear.

Like auto- and heterotrophs, mixoplanktonic diversity increased when resources were available (Edwards, 2019); yet, abundant light, nutrients or prey seemed to favor more strict autotrophic and heterotroph diversity (Mitra and Flynn, 2010). Our study revealed, in general, slight differences in α- and β-diversity temporal patterns for the studied microbial trophic groups. Schneider et al. (2020b) showed that mixoplankton had an advantage over strict autotrophs and heterotrophs in low biomass environments, nutrient-depleted areas and/or seasonally stratified water columns. The BCZ, however, is considered to have a well-mixed water column throughout the year (Masquelier et al., 2011). Therefore, when correlations of mixoplankton α-diversity with environmental data were investigated, species richness and Shannon diversity showed a stronger positive correlation with nutrient rich waters for CM than NCM, outlining the differences in nutrient acquisition strategies by each of the MFTs (Flynn et al., 2019).

More CM species prospered in nutrient rich and low Chl-a (proxy for biomass) environments and they did co-occur in time with strict autotrophs (i.e., Phaeocystis globosa). In contrast, strong light incidence and high temperature environment (low biomass and nutrient poor waters) favored NCM species diversity over CM. Overall, nutrient increase drove diversity and species richness within all trophic groups, particularly in autumn and winter, with DIP concentration changes being the most important factor driving this variability. The decrease in bacterial biomass was linked both to MFTs and heterotrophic richness α-diversity decrease, however, autotrophs were not correlated to bacterial biomass. Some of the above-mentioned correlates might be drivers of MFTs diversity and richness. In agreement to recent published modeling studies (Anschütz and Flynn, 2020; Leles et al., 2020) our molecular dataset demonstrated a niche separation between CM and NCM according to nutrient, prey (bacterial biomass) and light availabilities.

The multiple environmental contextual variables shaping microbial community shifts were located in similar places for the four groups, indicating that such seasonal changes in β-diversity for autotrophs, heterotrophs and mixoplankton, were controlled by the same physicochemical processes (temperature, light, and nutrient availability). These processes are directly affected by nutrient influx (from rivers, currents), remineralization or interspecific competition for resources (Da̧browska et al., 2020). Similar ecological analysis across different oceanic regions (Hörstmann et al., 2021) has already proved useful for the delineation of oceanographic provinces (ecosystem boundaries) and determination of important features generating microbial diversity changes. In our study, the time-series framework allowed to stress the importance of seasonality in BCZ for these community changes.

Ocean currents can drive overlaps in β-diversity signals: In an effort to understand the factors driving the differences between both MFTs (CM and NCM), auto- and heterotrophic plankton contrasting β-diversity patterns on our seasonally arranged RDA analysis for the microbial trophic groups associated under the same environmental conditions, results displayed for trophic groups.

We noted overlapping of in β-diversity signals: the light differences among the different trophic groups shows the value of microbial diversity data in holistically describing, especially temporal series seasonally arranged.