The study area straddles the Strait of Dover and extends from the 50.5° N to 52.0° N, comprising the extreme north-eastern part of the English Channel and the Southern Bight of the North Sea (SBNS, Figure 1 ). The main water inflows include Atlantic water masses entering the English Channel in the southwest and passing through the Strait of Dover at high velocity (mean annual flow velocity = 0.16x10⁶m³ s^-1; Winther & Johannessen, 2006), and a current along the UK east coast with a NE-SW direction originating from the northern North Sea ( Figure 1 ).

The study area constitutes a temperate, shallow, permanently mixed system, presenting some zones of intermittent stratification such as the offshore north-eastern part of the Belgian sector and some parts of the offshore waters of the Dutch sector (De Boer et al., 2009; Van Leeuwen et al., 2015). The main freshwater inputs into the southern North Sea come from the Rhine and Meuse, and to a much lower extent from the Scheldt ( Figure 1 ), which constitute together the main freshwater discharge of the zone, and called the Rhine-Meuse-Scheldt system throughout this article. On the coast of England and the coast of France in the English Channel, the Thames ( Figure 1 ) and the Seine respectively constitute the main freshwater inputs (Desmit et al., 2018). For the Channel coast of France, the Somme estuary and smaller estuaries are responsible for a coastal Region of Freshwater Influence (ROFI) known as the “coastal river flow” (Brylinski et al., 1991).

