Seasonal temporal dynamics of marine protists communities in tidally mixed coastal waters

Major seasonal community reorganizations and associated biomass variations are landmarks of plankton ecology. However, the processes determining marine species and community turnover rates have not been fully elucidated so far. Here, we analyse patterns of planktonic protist community succession in temperate latitudes, based on quantitative taxonomic data from both microscopy counts and ribosomal DNA metabarcoding from plankton samples collected biweekly over 8 years (2009-2016) at the SOMLIT-Astan station (Roscoff, Western English Channel). Considering the temporal structure of community dynamics (creating temporal correlation), we elucidated the recurrent seasonal pattern of the dominant species and OTUs (rDNA-derived taxa) that drive annual plankton successions. The use of morphological and molecular analyses in combination allowed us to assess absolute species abundance while improving taxonomic resolution, and revealed a greater diversity. Overall, our results underpinned a protist community characterised by a seasonal structure, which is supported by the dominant OTUs. We detected that some were partly benthic as a result of the intense tidal mixing typical of the French coasts in the English Channel. While the occurrence of these microorganisms is driven by the physical and biogeochemical conditions of the environment, internal community processes, such as the complex network of biotic interactions, also play a key role in shaping protist communities.

The DNA extraction and generation of metabarcodes were performed using the exact 242 same procedure for all samples. Samples were first incubated 45min at 37°C with 100µL 243 lysozyme (20mg/mL), and 1h at 56°C with 20µL proteinase K (20mg/mL) and 100 µL SDS 244 20%. Nucleic acids were then extracted using a phenol-chloroform method (Sambrook et al. 245 1989), and purified using silica membranes from the NucleoSpin® PlantII kit (Macherey-246 Nagel, Hoerdt, France). DNA was eluted with 100µL Tris-EDTA 1x pH8 buffer and  based on abundance and occurrence) was checked and refined manually by BLASTing them 297 (https://blast.ncbi.nlm.nih.gov/Blast.cgi) against the SILVA reference database 298 (https://www.arb-silva.de/). The origin and assignations of the best blast sequences (most of 299 which were 100% similar to our sequences) and of the corresponding strains or isolates were 300 carefully examined before taking the final taxonomic assignation decision (Table S1).    Variance partitioning analyses allowed to filter out the variations due to temporal structures, 344 or autocorrelation, which accommodate the use of statistical tests to further assess which 345 environmental variables can influence community dynamics and species composition. All 346 parameters were first tested for collinearity, then successively used in a forward selection to 347 identify those significant to be tested for the study. To interpret temporal variations, we 348 calculated Spearman's rank correlation coefficients between the environmental parameters 349 and the eigenvalues of the first three axes of the RDA.  360 At the SOMLIT-Astan time-series station (Fig. 1), as expected in temperate marine

465
A closer examination of the diatoms species dynamics was achieved by calculating 466 mean monthly relative abundances of diatoms OTUs or cells and selecting the 5 most abundant taxa for each month in the microscopic and metabarcoding datasets (Fig. 5B, C). 468 The resulting pools of cells or reads selected accounted for >75% and 70% of all counts/reads 469 in these two datasets, respectively. In microscopic counts, autumn and winter assemblages  The use of a dbMEM analysis to decompose the temporal patterns of the community 500 allowed us to detect and investigate the environmental and biological processes involved in 501 the control of protist assemblages' dynamics at different timescales (Fig. 6). The dbMEM 52.2% of the OTUs variability in community composition, respectively. As expected, 505 seasonality -expressed in the first 2 constrained axes of the RDA -explained most of the 506 observed temporal variability (RDA1: 19.8-17.8% and RDA2: 11.5% and 9.3%, for 507 morphological and metabarcoding datasets respectively; Fig. 6B, D). For both datasets, the 508 winter and summer assemblages on the one hand, and the autumn and spring assemblages on 509 the other, were clearly distinguished on axes 1 and 2. Spring assemblages showed more 510 interannual variability, especially when the morphological dataset was considered (Fig. 6A).

511
The annual cycle was better delineated when the metabarcoding dataset was considered (Fig.  6C). For both datasets, the taxa/OTUs with the highest RDA1 and RDA2 scores 513 corresponded to dominating species (section 3.3 and Fig. 5) and displayed clear seasonal 514 variations in terms of cells or reads abundances (See Table 1 and Fig. 7A, B). For the 515 morphological datasets, the pelagic chain forming Guinardia delicatula and Thalassiosira 516 levenderi/minima and the benthic or tychopelagic taxa Fragilaria/Brockmanniella, Paralia 517 sulcata and Psammodictyon pandutiforme had the highest scores for RDA1 and/or RDA1.

518
For the metabarcoding dataset, the OTUs with the highest RDA1 and 2 scores also included 519 G. delicatula, but pointed as well to nanoplanktonic diatoms (such as Minidiscus comicus), 520 and to species belonging to other phyla or classes such as the Dinophyceae and Cercozoa, all 521 displaying strong seasonality ( Fig. 7B and Table 1).

522
The axis 3 of the RDA (4.8 and 3.9 % of the variance explained for the morphological  To investigate the environmental factors that primarily drive seasonal protist 532 assemblages, we calculated Spearman rank correlation coefficients between the potential 533 explanatory variables and the first 3 axes of the RDA (Fig. 8A, B). Here, we considered the environmental variables selected by forward selection for both datasets, namely: temperature,   Our dbMEM analyses allowed us to uncover a clear and persistent biannual rhythm (axis 3 of 600 our RDA) in the protist community dynamics. This biannual pattern was also observed in the 601 monthly turnover rates of species over the 8 years (Fig. S4). However, no clear relationship 602 was found with annual variation in tidal amplitude which was not among the environmental 603 factors selected for correlation analyses (Fig.8). 604 We found a strong correlation between the first two axes of the RDA and temperature, 605 macronutrients and light (Fig. 8)    Bi-annual variations in the protist community dynamics were also identified from 727 analyses of both dataset (Fig. S4 and Fig. 6). In some ecosystems, rhythmic depletions of  Intrinsic, plankton self-organization processes are also involved in these annual oscillations.
In environments such as the coastal waters of the EC that support one of the busiest shipping     Location of the study area. The SOMLIT-Astan sampling station (48:46'49'' N; 3:58'14'' W) is located in the Western English Channel, 3.5 km from the coast. The water column at this site is 60 m deep and is never stratified due to intense tidal mixing. The site is strongly impacted by storms in Winter. .

Figure 2
Monthly variations of the hydrological and meteorological parameters at the SOMLIT-Astan time-series station in the period 2009-2016. Sampling was carried out at high neap tides. PAR is the photosynthetically available radiation calculated as the average light received during the 8 days before sampling. Kd490 is intended as the diffuse attenuation coefficient for downwelling irradiance at 490 nm. Interannual variations of all parameters presented in this figure can be found in figure S2. The lag values between samples, for each box plot correspond to a number of years (facet labels, from 0 to 7) plus a number of months (x axis of each facet, expressed as ranges). For example, the lag between samples considered for the first box plot is 0 years and 0 to 1 months and the lag between samples considered for the last box-plot in 7 years and 11 to 12 months. Panels (A) and (C) are based on the morphological dataset (cell counts) while graphs (B) and (D) are based on the metabarcoding dataset.

Figure 4
Low-taxonomic resolution contribution of protists at the SOMLIT-Astan time-series station over the period 2009-2016. Abundance of (A) the 12 main phytoplankton classes for the morphological dataset; and (B) the 52 main phyla -or classescalculated from the metabarcoding dataset. The tree maps show the overall contributions of the main phyla or classes to the (C) total species or (D) OTU richnesses.

Figure 5
Typical seasonal variations of the dominant OTUs and overall contribution of the major diatoms species to the protist assemblage at the SOMLIT-Astan sampling station over the period 2009-2016. The histograms show the contributions (A) to total DNA reads abundance of the 32 dominating OTUs (accounting for 51.5% of all reads), (B) of the main diatoms to total diatoms abundances (microscopy count of plankton >10um) and (C) of the main diatoms to total diatom reads abundances. All microscopy counts and OTUs were assigned at the highest taxonomic level. Species selected were the 10 most abundant (5 for diatoms) for at least one month, taking into account mean monthly abundances.  Monthly mean abundance (2009-2016), at the SOMLIT-Astan sampling site, for (A) the morphological species and (B) molecular OTUs as a function of the first three RDA axes (see Figure 6). For each RDA axis the (A) 10 species and (B) 10 OTUs with the highest score were selected. The * indicates dominant OTUs (reported in Fig.5).

Figure 8
Spearman correlation calculated between the environmental variables and the RDA axes (A, C), and variance partitioning analyses between environmental drivers and dbMEM (B, D). Spearman correlations were computed between each axes of the RDA and each environmental parameter selected for (A) morphology and (B) metabarcoding. Variance partitioning between selected environmental variables and dbMEM was also calculated for (C) morphology and (D) metabarcoding data, respectively. Table 1 Species and OTUs driving the seasonal oscillation showed in the RDA axes 1 and 2 and the biannual broad scale oscillations observed in axes 3. The 10 species/OTUs with the highest scores in each relative axes of the RDA were selected. The scores and the relative contributions to total abundance of the resulting list of species/OTUs are shown. For metabarcoding data, see Material and methods section and table S1 for assignations details.