As part of the H2020 JERICO-NEXT project (https://www.jerico-ri.eu/), an oceanographic campaign was carried out on board the RV “Simon Stevin” between 8^th and 12^th of May 2017, covering a short temporal window of the productive spring bloom season. It involved five research institutions, two from the Netherlands (RWS, NIOZ), two from Belgium (Ghent University (UGENT) and Flanders Marine Institute (VLIZ)) and one from France [Oceanology and Geosciences Laboratory (LOG)].

A total of 29 stations ( Figure 1 ), with a minimum and maximal distance (Euclidian) of 11 and 240 km respectively from each other, were sampled for environmental parameters (defined as abiotic factors; depth, water velocity, temperature, salinity, SPM and inorganic dissolved nutrients concentrations and ratios), bacteria, protists, pigments and zooplankton. Several sampling and analytical methods for characterizing the plankton compartment were applied, resulting in the consideration of five ecosystem components in addition to the abiotic factors: bacterial diversity analysed by amplicon sequencing, protist diversity (at the genus level) analysed with the FlowCAM, pico- to micro-plankton groups abundance analysed by flow cytometry (CytoSense), pigment concentrations analysed by HPLC and zooplankton taxonomic group abundances analysed with the ZooScan. The sampling and analytical methodologies are presented for each ecosystem component separately in the following sections.

At every station, discrete water samples were taken at 3 meters depth by means of six Niskin bottles (5 L) attached to a CTD-carrousel (Seabird SBE25plus). For each sample, 200 mL of sea water were then filtered through a 47 mm, 0.2 µm cellulose-acetate filter (Sigma-Aldrich). The filtrates were analysed using a SEAL QuAAtro analysis for concentrations of dissolved inorganic silica, phosphorous, ammonium, nitrite, and nitrate (expressed in µmol L^-1). Total nitrogen was calculated by computing the sum of ammonium, nitrite and nitrate. Full protocols are described in Mortelmans et al. (2019a). The SPM data were generated using a generic multi-sensor algorithm (Nechad et al., 2010) applied to the red remote sensing reflectance band (665 nm) available on the CMEMS server (OCEANCOLOUR_ATL _OPTICS_L3_REP_ OBSERVATIONS_009_066 ¹ ). Additional parameters originating from an underway data acquisition system were selected based on the time of arrival at each station, and include depth, temperature and salinity (SBE21 sensor). Associated current velocities were obtained based on the sdmpredictors R-package using mean surface current velocities (BO2_curvelmean_ss) from bio-ORACLE v2.0 >(Tyberghein et al., 2012; Assis et al., 2018).

A subsample of seawater from the sampler carousel (3 m depth) was filtered over a 25 mm 0.22 µm polycarbonate filter (mixed cellulose ester membrane GSWP filter, Merck) until saturation. After filtration the filter was stored in a 1.5 mL Eppendorf tube, snap frozen and transported to the laboratory in liquid nitrogen, where it was transferred for storage in the -80°C freezer until analysis. Extraction and isolation of genomic DNA included a beat-beating method with phenol-chloroform extraction based on the protocol of Zwart et al. (1998).

For each sample, bacterial amplicon libraries were constructed. For bacteria, the V1-V3 regions of the 16S SSU rRNA gene were amplified using the forward primer pA (5′-AGAGTTTGATCCTGGCTCAG-3′) (Edwards et al., 1989) and reverse primer BKL1 (5′-GTATTACCGCGGCTGCTGGCA-3′) (Cleenwerck et al., 2007). Amplifications were performed in duplicates with Polymerase Chain Reaction (PCR), using 2.5 µL PCR reaction buffer, 2.5 µL dNTP (2 mM) (Life technologies Inc.), 0.25 µL Fast Start High fidelity Taq polymerase (Roche Inc.), 2 µL of 16s and 18s SSU rRNA primer (0.25 μM) and 1 μL of extracted DNA. Sterilized HPLC grade water was added to obtain a final volume of 25 μL.

The PCR-program started with a DNA denaturation step of 96°C for 5 min, followed by 35 cycles of denaturation at 96°C for 1 min, annealing at 52-57°C for 1 min, extension at 72°C for 3 min, and a final elongation at 72°C for 20 min. The PCR products were purified with Agencourt AMPure XP beads (Beckman Coulter Inc.), duplicates were pooled after quality control with a BioAnalyzer (Agilent Inc.) and Qubit (Thermo Fisher Inc.) and the amplicon libraries were barcoded using the NEXTERA XT DNA kit (Illumina Inc.). High throughput sequencing was performed using a 300bp paired-end Illumina MiSeq machine (MiSeq, Edinburgh genomics). The forward and reverse reads were merged using Pear (Paired-End reAd merger) v.0.9.11. The UPARSE pipeline (Edgar, 2013) was used for dereplication, removal of singletons, removal of chimera’s and de‐novo clustering to the 97% similarity level and to transform the raw sequences to Operational Taxonomic Units (OTU). The Bayesian classifier of Mothur v.1.39.5 was used, with a cut-off of 80 to blast the bacteria against the SILVA Ribosomal Reference database (SILVA SSU Parc, version 123 (Yilmaz et al., 2014). Unclassified OTU’s were removed from the dataset. A threshold of 1x10^–5 was used to create a binary outlier variable to remove species with low variances. The Trimmed Mean of Mu;-values (TMM) was used to calculate standardization factors to normalise the data (Robinson et al., 2010; Robinson & Oshlack, 2010) and systematic variability (false positive) was removed. The bacterial groups obtained were expressed in standardized number of reads per OTU per sample and the whole component is referred to as “Bacteria” in the rest of the present paper.

The larger size fraction of the protists (>55 µm) was collected by filtering 50 L of surface water through a 55 µm Apstein net. The content of the cod end was retrieved and preserved in 2-5% final concentration acid Lugol solution and stored in dark conditions at 4°C. The samples were brought and stored at VLIZ marine station before being analysed with the FlowCAM VS-4 (Fluid Imaging Technologies), which combines the technology of flow cytometry, camera and microscopy (Álvarez et al., 2011). Before processing, the samples were sieved over a 300 µm mesh to remove large colonies (which clog the FlowCAM flow cell) and diluted when necessary. Thus, it is important to stress that the results will not comprise plankton colonies and individuals larger than 300 µm or protists smaller than 55 µm.

For each sample, 15 mL was then processed at a flow rate of 1.7 mL min^-1, using a FC300 flow cell and a 4X objective. The AutoImage mode captures 20 frames per second, imaging every particle between 70 to 300 Equivalent Spherical Diameter (ESD). The Region of Interest (ROI) per frame were semi-automated identified with the auto-classification tool of VisualSpreadsheet software, using the filters based on a reference library learning set created for the Belgian Part of the North Sea (BPNS) within the LifeWatch framework (Amadei Martínez et al., 2020).

The learning set consists of 26 libraries, each with their own filter. Then, the classification was manually validated to remove the errors of the automatic prediction (using books of Tomas, 1997; Kraberg et al., 2010 and Alfred Wegener Institute for Polar and Marine Research (AWI), 2020). Finally, abundances (cells per Liter) were calculated, dividing the total number of counts per taxon by the fluid volume imaged, the volume of water filtered and the dilution factor (Amadei Martínez et al., 2020). Only plankton genera were considered in the statistical analysis since the FlowCAM also provides counts for non-organic particles such as detritus (pieces of plants, plastics, etc.), air-bubbles and camera artefacts. In addition and in relation to the scope of the initial project, only autotrophic groups have been considered in the analysis. Thus, apart from two genera not strictly autotrophic considered (Protoperidinium and Tripos, being heterotrophic and mixotrophic respectively), the nano-and micro-heterotrophic groups are not part of the data-set. This FlowCam data-set is referred to as “Protist-FlowCAM” in the rest of the manuscript.

A CytoSense^® (Cytobuoy b.v., the Netherlands) automated flow cytometer (FCM) was connected to the RV Simon Stevin underway data acquisition system, which pumps sea water at 3 m depth, for in situ protists measurements. This device is a “pulse shape-recording” flow cytometer, which records the complete scatter and fluorescence pulse shape of each particle (1-800 µm) that passes the laser beams. Particles are pumped with a calibrated peristaltic pump to pass the laser beams in a laminar flow, ensuring a single-cell/particle analysis of the samples.

The flow cytometer is equipped with two lasers, a blue (488 nm – 50 mW solid-state laser (Coherent Inc.) and a red one (635 nm – 50 mW solid-state laser (Coherent Inc.)). To capture the scattered and fluorescent light, the flow cytometer records the forward scatter signal (FWS) through a PIN photodiode and the sideward scatter (SWS), fluorescence orange (FLO) (536-601 nm), fluorescence yellow (FLY) (601-668 nm) and fluorescence red (FLR) (668 – 734 nm) on PhotoMultiplier Tubes (PMT). A measurement protocol with a pump speed of 2.1 µL s^-1 was applied for six minutes, resulting in an analysed volume of on average 500 µL per sample (favouring small cells counting). The sensitivity of the PMT was set for SWS = 60, FLO = 80, FLY = 80, FLR1 (FLR 488 nm laser) = 95, FLR2 (FLR 635 nm laser) = 95. A trigger was set on the SWS (29 mV) to eliminate background noise and unwanted particles. To remove unwanted non-fluorescent particles an additional trigger was applied, respectively on maximum fluorescence red on the blue laser (FLR1 max 6).

Prior to the clustering analysis the generated dataset required pre-processing to remove signals from unwanted particles. All pre-processing was carried out with the Cytoclus 3 v3.7.4.14 software (Cytobuoy b.v., the Netherlands). Protist clusters were defined based on their light scattering and fluorescence properties using Easyclus version 1.28 (Easyclus^© v1.28, Thomas Rutten Projects, the Netherlands). Two different clustering tools were combined. The first one was the lasso tool, which used a training dataset to define polygons (lasso’s) around the Synechococcus and Cryptophytes clusters, thus defining the polygons to be used to cluster the rest. The defined selection sets of the lasso tool were then combined with the size fractionation (pico <3 µm; nano 3-20 µm; micro >20 µm), based on the length measurement obtained from the FWS signal, in the fixed clustering tool. Following this approach, a total of 5 clusters were defined: pico-red, nano-red, micro-red, pico-Synecho and nano-Crypto, which correspond to standardized flow cytometer cluster names being respectively eukaryote pico-phytoplankton, eukaryote nano-phytoplankton, micro-phytoplankton, Synechococcus and Cryptophytes (www.bodc.ac.uk) ( Table 1 ). The abundance of individual clusters was expressed in cells L^-1. This component is referred to as “Protist - FCM” in the rest of the paper.

The water samples for the pigment analysis were obtained from the same six Niskin bottles (5 L) deployed at 3 m depth and used for the nutrient analysis. Using a vacuum pump, the water was filtered on a 47 mm, 0.4 µm glass fibre filter (Whatman GF/F) up to saturation. Immediately after filtration, the filter was stored in a 1.5 mL Eppendorf tube, snap-frozen in liquid nitrogen and kept frozen at -80°C until analysis. Pigments were extracted in 90% HPLC grade acetone and sonicated for 30 seconds at 40 Hertz. After extraction, the extracts were immediately analysed by reverse phase HPLC following the protocol described by Van Heukelem and Thomas (2001), using an Agilent 1100 series HPLC system, with an Agilent Eclipse XDB-C8 column. Marker pigments of key phytoplankton groups were identified based on their retention time and absorption spectra and were validated using pure pigment standards (DHI Denmark). The final concentrations are expressed in µg per L. and we refer to this component as “Pigment”. The different pigments analysed can be seen listed in Table A in the Supplementary Material .

Zooplankton was sampled at each station with a 200 µm-mesh size WP-2 net, with a flow-meter attached to the frame, and deployed vertically from near-bottom to surface. Organisms collected in the cod-end were immediately preserved on-board with 7% buffered-formaldehyde and then stored, after the campaign, at the Marine Station Ostend, Belgium (MSO). Each sample was analysed with a ZooScan and processed with the Zooprocess software for semi-automated zooplankton identification (Gorsky et al., 2010). For more details on these specific steps, see Mortelmans et al. (2019b). The learning set used consists of 22 zooplankton taxa (at the class, phylum or order level) and one detritus group. These groups have been specifically created for the BPNS within the framework of the LifeWatch program (Mortelmans et al. 2019b). Noctiluca is a heterotrophic dinoflagellate of very large size, which is not adequately quantified by the method used for the protist-FlowCAM here. Since the ZooScan better quantifies the abundances of this group, Noctiluca has been considered within the zooplankton component in this paper. The final classification was manually validated to ensure the data quality. Zooplankton counts were then expressed as abundance (individuals per L.) using the sampled volume calculated from the flow meter data.

The different plankton compartments considered in this study and their sampling and analysis details are summarized in Table 2 .

The six ecosystem components (the abiotic component and the five biotic ones) were simultaneously analysed by focusing on the structural commonalities between the six “sampling stations × variables” data tables.

While some overlap exists in terms of size range amongst some biotic components, each of them represents a different and complementary ecological information justifying the inclusion of all components in a single analysis. The data set can be viewed as six sets of variables returning six spatial distributional patterns with their respective specificities, but also with potential commonalities.

In order to simultaneously handle the spatial variability of the six ecosystem components, the STATIS method was used (Structuration des Tableaux A Trois Indices de la Statistique; Abdi et al., 2012). STATIS is a multi-table ordination technique, which goes beyond traditional multi-variable analyses such as Principal Component Analysis (PCA) since it takes into account the importance of both individual variables and data tables on the multivariate axes (Thioulouse et al., 2018). Prior to the analysis, the data were appropriately transformed. Abiotic descriptors, given their different measurement units, were standardized (centred and reduced). Ranges of biotic descriptors strongly differed between tables so that non-zero values (presences) were rescaled between 1 and 5 within tables (intervals of 0.2 quantiles), after which they were centred.

The STATIS method proceeds in two main steps. Firstly, it builds a matrix of correlations between tables using the vectorial correlation RV (Robert & Escoufier, 1976), a multivariate equivalent of Pearson’s r-correlation coefficient. This matrix is then diagonalized to generate a system of axes called the “interstructure”, enabling the ordination of tables as in PCA, showing the degree of structural similarity between ecosystem components. Secondly, the first interstructure axis score, indicating the average correlation strength of each table with the others, is used to weight tables, in order to give them more or less importance in a second system of axes, called the “compromise”. Compromise axes are built from the sum of the six station-vector correlation matrices, maximizing the overall covariation between tables and providing synthetic scores of stations as an average spatial pattern associated to the projections of ecosystem component descriptors. Hence, the compromise axes were used to represent small- to larger-scale variations of the pelagic plankton system. A given ecosystem component may not be necessarily expressed on all compromise axes, and thus, we also considered the projections of the six separate PCA axes onto the STATIS compromise axes. The analysis was carried out in R 4.0.3 (R Core Team, 2020) with the package “ade4” (Chessel et al., 2004; Dray et al., 2007).

Figure 2 displays the spatial variations of the abiotic descriptors in the study area; complementarily, Table 3 provides statistical details on each descriptor. Chlorophyll a is also displayed as it is a core parameter measured in oceanographic studies. For conciseness, we will refer hereafter to “west coast” for the UK coast, and to “east coast” for the coast of continental Europe.

The depth pattern reflected the nearshore-offshore gradient. Current speeds were highest in the Strait of Dover, and tended to be overall lower in the shallower coastal zones. Salinity was lowest in the Belgian and Dutch coastal waters. Although the east coast was warmer, the thermal range of variation remained limited (10.3°C min – 12.2°C max). Nutrient concentrations and ratios were overall higher in shallow coastal waters, although high values were more homogeneously found in the western part, except for NH4 and SiO2, in higher concentrations along the eastern part. Whereas SPM concentrations decreased with depth according to a nearshore-offshore gradient, the pattern in Chlorophyll a was less clear, with lowest values isolated in the southwest and northeast zones of the study area.

Figure 3 displays the spatial variations of each biotic component total densities; complementarily, Table A in the Supplementary Material provides statistical details for the descriptors of each biotic component. Lowest bacterial densities were concentrated around the Strait of Dover and the west side, without much variations elsewhere ( Figure 3A ). Protist-FlowCam densities globally decreased from the nearshore to the offshore from both west and east sides ( Figure 3B ), whilst Protist-FCM densities exhibited opposite trends ( Figure 3C ). Pigment distribution ( Figure 3D ) was quite similar to Chlorophyll a distribution ( Figure 2M ), with lowest values isolated in the southwest and northeast zones. Zooplankton density exhibited the clearest pattern with highest densities concentrated along the east coast ( Figure 3E ).

The STATIS interstructure, i.e. the relationships between the ecosystem components, is displayed on Figure 4 and shows a general uni-dimensional pattern as all components were mostly correlated along the first axis. This indicated that the ecosystem components were not independent from each other and that they reflected a common spatial pattern. As displayed in Table 4 , correlations between ecosystem components were relatively homogeneous (weight around 0.40 for each component) with Bacteria and Pigment showing respectively the lowest weight (0.38 and 0.39) and Protist-FlowCAM the highest (0.45).

Three main axes emerged from the STATIS compromise, representing nearly 50% of the total variance of the data set ( Figure 5A ). The scores of the sampling stations for these three axes are mapped in Figures 5B–D .

This first axis of the STATIS compromise ( Figure 5B ), was spatially expressed as a longitudinal gradient opposing the western (coastal areas of the UK, high axis scores) to the eastern part (coastal areas of continental Europe, low axis scores) of the study area. The second axis of the STATIS compromise ( Figure 5C ) mainly reflected nearshore-offshore gradients with positive scores in the offshore zone. The third axis ( Figure 5D ), mainly opposed the southern from the northern parts.

The contributions of each ecosystem component to each compromise axis are displayed in Figure 6 . All components substantially contributed to the first axis (as seen on the correlation circles). Globally, the total ecosystem component covariation with the compromise was best for Protist-FlowCAM and Abiotic descriptors, unlike Bacteria and Zooplankton for which projections were limited to the first compromise axis (RV-coefficients, Figure 6 ), and in agreement with their respective contribution to the initial STATIS interstructure.

All ecosystem components contributed significantly to the first STATIS compromise axis (Figure 6). In comparison, bacteria, Protist-FCM and Zooplankton had a lower contribution to the second STATIS compromise ( Figures 6B,D, F ). The third axis resulted mostly from the covariation of Protist-FCM and Pigment ( Figures 6D, E ).

In Figure 6A , the PCA of Abiotic component returned two main axes that were mostly correlated with the two first STATIS axes (left correlation circle), whereas these two axes were weakly expressed on the third STATIS axis (right correlation circle, vertical axis). By contrast, Zooplankton ( Figure 6F ) shows a more complex PCA structure, expressed mostly on the first STATIS axis, the fourth Zooplankton axis (fourth eigenvalue) inducing little variance on the second STATIS axis; no Zooplankton axis was substantially expressed on the third STATIS axis.

Interplays between ecosystem component descriptors on the compromise axes are provided in Figure 7 for axes 1 and 2, and in Figure 8 for axes 2 and 3. In Supplementary Material, Table B provides correlations between descriptors and axes, and Supplementary Figures C-G and G provide the spatial distribution of each descriptor, respectively for Bacteria, Protist-FlowCAM, Protist-FCM, Pigment and Zooplankton.

The first axis was strongly explained by the nutrient concentrations and ratios (NO₂, NO₃, PO₄, N/P/Si and N/Si) ( Figures 7A and 8A ), consistently with the longitudinal trends visible in Figure 5B , which increased from the left to the right side of the axis (from the eastern to the western part, respectively). Among biotic components, size effects (multiple and positive covariances) were the most prominent features along this gradient. Most Protist-FlowCAM and Pigment descriptors (9 out of 16) covaried to the western part ( Figures 7C, E ); Bacillaria, Helicotheca and Pleurosigma/Gyrosigma were the most strongly characteristic taxa, followed by Ditylum, Odontella, Thalassiosira and other centric diatoms ( Figure 7C and Supplementary Figure D ); chlorophyll b, alloxanthin and beta-carotene, followed by diatoxanthin and diadoxinoxanthin were the most characteristic pigments ( Figure 7E ). Reversely, most Zooplankton descriptors positively covaried to the east coast, especially represented by Cumacea, Echinodermata and Branchiopoda, and with Harpaticoida more specific to the western part ( Figure 7F ). More symmetrically, Bacteria were represented in the western part by Gammaproteobacteria, and, to a lesser degree, Actinobacteria and Cyanobacteria; Deltaproteobacteria and Flavobacteria were the main representative of the eastern part ( Figure 7B ). Pico-Red, Pico-Syneco and Nano-Red were the dominant Protist-FCM clusters in the western part, opposed to Micro-Red, highly specific to the eastern part ( Figure 7D ). Nano-Crypto to a smaller extent was associated with the western part but only in the Strait of Dover.

The second axis, as a nearshore-offshore gradient ( Figure 5C ), showed increasing depth and current speed from the lower to the upper part of the axis ( Figure 7A ). Offshore stations (positive axis scores, Figure 5C ) included those located in the Strait of Dover and accurately aligned with the main SW-NE directed inflow of Atlantic water ( Figure 1 ). In addition, the negative axis scores were associated with high values of SiO₂ and NH₄ ⁺and SPM, nutrients, and less specifically high temperature, N/P ratio and nutrient contents (NO₂ and PO₄ ³; Figure 7A ) which corresponded to the characteristics of the stations in the coastal shallow areas ( Figure 2 ).

Most Protist-FlowCAM taxa, especially the pennate diatom Pseudo-nitzschia (due to its high abundance in the coastal BPNS; Supplementary Figure D ), and the centric diatoms Lauderia, Rhizosolenia, Thalassiosira and Lithodesmium ( Figure 7C ) were associated with the low axis scores and thus with the coastal areas. The dinoflagellates Protoperidinium and Tripos, heterotrophic and mixotrophic genus respectively, were the only taxa specific to the high axis scores, and thus to the offshore zone and the coastal stations in the Strait of Dover ( Figure 7C ). The main pigments associated with the coastal waters were chlorophyllide a, betacarotene, fucoxanthin and antheraxanthin in opposition to zeaxanthin characterizing the offshore zone ( Figure 7E ). The bacterial signature was limited to an offshore-nearshore opposition of Sphingobacteria and Gammaproteobacteria to Deltaproteobacteria ( Figure 7B and Supplementary Figure C ). Zooplankton was only expressed by its fourth axis ( Figure 6 ), and no clear trends could be distinguished on this second axis due to the absence of dominance in the offshore zone ( Supplementary Figure G ).

The third compromise axis ( Figure 8 ) was mostly limited to the expression of Protist-FCM and Pigment ( Figures 6D, E ; Figures 8D, E ). Within these components, limited positive covariances could be observed for nano-Crypto, 19-butanoyloxyfucoxanthin, chlorophyll c3, chlorophyllide a and diadinoxanthin ( Figures 8D, E and Supplementary Table B ) and seemed to create a slight south-north gradient along the axis.

The main regional gradient opposed the two main ROFIs and was strongly determined by both abiotic and biotic components, particularly nutrient concentrations and phytoplankton descriptors (Protist-FlowCAM, and Protist-FCM), with higher concentrations in the Thames ROFI. High nutrients and SPM concentrations were not limited to the coastal part near the Thames estuary, and actually extended to the stations offshore and NW of the estuary, most likely explained by the extent of the estuarine output, which forms an extensive plume as shown by satellite data for the study period (NASA Worldview, 2020 ² ). Given that abiotic descriptors are known to be determinant during the productive bloom season in the area (Lefebvre et al., 2011; Desmit et al., 2015; Desmit et al., 2019), particularly in relation to nutrient load, a reversed spatial pattern could have been expected, i.e. with higher concentrations in the Scheldt-Rhine-Meuse ROFI. Indeed, the nutrient load of the Scheldt-Rhine-Meuse system is on average 6 and 7 times higher in terms of N and P respectively than in the Thames (Pätsch et al., 2004; Desmit et al., 2019). In addition, the NE-SW current along the UK coast, originating from the northern North Sea and influencing the Thames ROFI, has not been reported to be nutrient rich (Blauw et al., 2012).

The different nutrient concentrations and ratios measured in the two ROFIs thus rather suggest variations in the uptake by the biotic compartment, in agreement with the strong relationships between nutrients and bloom productivity at this period (Reid et al., 1990; Nohe et al., 2020) and the potential variability in bloom timing and bloom species composition (Moschonas et al., 2017; Browning et al., 2020). On the western side, total Protist-FlowCAM concentrations were higher ( Supplementary Figure D ), particularly in the Thames ROFI, and was characterized by a mix of highly silicified taxa (e.g. Ditylum) and less silicified ones [Rhizosolenia (Rousseau et al., 2002), Bacillaria and Helicotheca (Kapinga & Gordon, 1992; Mann, 2006)], by specific diatom markers (diatoxanthin and diadoxinoxanthin), as well as by Gammaproteobacteria, which have been shown to be often associated to diatoms (Bidle & Azam, 2001; Klindworth et al., 2014; Wemheuer et al., 2014). In contrast, the east exhibited lower concentrations in protist descriptors, and the most characteristic ones, i.e. micro-red and to a lesser extent chlorophyll c3 and nano-red, tend to be used as indicators of the presence of Phaeocystis (Astoreca et al., 2009; Bonato et al., 2015; Li et al., 2021). This translates into different timing of the bloom period between the two ROFIs, in agreement with the large spatio-temporal variability characterisation of bloom season of the Southern part of the North Sea, which extends from February to October (Dulière et al., 2019).

In the western part, the descriptors indicated an exponential diatom bloom phase, as a result of silica uptake (Billen et al., 1991; Rousseau et al., 2002) and potential adaptation to low silicate levels (mainly under 2 µM), a yearly characteristic of this zone (Weston et al., 2008). Protists within the micro- and nano- size classes, given their concentrations and size, were most likely responsible for most of the bulk primary producers’ biomass in agreement with a clear bloom situation (Weston et al., 2008; Mills & Arrigo, 2010). A decrease in SPM, better light conditions and important dissolved nitrogen amounts are factors most likely explaining that picoeukaryotes were also present in higher concentrations in this zone (pico-red and pico-Synecho, Supplementary Figure E ) compared to the east side. Pico-plankton organisms have the ability to increase simultaneously with diatoms during diatom-dominated blooms (Barber & Hiscock, 2006) and their abundances largely outnumbered both the ones of the nano- and micro-phytoplankton groups (Protist-FCM data), a feature also confirmed for some stations of the same zone and period by Louchart et al. (2020). Low concentrations of their potential grazers, notably nano- and micro-dinoflagellates (Jeong et al., 2010), not measured in the present study, might have enhanced their presence, at least for the Thames ROFI since both nano- groups and pico-Synecho were in higher concentrations in the more southern UK coast. The bloom situation did not translate into the presence of mesozooplankton grazers in this west zone potentially due to the late diatom bloom timing rendering most mesozooplankton not able to meet their food requirement earlier in the season. Harpacticoid copepods were an exception, probably in relation to the ability of some species to withstand extreme conditions and undergo states of dormancy (Williams-Howze, 1997).

Late spring stage including the ending of a Phaeocystis globosa bloom would be the most straightforward explanation for nutrient depletion, with lower total Protist-FlowCAM concentrations in the east compared to the west (where only silicates were depleted) and despite continuous nutrient-rich river discharge in the Scheldt ROFI (also clearly indicated by the low salinities in this zone, cf. Figure 2). This is in agreement with the bloom onset occurring earlier than in the Thames ROFI (Gohin et al., 2019; Louchart et al., 2020), the latter having a peak generally in late May when the high SPM concentrations in this strong tidal mixing regime decrease (Blauw et al., 2018).

The Southern Bight of the North Sea, particularly along the coasts, is usually characterized by diatom blooms preceding, succeeding, or concomitant with a large bloom of Phaeocystis globosa that can represent 80% of the total bloom biomass (Breton et al., 2006; Lefebvre et al., 2011; Aardema et al., 2019). The presence of P. globosa suggested by our results (strong signature of micro-red, chlorophyll c3 and nano-red) was confirmed by the complementary study of Louchart et al. (2020) showing that the species was present in late spring 2017 all along the French coast by the Dover strait and the southern North sea. Most of the time, Phaeocytis globosa might have co-occurred with lower concentrations in diatoms compared to the west, with the exception of the coastal stations of the BPNS ( Supplementary Figure E ). The BPNS was the only zone in the east to be low in chlorophyll c3 ( Supplementary Figure F ) while high in micro-red concentrations related to a bloom of Pseudo-nitzschia (representing >70% of total Protist-FlowCAM abundance in the BPNS). This pennate diatom is characteristic of the intermediate assemblage and occurs before the summer season in this zone (Nohe et al., 2020). Correspondingly, the concentrations of mesozooplankton were particularly higher on the Belgian coast, potentially indicating grazing on this nano-colonial Pseudo-nitzschia (Harðardóttir et al., 2015; Hoffmeyer et al., 2020).

Interestingly, the higher mesozooplankton concentrations along the east coast do not only indicate that the earlier bloom situation was more adequate for the development of important contributor species of the meso-zooplankton communities compared to the west (Mackas et al., 2012; Mortelmans et al., 2021), but also indicates that mesozooplankton may benefit from grazing on P. globosa. In the most north-eastern side of the study area, the high characterisation by Cumacea ( Supplementary Figure G ) also corresponded to a strong characterisation of the zone by the indicators of P. globosa. This is in agreement with the results of Dauvin et al. (2008) which showed that Cumacea, a suprabenthic mesozooplankton group (feeding temporarily in the water-column), significantly increased its vertical migration into the water column during a P. globosa bloom in the eastern English Channel.

The second spatial gradient depicted through the analysis corresponded to the expected nearshore-offshore trend. It directly supports the implication of Atlantic influence that strengthens the ecological contrasts with the coast, except for the coastal part of the Strait of Dover to some extent, where the Atlantic influence generally dominates over local coastal processes (Lee, 1980; Huthnance, 1991; Dulière et al., 2019). In addition, and as expected, the Atlantic influx of water from the English Channel was characterized by low SPM and low concentrations of nutrients (Ruddick & Lacroix, 2006; Desmit et al., 2015; Dulière et al., 2019).

Along the coast, the continental input could be clearly seen through high concentrations of autotrophic protists within the nano- and micro-size classes (Protist-FlowCAM, and to a lower extent nano-red from the FCM), represented by centric diatoms, with the exception of the bloom of the diatom Pseudo-nitzschia in the BPNS. This was also confirmed by the fucoxanthin signature (Claustre, 1994; Jeffrey & Vesk, 1997) in both the east and the west coastal part of the study area ( Supplementary Figure F ). Diatoms were more characteristic of the west coast, as seen previously, and within the centric diatoms, only Rhizosolenia was present in both coastal areas. Except in the BPNS, this taxon dominated the relative abundance of the Protist-FlowCam in the coastal stations, which has been shown to be typical of the late bloom stage and/or summer assemblage (Weston et al., 2008; Nohe et al., 2020). Another coastal feature was the strong characterisation by chlorophyllide a, exhibiting higher concentration at a few stations in the BPNS and in the most north-eastern part ( Supplementary Figure F ). While chlorophyllide a has been usually considered as a physiological marker of chlorophyllase-containing diatoms (Barrett & Jeffrey, 1971; Bidigare, 1989), a recent study (Helmann, 2019) suggested that it should be used as an indicator of senescent, physiologically compromised phytoplankton due to stress. The bloom of potentially toxin-producing Pseudo-nitzschia (Trainer et al., 2012; Xu et al., 2015) in the BPNS and the presence of Phaeocystis colonies in the north-east could have been thus causing stress for the other phytoplankton species present.

In agreement with the onshore-offshore gradient in salinity, SPM and nutrients, the offshore stations were, in contrast with the coastal part, poorly characterized by the large protist size-class (Protist-FlowCAM, Figure 3 and Supplementary Figure D ). They were however highly characterized by higher concentrations of large heterotrophic/mixotrophic dinoflagellate taxa, Protoperidinium and Tripos. Dinoflagellates, in particular Protoperidinium, have been shown to be able to cope with highly turbulent systems (Hernández-Fariñas et al., 2014; Smayda et al., 2010), and can have a high grazing rate on autotrophic phytoplankton (Gribble et al., 2007; Grattepanche et al., 2011). While no important concentrations of autotrophic large cells were present in the offshore zone, dinoflagellates could have potentially taken advantage of the nano-plankton (nano-red), present in higher concentrations in most offshore stations, and of nano-Crypto and pico-plankton (both pico-red and pico-Synecho, cf. Supplementary Figure E ) in the southwestern offshore part corresponding mainly to the Strait of Dover. Part of the nano-red might actually correspond to the diatoms of the nano size class found with the FlowCAM (i.e. Chaetoceros, Navicula, centric diatoms; Supplementary Figure D ). In addition, cyanobacteria, indicated by zeaxanthin and matching the spatial distribution of pico-Synecho for the whole Strait of Dover, appeared as an important contributor to the offshore plankton pattern. These results corroborates with nano-and pico-size plankton being respectively an important proportion of the protist assemblage in the Atlantic offshore water and in oligotrophic waters more generally (Gibb et al., 2000; Marañón et al., 2003; Kostadinov et al., 2010; Marie et al., 2010; Masquelier et al., 2011). In spite of high nano- and pico-plankton concentrations in several offshore stations, the Protist-FCM component was not strongly expressed along this nearshore-offshore gradient (second STATIS axis, Figure 5c) due to a strong contribution of the pico-plankton groups to the east-west gradient (first axis, Figure 5b) but also to the differential spatial distribution of nano-Crypto according to a south-north opposition.

The south-north gradient drove a substantial part of the total spatial variation (STATIS axis 3, Figures 5a,d), with a strong characterisation of the south part by nano-Crypto which however did not match higher concentrations of its pigment marker, alloxanthin (Ansotegui et al., 2001). The nano-Crypto cluster, in this case, might has been rather composed of colonies of Synechococcus and Cyanobacteria, as supported by the distribution of their pigment marker zeaxanthin, and as they are known to be present in this zone ( Supplementary Figure F ; Mackey et al., 1996; Bonato et al., 2015). Some inconsistencies were also detected for the strong signature of the carotenoid 19-butanoyloxyfucoxanthin, defined as a clear marker of pelagophytes within the pico-size class (Irigoien et al., 2004; Raven, 2012), while its distribution did not exactly match the one of the pico-plankton groups. On the other hand, the characterisation by chlorophyll c3 for this south part, more particularly on the coast of France and in the north-easternmost zone, matched the one of the stress related pigment chlorophyllide a, most likely indicating the characteristic presence of Phaeocystis at this period (Astoreca et al., 2009; Bonato et al., 2015; Louchart et al., 2020; Li et al., 2021).

This study showed that the plankton spatial distribution during the productive spring period in the study area did not follow a unique pattern but was characterized by an overriding imprint of both local and regional drivers. High spatial heterogeneity in plankton distribution has been previously highlighted in this study area (Blauw et al., 2018; Aardema et al., 2019). However, while a nearshore-offshore gradient under Atlantic influence could be expected, spring bloom timing appeared to be a dominant driver of plankton community structure differing between the two main ROFI of the study area, bringing complementary insights on this spatial heterogeneity. By simultaneously considering several components of the plankton, our results also point to potential interactions within the plankton realm that are rarely assessed, especially at spatial scales beyond a unique sampling station.

However, while this study included five different plankton components, some plankton size-classes and groups are still missing from the analysis, limiting the extent of the data interpretation. For instance, we could infer the presence of Phaeocystis from proxies (micro-red, chlorophyll c3 and nano-red), but it would have been more adequate to assert its presence through microscopic counts, especially to provide quantitative data for such important protist biomass contributor. Blooms of colonial large protists such as P. globosa are often problematic for semi-automatic imaging devices (clogging effect) and thus, complementary classical techniques based on microscopy should complement the data analysis to allow a full coverage of the nano- and micro-protist size-classes. This is in general true for most plankton size-classes since the calibration settings and the volume of water sampled considered for semi-autonomous devices create some limitations (as an example, the calibration used for the flow cytometer, in this case study, favoured small organisms with the use of small sample volumes of 500 µL). For instance, taxonomic information for the nano-size classes complementing the data from the flow cytometer would have allowed us to better discuss potential inter-plankton relationships helping the interpretation of axis 3 and discrepancies between some descriptors. A full coverage of the plankton size-spectra would be also beneficial to address the challenging exercise in disentangling the complex interaction of abiotic and biotic factors driving the plankton distribution (Lima-Mendez et al., 2015; Bunse et al., 2016; Striebel et al., 2016). Nevertheless, in spite of missing groups, our results showed that usually neglected small components can have a preponderant structuring role as evidenced by the dominant weight of Protist-FlowCam in the STATIS analysis.

Despite some limitations, our joint analysis of abiotic factors and several plankton compartments enabled us to highlight interesting issues, particularly in the frame of monitoring. Chlorophyll a, often used as a proxy for phytoplankton biomass (Garmendia et al., 2013; Harvey et al., 2015), did not appear as an important factor contributing to the main spatial plankton distribution pattern during this period of intense production for instance. The use of chlorophyll a alone could have thus not enabled to highlight all the regional differences revealed by our analysis, with the exception of the expected nearshore-offshore gradient in protist biomass. While P. globosa is also generally considered detrimental for zooplankton grazers (Weisse et al., 1994; Rousseau et al., 2000; Nejstgaard et al., 2008), the present results suggests a positive interaction with higher planktonic trophic levels. As such, and more generally, our observations underscore that mere monitoring of bulk phytoplankton biomass is often insufficient and even inaccurate to inform about the plankton status. This strengthens the general requirement for a multi-faceted plankton monitoring which includes information on the plankton composition and which consider complementary techniques including the classical ones based on microscopy (Alvarez-Fernandez & Riegman, 2014; Aubert et al., 2017; McQuatters-Gollop et al., 2017).

Finally, only a snapshot of the seasonal cycle was considered here, and the question of temporal variability remains to be explored. The smallest planktonic size-classes (bacteria to nano-plankton) have different population dynamics compared to larger plankton (micro-phytoplankton and zooplankton) and are often neglected in monitoring programs, particularly in the frame of management strategies (Aubert et al., 2017). As they clearly constitute an important component of the plankton as demonstrated by our study, their integration in plankton studies could improve our understanding of ecosystem dynamics (Aubert et al., 2017; Chust et al., 2017; Lombard et al., 2019), especially in areas of high spatio-temporal variations such as the Southern Bight (Blauw et al., 2018; Ivanov et al., 2020). A more in depth understanding of the whole plankton dynamics, considering plankton inter-relationships, will definitely help to improve our predictability analysis in areas such the Southern Bight known to experience important bloom dynamics shifts (Aardema et al., 2019; Nohe et al., 2020; Mortelmans et al., 2021).