The contrasted phytoplankton dynamics across a frontal system in the southwestern Mediterranean Sea

Numerical simulations have shown that finescale structures such as fronts are often suitable places for the generation of vertical velocities, transporting subsurface nutrients to the euphotic zone and thus modulating phytoplankton abundance and community structure. In these structures, direct in situ estimations of the phytoplankton growth rates are rare; although difficult to obtain, they provide precious information on the ecosystem functioning. Here, we consider the case of a front separating two water masses characterized by several phytoplankton groups with different abundances in the southwestern Mediterranean Sea. In order to estimate possible differences in growth rates, we measured the phytoplankton diurnal cycle in these two water masses as identified by an adaptive and Lagrangian sampling strategy. A size-structured population model was then applied to these data to estimate the growth and loss rates for each phytoplankton group identified by flow cytometry, showing that these two population parameters are significantly different on the two sides of the front and consistent with the relative abundances. Our results introduce a general method for estimating growth rates at frontal systems, paving the way for in situ exploration of finescale biophysical interactions.

1 Introduction

Phytoplankton forms the basis of the marine food web (Sterner and Hessen, 1994), and it is responsible for half of the primary production of the planet (Field et al., 1998), while its biomass is only ≤ 1 % of the global biomass (Winder and Cloern, 2010). Thanks to photosynthesis, phytoplankton fuels the ocean and the atmosphere in free O₂, and it fixes and exports the CO₂ into the ocean depth (Field et al., 1998; De La Rocha and Passow, 2007). This process called biological carbon pump is critical for global ocean sequestration of carbon and therefore for the modulation of atmospheric CO₂. The biological carbon pump is modulated by the size structure of the phytoplankton community. Small or large phytoplankton species are associated with different efficiencies for particle export, remineralization and transfer to the deep ocean (Boyd and Newton, 1999; Guidi et al., 2009; Hilligsøe et al., 2011; Mouw et al., 2016, etc). That is why it is fundamental to understand the factors that rule phytoplankton abundance and diversity.

Here, we use the term “finescale” to refer to ocean dynamical processes induced by mesoscale interactions and frontogenesis (Capet et al., 2008b, a; McWilliams, 2016; McGillicuddy, 2016; Lévy et al., 2018). Finescale structures are characterized by a small Rossby number, a horizontal scale of the order of 1–100 km and a short lifetime (days–weeks). Numerical simulations and remote sensing observations have demonstrated that finescale lifetime is often similar to the phytoplankton growth timescale, suggesting that finescale processes can affect and modulate the phytoplankton community. Different physical processes associated with finescale structures are able to generate vertical velocities, such as deformations of the flow and spatial inhomogeneities (Giordani et al., 2006), eddy perturbation (Martin and Richards, 2001; Pilo et al., 2018), linear Ekman pumping (McGillicuddy et al., 1998; Gaube et al., 2015), or eddy–wind interactions (McGillicuddy et al., 2007). Previous studies have well established that vertical motions impact biogeochemistry (Mahadevan and Tandon, 2006; Mahadevan, 2016; McGillicuddy, 2016). Upward vertical velocities drive deep nutrients into the euphotic layer and also move the phytoplankton cells along the water column, resulting in changing light conditions. However, most of the in situ studies related to the physical–biological coupling at finescale have focused on extreme situations occurring in coastal upwelling regions (Ribalet et al., 2010) or in boundary currents (Clayton et al., 2014, 2017), where intense fronts and dramatic contrasts in water properties are found that are, however, not representative of the global ocean. Indeed, vast oceanic regions are dominated by weak fronts continuously created, moved and dissipated, which separate different water masses with similar properties. The ephemeral nature of these finescale structures makes them particularly challenging to sample. As a consequence, some recent cruises have used remote sensing and numerical simulations to define the sampling strategy, allowing one to target and measure finescale features with physical sensors at high frequency (Shcherbina et al., 2015; Pascual et al., 2017; Petrenko et al., 2017). Concerning the biological variables, although progress in the understanding of the phytoplankton cell cycle has been obtained from incubation, sample manipulation (Worden and Binder, 2003) and models (Geider et al., 1997; MacIntyre et al., 2000), performing in situ measurements at high frequency and resolution is a necessity to better understand these biological processes and their responses to the environment.

An efficient solution is to lead Lagrangian cruises using automated flow cytometers sampling at high frequency in order to resolve the phytoplankton diurnal cycle in situ, which is challenging using more conventional methods such as cultures or counting by optical microscopy (Thyssen et al., 2008; Fontana et al., 2018). This is the solution chosen by the PROTEVSMED-SWOT cruise. This cruise was performed in the southwestern Mediterranean Sea, south of the Balearic Islands (Dumas, 2018; Garreau et al., 2020), with the aim to study the physical and biological coupling at finescale. This area is characterized by the presence of both fresh surface waters coming from the Atlantic (AW) and more saline waters from the Mediterranean region (Millot, 1999; Millot et al., 2006). AW enters the Mediterranean Sea through the Strait of Gibraltar and then forms a counterclockwise circulation along the continental slope of the western Mediterranean basin, caused by the combination of the Coriolis effect and the topographical forcing (Millot, 1999; Millot and Taupier-Letage, 2005; Millot et al., 2006). In the southwest part of the basin, this circulation is dominated by the Algerian Current (AC), which can form meanders and mesoscale eddies due to baroclinic and barotropic instabilities (Millot, 1999). These eddies spread over the basin and join the study area south of the Balearic Islands, carrying with them the newly arrived AW, known as younger AW. In this region, the younger AW encounters the older AW, sometimes also called resident AW (Balbín et al., 2012) or local AW (Barceló-Llull et al., 2019). The older AW is AW modified by cooling and evaporation during its progression along the northern part of the western Mediterranean basin. The encounter between these two AWs often generates finescale frontal structures (Balbín et al., 2014). To our knowledge, except for the works of Balbín et al. (2012, 2014) and the glider experiments of Cotroneo et al. (2016) and Barceló-Llull et al. (2019), very few studies have been performed in this region, and these frontal finescale structures have scarcely been sampled due to the difficulty of performing in situ experiments over these short-lived and small features. During the PROTEVSMED-SWOT cruise, a Lagrangian adaptive sampling strategy was performed across a moderately energetic front separating two distinct AWs at different stages of mixing (Tzortzis et al., 2021). The AW located south of the front is characterized by absolute salinity (S_A) between 37 and 37.5 g kg⁻¹, corresponding to the younger AW recently entered into the Mediterranean Sea. However, north of the front, the AW referred to as the older AW is characterized by a higher S_A (37.5 to 38 g kg⁻¹). Tzortzis et al. (2021) have also observed contrasted phytoplankton abundances in these two water masses, with the smallest phytoplankton such as Synechococcus dominating south of the front in the younger AW, while microplankton was more abundant north of the front in the older AW. As a consequence, our previous study constitutes an important improvement in the understanding of the role of frontal structures at finescale on phytoplankton distribution in a moderately energetic ocean. Nevertheless, open questions remain concerning the mechanisms generating this observed distribution. Is it exclusively driven by the dynamics of the ocean currents? What is the role of biological processes? In the present study, we attempt to shed some light on these questions and in particular on the patterns of phytoplankton abundances observed by automated flow cytometry during the PROTEVSMED-SWOT in the frontal structure, using the size-structured population model of Sosik et al. (2003). The specific objective of our study is to assess whether the observed contrasted abundances across the front were due to different growth and loss rates. Using high-frequency flow cytometry measurements across the front dividing two water masses, we were able to separately analyze each phytoplankton functional group and reconstruct their biovolume dynamics over a diel cycle in each water mass.

Figure 1 (a) Route of the R/V Beautemps-Beaupré during the PROTEVSMED-SWOT cruise. The purple box encloses a (b) zoom of the sampling region with overlaid chlorophyll a concentration (µg L⁻¹) of 11 May 2018. In panel (b) the dotted black line represents the route of the ship and the bold black line represents the route of the Lagrangian sampling across the older AW (delimited by the box in dark blue) and the younger AW (delimited by the box in light blue). (c) Dates of the transects across the older AW and the younger AW, used to reconstruct a 24 h period in each water mass.

2 Materials and methods

2.1 The sampling strategy

The PROTEVSMED-SWOT cruise, dedicated to the study of finescale dynamics, was conducted in the south of the Balearic Islands between 30 April and 18 May 2018 on board the R/V Beautemps-Beaupré (Fig. 1a). This cruise followed an adaptive Lagrangian strategy to measure at high spatial and temporal resolution several physical and biological variables with both in situ sensors and analysis of the sea surface water intake. The vessel route was designed ad hoc on the basis of a daily remote sensing dataset provided by the Software Package for an Adaptive Satellite-based Sampling for Oceanographic cruises (SPASSO, https://spasso.mio.osupytheas.fr, last access: 17 July 2023). SPASSO used altimetry-derived currents from the Mediterranean regional product (nrt_med_allsat_phy_l4) AVISO (Archiving, Validation and Interpretation of Satellite Oceanographic data, https://www.aviso.altimetry.fr, last access: 17 July 2023) and ocean color observations. Chlorophyll a concentrations ([Chl a], level 3, 1 km resolution, MODIS Aqua and National Polar-orbiting Partnership (NPP) Visible Infrared Imaging Radiometer Suite (VIIRS) sensors combined (after 27 May 2017) into a new product called MULTI) were provided by CMEMS (Copernicus Marine Environment Monitoring Service, https://marine.copernicus.eu, last access: 17 July 2023). In addition, Collecte Localisation Satellites (CLS) provided the surface Chl concentration composite products, with the support of the Centre National d'Etudes Spatiales (CNES). They were constructed using a simple weighted average over the previous 5 d of data gathered by the Suomi NPP VIIRS sensor. SPASSO generated maps of dynamical and biogeochemical structures in both near real time (NRT) and delayed time (DT). Maps of [Chl a] allowed us to identify two water masses, characterized by distinct [Chl a] values and separated by a zonal front at around 38^∘30^′ N. This front was also detected using in situ horizontal velocities, temperature and salinity, as described in Tzortzis et al. (2021). These two water masses were sampled along a designated route of the ship, represented in black in Fig. 1b. Special attention was paid to adapting the temporal sampling in order to measure the phytoplankton diel cycle in each water mass. This was achieved by continuously sampling across both water masses along the transects. While the ship did not remain in each water mass for 24 h, day-to-day variability remained low, and measurements from several days were combined into one diel cycle (Fig. 1c). The shape depicted by the ship's track led us to call these areas the north–south (NS) hippodrome (bold black line in Fig. 1b), performed between 11 and 13 May 2018.

2.2 In situ measurements

During the cruise, the irradiance (wavelengths between 400 and 1000 nm) was measured by a CMP6 pyranometer (Kipp and Zonen; https://www.campbellsci.fr/cmp6, last access: 17 July 2023). Temperature and salinity were measured by a thermosalinograph (TSG). The TSG was equipped with two sensors: a CTD Sea-Bird Electronics (SBE) 45 sensor installed in the wet lab, connected to the surface water and which continuously pumped seawater at 3 m depth, and a SBE 38 temperature sensor installed at the entry of the water intake. The TSG measurements were taken every 30 min, which corresponds to around 2 km spatial resolution at typical ship speeds. The data were converted into conservative temperature (Θ) and absolute salinity (S_A) using the TEOS-10 standards of McDougall et al. (2012). To automatically sample and analyze phytoplankton cells, an automated CytoSense flow cytometer (CytoBuoy, b.v.; (Dubelaar et al., 1999; Dubelaar and Gerritzen, 2000)) was installed on board and connected to the seawater circuit of the TSG. The flow cytometer sampled the seawater in a dedicated small container called a “subsampler”. The subsampler isolates the seawater every 30 min, which allowed us to ignore the movement of the ship while the flow cytometer performed its analysis. Between two consecutive samples the subsampler was flushed continuously by the seawater circuit of the ship in order to clean and renew the seawater. A sheath fluid made of 0.1 µm filtered seawater stretched the sample in order to separate, align, center and drive the individual particles (i.e., cells) through a laser beam (488 nm wavelength). Several optical signals were recorded when each particle crossed the laser beam: the forward angle light scatter (FWS) and 90^∘ sideward angle scatter (SWS), related to the size and the structure (granularity) of the particles. Two distinct fluorescence emissions induced by the light excitation were also recorded, a red fluorescence (FLR) induced by chlorophyll a content and an orange fluorescence (FLO) induced by the phycoerythrin pigment content. The CytoUSB software (CytoBuoy b.v.) was used to configure and control the flow cytometer and set two distinct protocols. The first protocol (FLR6) was dedicated to the analysis of the smaller phytoplankton using a red fluorescence (FLR) trigger threshold fixed at 6 mV and an analyzed volume set up at 1.5 mL. The second protocol (FLR25) targeted nanophytoplankton and microphytoplankton with an FLR trigger level fixed at 25 mV and an analyzed volume of 4 mL. The FLR trigger was used to discriminate the red fluorescing phytoplanktonic cells from other particles (such as heterotrophic prokaryotes, nanoflagellates, ciliates, etc.). Recorded data were analyzed with the CytoClus software (CytoBuoy b.v.) which retrieves information from the four pulse shape curves (FWS, SWS, FLO, FLR) obtained for every single cell. These curves were then projected into distinct two-dimensional planes (cytograms) by computing the curves' integrals. Using a combination of various cytograms (e.g., FWS vs. FLR, FLO vs. FLR) allows us to determine optimal cell clusters, i.e, cells sharing similar optical properties (see Fig. A1). The identification of phytoplankton functional groups is described in the Appendix. These clusters have been demonstrated in the literature to represent phytoplankton functional groups (PFGs) (Dubelaar and Jonker, 2000; Reynolds, 2006; Thyssen et al., 2008; Edwards et al., 2015; Thyssen et al., 2022). Finally, the PFG abundance (cells per milliliter), mean light scatter and fluorescence intensities were extracted from each sample.

2.3 The size-structured population model

We used the size-structured population model described by Sosik et al. (2003) and adapted by Dugenne et al. (2014) and Marrec et al. (2018) to estimate the in situ growth rates of every phytoplankton group identified by the CytoSense flow cytometer in the older AW and the younger AW. Before applying the model, we reconstructed a daily cycle of 24 h in the two water masses for each phytoplankton group. We use the term reconstruction because the ship did not spend 24 h in a row in each water mass but sailed along two routes, each forming a sort of racetrack passing alternately through the two water masses (Fig. 1b and c). By eliminating the dates and keeping the associated sampling times, the 24 h diel cycle can be reconstructed for each water body (Fig. 1c). This relies on the hypothesis that the phytoplankton community and dynamics remained similar over the 2 d and that the hydrology and physics for each water mass remained alike during sampling. We also reconstructed the 24 h irradiance in the two water masses (Fig. A2) because one of the most important parameters of this model is irradiance, since cell growth is dependent on light exposure due to photosynthesis.

The model of Sosik et al. (2003) uses the phytoplankton cell volume (biovolume) derived from cell light scatter intensities (FWS) as input (Eq. 1). Biovolumes were estimated using coefficients previously obtained by measuring a set of silica beads with the flow cytometer following the same settings used for phytoplankton analysis. The coefficients β₀ and β₁ used to convert FWS (arbitrary units, a.u.) to biovolume v (µm³) were derived from a log–log regression between FWS and silica bead volumes. These methods come from the studies of Koch et al. (1996) and Foladori et al. (2008).

$$\begin{array}{}\text{(1)}& v=\text{exp}\left({\mathit{\beta }}_{\mathrm{0}}\right)\times {\mathrm{FWS}}^{{\mathit{\beta }}_{\mathrm{1}}},\end{array}$$

with, in our case, β₁ = 0.9228 and β₀ = −5.8702.

In the size-structured population model, cells are classified into several size classes according to their dimensions at time t. Classes are logarithmically spaced as follows: for i in $\mathrm{1},\mathrm{2},\mathrm{\dots },m$ $v_i=v_{\mathrm{1}}{\mathrm{2}}^{(i-\mathrm{1})\mathrm{\Delta }v}$, where Δv is constant and chosen to ensure that size classes cover the entire observed biovolume v, from v₁ to v_m (Fig. 2). For Synechococcus, Δv = $\mathrm{1}/\mathrm{6}$ with Δv constant and m = 40 so that the model size classes encompassed our full measured size distributions (0.0279–2.5209 µm).

Figure 2 Cell cycle stages in the size-structured population model. Cells may grow to the next size class (γ) or be at equilibrium $(\mathrm{1}-\mathit{\gamma }(t\left)\right)(\mathrm{1}-\mathit{\delta }(v,t\left)\right)$. Above a particular size, cells are large enough to divide into two daughter cells with probability (δ). Figure adapted from Sosik et al. (2003).

At any time t, the number of cells in size classes N (and w its corresponding normalized distribution) was projected to t+dt via matrix multiplication (Eq. 2):

$$\begin{array}{}\text{(2)}& N(t+\mathrm{dt})=\mathbf{A}\left(t\right)N\left(t\right)\text{ and }w(t+\mathrm{dt})={\displaystyle \frac{\mathbf{A}\left(t\right)N\left(t\right)}{\sum \mathbf{A}\left(t\right)N\left(t\right)}}.\end{array}$$

We chose dt = 10 min (i.e., $\mathrm{10}/\mathrm{60}$ h), as in Sosik et al. (2003) and Dugenne et al. (2014), because cells of a specific phytoplankton group are unlikely to grow more than one size class over such a small time duration.

A(t) is a tridiagonal transition matrix that contains

1. γ: the probability of cellular growth

2. δ: the probability of cells entering mitosis

3. the cell stasis, i.e., the probability for cells to maintain their state (i.e size) in equilibrium during the temporal projection.

Probability of cellular growth

The probability of cells growing to the next size class (γ) depends only on the light intensity (irradiance) necessary for photosynthesis, expressed as (Eq. 3)

$$\begin{array}{}\text{(3)}& \mathit{\gamma }\left(t\right)={\mathit{\gamma }}_{\text{max}}\cdot (\mathrm{1}-\text{exp}(-E\left(t\right)/E^{\ast }\left)\right),\end{array}$$

with γ_max representing the maximum proportion of cells growing (dimensionless quantity), E the irradiance ($\mathrm{\mu }\mathrm{E}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and E^∗ the irradiance normalizing constant ($\mathrm{\mu }\mathrm{E}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$).

Probability of cells entering mitosis

According to Dugenne et al. (2014), δ expresses a proportion (between 0 and 1) modeled by the combination of two normal distributions (𝒩). One is linked to the cell size, the other is linked to the time of cell division. Both imply an optimum, reached at $\stackrel{\mathrm{‾}}{v}$ and $\stackrel{\mathrm{‾}}{t}$, respectively, for cell division, above which the cell size and the timing of division is suboptimal (Eq. 4).

$$\begin{array}{}\text{(4)}& \mathit{\delta }(t,v)={\mathit{\delta }}_{\text{max}}\mathcal{N}(\stackrel{\mathrm{‾}}{v},{\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}})\mathcal{N}(\stackrel{\mathrm{‾}}{t},{\mathit{\sigma }}_{\mathrm{t}}^{\mathrm{2}}),\end{array}$$

where δ_max represents the maximum proportion of cells entering mitosis (dimensionless quantity), $\stackrel{\mathrm{‾}}{v}$ the mean of the size normal distribution (µm³), σ_v the standard deviation of the size normal distribution (µm³), $\stackrel{\mathrm{‾}}{t}$ the mean of the time normal distribution (h) and σ_t the standard deviation of the time normal distribution (h).

Cell stasis

A third functional proportion is included in the transition matrix A(t) to represent cell stasis. Since this function illustrates a non-transition, it is modeled by the proportion of cells that neither divided nor grew between t and t+dt (Eq. 5).

$$\begin{array}{}\text{(5)}& [\mathrm{1}-\mathit{\gamma }(t\left)\right][\mathrm{1}-\mathit{\delta }(t,v\left)\right]\end{array}$$

Table 1 Model parameters being optimized.

Optimal parameters

In the set of parameters, θ is estimated by the maximum likelihood function, assuming errors between observed w and predicted $\widehat{w}$ normalized size distributions (Eqs. 6–8). Their standard deviations are estimated by a Markov Chain Monte Carlo approach (Geyer, 1992; Neal, 1993) that samples θ from their prior density distribution, obtained after running 200 optimizations on bootstrapped residuals to approximate the parameter posterior distribution using the normal likelihood (the likelihood function represents the probability of random variable realizations conditional on particular values of the statistical parameters).

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathit{\theta }=[{\mathit{\gamma }}_{\text{max}},E^{\ast },{\mathit{\delta }}_{\text{max}},\stackrel{\mathrm{‾}}{v},{\mathit{\sigma }}_{\mathrm{v}},\stackrel{\mathrm{‾}}{t},{\mathit{\sigma }}_{\mathrm{t}}]=\text{argmin}\left(\sum \left(\mathit{\theta }\right)\right),\text{(7)}& {\displaystyle }\sum \left(\mathit{\theta }\right)=\sum _{\mathrm{t}}^{t+\mathrm{dt}}\sum _{i=\mathrm{1}}^m\left(w\right(t)-\widehat{w}(t,\mathit{\theta }){)}^{\mathrm{2}},\text{(8)}& {\displaystyle }\widehat{N}(t,\mathit{\theta })=\mathbf{A}(t-\mathrm{dt},\mathit{\theta })N(t-\mathrm{dt}),\end{array}$$

where $\widehat{w}$ is computed from $\widehat{N}$ following Eq. (2). The fit of the model is quantified using two numbers: the loss rate (∑(θ), lower indicates better fit) and the correlation between the observed and modeled mean biovolumes ${\stackrel{\mathrm{‾}}{v}}_{\text{obs}}$ and ${\stackrel{\mathrm{‾}}{v}}_{\text{mod}}$ over the diel cycle ($\text{corr}({\stackrel{\mathrm{‾}}{v}}_{\text{obs}},{\stackrel{\mathrm{‾}}{v}}_{\text{mod}})$, higher indicates better fit). Table 1 provides the model parameters being optimized.

Growth rate and loss rate

Once optimal parameters are identified, the model estimates a population intrinsic growth rate μ_size, and a specific loss rate l, integrated over a 24 h period. The method uses the fact that the observed size distribution N is the result of both growth and loss processes, while the time projection of the initial size distribution N(0) using the model $\widehat{N}$ is only the result of growth processes. The growth rate is calculated at each time step following Eq. (9) and integrated over 24 h. A total of 200 iterations by a Markov Chain Monte Carlo were run to estimate the standard deviation of group-specific growth rates.

$$\begin{array}{}\text{(9)}& {\mathit{\mu }}_{\text{size}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{dt}}}\text{ln}\left({\displaystyle \frac{{\sum }_{i=\mathrm{1}}^m{\widehat{N}}_i(t+\mathrm{dt})}{{\sum }_{i=\mathrm{1}}^m{\widehat{N}}_i\left(t\right)}}\right),\end{array}$$

where i is the ith size class, $\widehat{N}$ the predicted size distribution (cells cm⁻³), m the number of size classes, dt the time step (h) and μ_size the growth rates (d⁻¹).

An independent growth rate estimation was obtained as ${\mathit{\mu }}_{\text{ratio}}=\text{ln}({\stackrel{\mathrm{‾}}{v}}_{\text{max}}/{\stackrel{\mathrm{‾}}{v}}_{\text{min}})$, where ${\stackrel{\mathrm{‾}}{v}}_{\text{min}}$ and ${\stackrel{\mathrm{‾}}{v}}_{\text{max}}$ are the minimum and maximum of the mean observed biovolume ${\stackrel{\mathrm{‾}}{v}}_{\text{obs}}$ over the diel cycle (Marrec et al., 2018). μ_ratio represents a minimum estimate of the daily growth rate that would be observed if cells synchronously only grew from the time ${\stackrel{\mathrm{‾}}{v}}_{\text{min}}$ is observed (typically dawn) to the time ${\stackrel{\mathrm{‾}}{v}}_{\text{max}}$ is observed (typically dusk) and only divided while $\stackrel{\mathrm{‾}}{v}$ decreases. Since the model allows for any cell to grow, divide or be at equilibrium over the entire integration period (asynchronous populations), μ_size is expected to be higher than μ_ratio. In practice, μ_ratio is sensitive to noise in the data and is only provided here as an alternative estimate of the growth rate that does not rely on the model.

The population loss rate l is obtained by the difference between the intrinsic growth rate μ_size(t) and the temporal change in logarithmic observed size distribution N, which represents the net growth rate $r\left(t\right)={\mathit{\mu }}_{\text{size}}\left(t\right)-l\left(t\right)$ so that

$$\begin{array}{}\text{(10)}& \stackrel{\mathrm{‾}}{l}=\stackrel{t}{\int }{\mathit{\mu }}_{\text{size}}\left(t\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{dt}}}\text{ln}{\displaystyle \frac{N(t+\mathrm{dt})}{N\left(t\right)}}.\end{array}$$

Figure 3 Temporal evolution of sea surface conservative temperature (Θ) in degrees Celsius (^∘C), absolute salinity (S_A) in grams per kilogram (g kg⁻¹) and phytoplankton abundances in cells per milliliter (cells mL⁻¹) from 11 May 00:00 to 13 May 12:00 UTC. The vertical colors correspond to the two water masses separated by the front (see Fig. 1).

3 Results

3.1 Spatio-temporal distribution of phytoplankton abundances in the two water masses

The sampling strategy adopted during PROTEVSMED-SWOT enabled us to sample two water masses with different properties. The map of the satellite-derived surface [Chl a] shows a higher concentration in the northern part of the sampling route, corresponding with the older AW, than in the southern part, corresponding with the younger AW (Fig. 1b). Figure 3 shows the properties of the sea surface water as a function of time (from 11 May 00:00 to 13 May 12:00 UTC) along the sampling route. The older AW is characterized by a colder temperature and higher salinity than the younger AW. Figure 3 also displays the abundances of each phytoplankton group over these two water masses. Synechococcus and Pico2 are the most abundant. They present a clear surface distribution pattern, with high abundances in the warm and low-salinity water, corresponding with the young AW. A similar distribution is observed for Pico1, Pico3 and RNano but with lower abundances than Synechococcus and Pico2. The abundances of SNano, PicoHFLR and cryptophytes show less contrast along the cruise than the previous groups; nonetheless, the highest abundances can be distinguished in the younger AW, in particular during the second and third passage (transect) across this water mass. Finally, microphytoplankton is the less abundant group, but it clearly shows a contrast between the two water masses, opposite to that of the other phytoplankton groups.

Figure 4 The background color represents the Synechococcus cell-size distribution (i.e., biovolume in µm³) observed (a, b) and predicted by the model (c, d) in the older AW (a, c) and in the younger AW (b, d) over 24 h. The black dots represent the mean of the biovolume (${\stackrel{\mathrm{‾}}{v}}_{\text{obs}}$ and ${\stackrel{\mathrm{‾}}{v}}_{\text{mod}}$), and the red line represents the irradiance ($\mathrm{\mu }\mathrm{E}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$).

Figure 5 Same as Fig. 4 but for RNano.

Figure 6 Same as Fig. 4 but for SNano.

Table 2 Means of biovolumes observed (${\stackrel{\mathrm{‾}}{v}}_{\text{obs}}$) and modelized (${\stackrel{\mathrm{‾}}{v}}_{\text{mod}}$) in µm³, growth rates (μ_size, μ_ratio in d⁻¹), and loss rates (l, in d⁻¹) for the phytoplankton groups, in the older and younger AW, as well as the model fit parameters (see Sect. 2.3).

3.2 Phytoplankton cellular growth and division in the two water masses

The phytoplankton diurnal cycle was reconstructed in the two water masses using the size-structured population model originally developed by Sosik et al. (2003). Figures 4–6 represent the phytoplankton size distribution (i.e., biovolume) observed in situ and predicted by the model over 24 h for Synechococcus, RNano and SNano, respectively. From the predicted biovolume it is possible to derive specific growth (μ_size) and loss (l) rates, summarized in Table 2 for the different phytoplankton groups in the two water masses, along with the metrics of model performance. We also attempted to model the diurnal cycle for the picophytoplankton groups, i.e., Pico1, Pico2, Pico3 and PicoHFLR. However, their very noisy size distributions prevented us from obtaining reliable growth rate estimates. Similarly, microphytoplankton and cryptophytes were not abundant enough to allow a reliable determination of their abundances and cell cycles. These cytometric groups are thus not considered further in this study.

For Synechococcus, in the older AW the prediction of the model (i.e., predicted biovolume) is similar to the observed size distribution (i.e., observed biovolume). Both display a day-long large size-class distribution centered approximately on 0.3 µm³. In the younger AW (Fig. 4a and c) the distributions of observed and predicted biovolume are narrower than in the older AW and centered approximately on 0.2 µm³ (Fig. 4b and d). As a consequence, the older AW is populated by larger cells of Synechococcus (mean observed biovolume ${\stackrel{\mathrm{‾}}{v}}_{\text{obs}}$ = 0.38 ± 0.04 µm³) than in the younger AW (mean biovolume ${\stackrel{\mathrm{‾}}{v}}_{\text{obs}}$ = 0.21 ± 0.04 µm³) (Table 2). Growth and loss rates also differ between the two water masses. In the older AW, the large cells of Synechococcus have a growth rate of μ_size = 0.24 ± 0.91 d⁻¹ and a loss rate of l = 0.36 d⁻¹, whereas in younger AW the smaller cells are characterized by higher growth (μ_size = 0.68 ± 1.56 d⁻¹) and loss (l = 0.48 d⁻¹) rates.

Relative to Synechococcus, cell-size distribution and growth and loss rates are less contrasted between the older and younger AW for SNano (Fig. 6) and even more so for RNano (Fig. 5). The mean observed RNano biovolumes are similar in the older and younger AW (63.5 ± 2.67 µm³ and 61.2 ± 5.23 µm³, respectively) (Table 2). For SNano, similar to Synechococcus, the older AW is predominantly composed of larger cells (${\stackrel{\mathrm{‾}}{v}}_{\text{obs}}$ = 85.0 ± 1.98 µm³) than in the younger AW (${\stackrel{\mathrm{‾}}{v}}_{\text{obs}}$ = 63.8 ± 4.45 µm³). For both nano groups, growth rates are generally very low in both water masses (μ_size < 0.1 d⁻¹). Loss rates are higher than growth rates, except for RNano in the younger AW (negative loss rate implying an external input of cells such as by advection). However, the corresponding optimization factor is the highest observed across the six modelizations, indicating that this result is subject to caution.

4 Discussion

4.1 The phytoplankton diurnal cycle

Although it has been clearly demonstrated that phytoplankton plays a fundamental role in the ocean ecosystem functioning (Watson et al., 1991; Field et al., 1998; Allen et al., 2005), numerous questions remain about its population dynamics in relation to finescale structures.

Coupling high-resolution in situ flow cytometry measurements in two contrasted water masses with the size-structured population model developed by Sosik et al. (2003) allowed us to characterize the structure of phytoplankton and to reconstruct its diel cycle of cell growth and division on both sides of a finescale front. The growth and loss rates (μ_size and l) found for Synechococcus are of the same order of magnitude as those obtained by Marrec et al. (2018) in the northwestern Mediterranean Sea using the same method. In Sect. 3.2, we showed that the largest cells of Synechococcus were found in the older AW. These Synechococcus cells are characterized by a larger range of biovolume and lower growth and loss rates than those located in the younger AW (Table 2). They are on average larger than in the younger AW, as they grow slower at the population scale and divide less. Conversely, in the younger AW the distribution of the Synechococcus biovolume is narrower, which could be explained by cells being more active, more homogeneous in terms of size (biovolume) and better synchronized, leading to a smaller spread of the cell biovolume (Fig. 4b and d) with a dominance of small Synechococcus cells (Fig. 3). This also explains why higher abundances of Synechococcus are found in the younger AW (Fig. 3). Interestingly, the resulting net growth rate (growth minus loss) is negative in the older AW and positive in the younger AW.

Figure 7 The contrasted distribution of phytoplankton in the frontal area. The circles represent the abundances of the several phytoplankton groups in the two water masses separated by the front. The boxes indicate the biovolume observed (${\stackrel{\mathrm{‾}}{v}}_{\text{obs}}$) and the growth rates (μ_size) for each phytoplankton group, as estimated from the model. Figure adapted from Tzortzis et al. (2021).

Results are more difficult to interpret for the nanoplankton groups RNano and SNano, expected to be mostly dominated by diatoms in the Mediterranean Sea (Marty et al., 2002; Siokou-Frangou et al., 2010; Navarro et al., 2014; El Hourany et al., 2019), especially in frontal systems (Claustre et al., 1994). RNano and SNano diel cycles are not as well defined as for Synechococcus, leading to very small estimates of growth rates by the model. Optimization factors (linked to the mean squared difference between observed and predicted normalized size distributions) are relatively high and/or temporal correlations between observed and predicted mean biovolume relatively low, indicating that these results must be considered with caution. Nevertheless, our results suggest much lower growth and loss rates for nanoplankton than for Synechococcus and potentially higher growth rates in the younger AW, similar to Synechococcus (excluding the likely unrealistic loss rate obtained for RNano in the younger AW).

4.2 Influence of the frontal system on the phytoplankton dynamics

Our previous article (Tzortzis et al., 2021) provided a description of the hydrodynamics and the hydrology of the region. In the following, we attempt to establish the potential link between the characteristics of the two AWs separated by the front, the physical forcings associated with this frontal structure and the particular distribution of phytoplankton in terms of cell size and abundances. Figure 7 summarizes the physical forcing evidenced in this frontal area in the previous publication during the PROTEVSMED-SWOT cruise, superimposed with the biovolumes and the abundances of the different phytoplankton groups sampled in situ by the automated flow cytometer.

The older AW is characterized by larger cells of Synechococcus and nanophytoplankton with low abundances, low intrinsic growth rates and negative net growth rates, suggesting an older, declining population, whereas the younger AW is dominated by small cells with high abundances and at least for Synechococcus high intrinsic growth rates and a positive net growth rate, suggesting a slightly growing or stable population (nanoplankton results in the younger AW are subject to caution, as optimization factors are relatively high). Furthermore, microphytoplankton (i.e., the largest type of phytoplankton) is more abundant in the older AW than in the younger AW. The early experimental works of Marshall and Orr (1928), Jenkin (1937) and Huisman (1999) have well established that light and nutrients are essential for phytoplankton growth. The reconstruction of the circadian cycle indicates that irradiance was similar in the two water masses (Figs. 4–6, red lines), with corresponding daily total irradiances of 286 and 299 µE m⁻² for the older AW and the younger AW, respectively (Fig. A2). That is why the availability of light seems not to be the principal cause explaining the difference of phytoplankton dynamics and its distribution in the two AWs. An other possible explanation is that these two water masses are characterized by different nutrient concentrations, thus favoring certain phytoplankton groups. Bethoux (1989) and Schroeder et al. (2010) have observed that the older AW is slightly more enriched with nutrients than the younger AW because the older AW receives nutrient inputs from the continent (river discharges, rain, wind) during its circulation across the Mediterranean basin. Unfortunately, it was not possible for both technical and funding reasons to perform nutrient measurements during the 2018 cruise, so we cannot conclusively assess nutrient patterns during the cruise. Assuming that the nutrient distribution across the two water masses was similar to what was previously measured by Bethoux (1989) and Schroeder et al. (2010), we propose that higher nutrient concentrations in the older AW explain the observed phytoplankton cell-size and abundance distributions. Our hypothesis is supported by similar observations by Jacquet et al. (2010) and Mena et al. (2016), who also found the highest abundances of the small phytoplankton (Synechococcus and picophytoplankton) in the most oligotrophic waters, i.e., the younger AW. Furthermore, previous studies have shown that the proportion of picophytoplankton in the total phytoplankton biomass is higher in oligotrophic regions than in mesotrophic or eutrophic regions (Zhang et al., 2008; Cerino et al., 2012). Indeed, their better surface to size ratio due to their small size confers them with a better capacity to inhabit areas with very low nutrient concentrations compared to larger phytoplankton (Kiørboe, 1993; Mara nón, 2015). Since our study area is always oligotrophic (Moutin et al., 2012), a small variation of the nutrient concentration (typically ≤ 0.1 µM of nitrate) is sufficient to generate a higher abundance of picophytoplankton. Some studies have attempted to link hydrological condition and the phytoplankton dynamic (e.g., Qasim et al., 1972; Brunet et al., 2006; Mara nón et al., 2012). However, their results showed that the influence of these hydrological parameters on the phytoplankton growth and distribution was difficult to estimate, compared to the effects of nutrient availability and radiation exposure.

Other physical processes occurring at the front can explain the different dynamics of phytoplankton groups. The works of Lévy et al. (2001), Pidcock et al. (2016) and Mahadevan (2016) have highlighted that the availability of light and nutrients is driven by physical dynamics such as vertical velocities. The computation of the vertical motions in the frontal area, as represented in Fig. 7 (see also Fig. A3), shows the presence of upwellings and downwellings in the frontal area. However, due to the lack of nutrient measurements during the cruise, we are not able to quantify the impact of these vertical velocities.

4.3 Limitations of the study and recommendations

Growth and loss rates were estimated using the size-structured population model originally developed by Sosik et al. (2003), which was fitted to a measured diel cycle of cell-size distributions. More precisely, the rates were calculated based on the fitted size distribution predicted by the model and its comparison with the observed size distribution. Because of this, results are sensitive to noise in the measured size distributions. We could not obtain reliable results for the picophytoplankton groups due to noisy distributions because they probably contained several taxa with differing dynamics (Siokou-Frangou et al., 2010; Le Moal et al., 2011). To take into account this constraint, in future experiments, sorting by flow cytometry and identification with a microscope and/or genetics analysis should be planned to identify taxa in the various phytoplanktonic groups defined by flow cytometry. Nevertheless, these techniques are not easily applicable to large numbers of samples contrary to automated flow cytometry, which means that a careful selection of samples will be necessary.

For the taxonomic groups where reasonable size distributions could be estimated over a diel cycle, the model fit was evaluated using two metrics: the optimization loss rate ∑(θ) and the correlation between the observed and modeled mean biovolumes over the diel cycle $\text{corr}({\stackrel{\mathrm{‾}}{v}}_{\text{obs}},{\stackrel{\mathrm{‾}}{v}}_{\text{mod}})$ (Table 2). These metrics, as well as visual comparisons of the modeled and observed size distributions (Figs. 4–6), indicate differing degrees of confidence in our estimated growth and loss rates, with the highest confidence obtained for Synechococcus. In future experiments, these rate estimates could be improved by more accurately measuring the phytoplankton diel cycle (i.e., by continuously sampling the same water mass over 24 h rather than by compiling several days to reconstruct a diel cycle). Furthermore, coupling flow cytometry with NanoSIMS analysis, as in the works of Bonnet et al. (2016) and Berthelot et al. (2019), could also be useful to get independent estimates of the growth rates, although the cost and the successive incubations required by this methodology are not adapted to high-frequency sampling of finescale ocean structures.

Overall, while estimating growth and loss rates by fitting a model to automated flow cytometry data remains limited by our capacity to accurately resolve the size distribution of independent phytoplankton groups over a full diel cycle, the method used here has several advantages. Other methods involve incubations that are dependent on the accurate reproduction of the marine environment in incubators. By contrast, automated flow cytometry as applied here measures the temporal evolution of phytoplankton in situ. The automated CytoSense flow cytometer, deployed underway, requires little maintenance or manipulation during the cruise contrary to time-consuming sampling and incubations. As such, while growth and loss rates obtained from automated flow cytometry would benefit from independent validation, they have the potential to provide in situ estimates of biological rates that are traditionally difficult to measure.

5 Conclusion and perspectives

Phytoplankton structure and dynamics are a complex result of many interacting biological and physical phenomena. Finescale structures, and in particular fronts, generate vertical velocities which displace phytoplankton cells and nutrients in the water column, thus influencing phytoplankton communities. These mechanisms are only partially understood because the spatial scale of these structures and their ephemeral nature makes them particularly difficult to monitor in situ. In particular, the specific growth rates for the various phytoplankton groups, while being a key quantity for explaining the community structure, are challenging to obtain. In this study, we addressed this problem by monitoring the dynamics of several phytoplankton groups in two distinct water masses both in terms of hydrology and phytoplankton abundances. In the Mediterranean Sea, the low nutrient content is indeed the perfect condition when addressing this question because even weak horizontal or vertical nutrient redistributions associated with the finescale circulation are likely to result in a biological response (Talmy et al., 2014; Hashihama et al., 2021).

The originality of our work resides in the combination of (i) a Lagrangian sampling strategy to adaptively track a frontal region for several days, (ii) a high spatiotemporal resolution of the phytoplankton community thanks to flow cytometry and (iii) the use of a size-structured population model for reconstructing the diurnal cycle of several phytoplankton groups and for identifying contrasted dynamics in the two water masses.

For Synechococcus and nanophytoplankton, we found a higher cell size in the older AW located north of the front, associated with lower abundances. A possible explanation is that the older AW is more enriched in nutrients than the younger AW, thus favoring larger cells. This remains a hypothesis because of a lack of nutrient data. Besides the employment of a Lagrangian adaptive strategy, another novelty of our study is that we applied the Sosik et al. (2003) model on several phytoplankton groups identified by flow cytometry, whereas previous studies only applied it to Synechococcus and Prochlorococcus (Ribalet et al., 2010; Hunter-Cevera et al., 2014; Marrec et al., 2018; Fowler et al., 2020) or to certain types of diatoms (Dugenne et al., 2014). We obtained good results for Synechococcus and nanophytoplankton. Our analysis has been less conclusive for the other phytoplankton groups because of the presence of confounding effects and in particular in the limitation of flow cytometry in taxonomic resolution. A further limitation of our study has been the lack of concomitant nutrient data, which hindered the possibility of testing a mechanistic hypothesis, leading to the emergence of type dominance.

Built over the experience that we acquired with the study presented here, the recent BioSWOT-Med cruise has been lead in the southwestern Mediterranean in spring 2023 to address these limitations. This cruise added to the methodological approach presented here with high-resolution, high-precision nutrient measurements (necessary considering the oligotrophy of the Mediterranean Sea), as well as a better taxonomic resolution based on DNA metabarcoding and microscopy. The BioSWOT-Med cruise also included zooplankton and virus sampling aimed at estimating zooplankton grazing and viral lysis on the different phytoplankton groups to better characterize the cell loss. In the long term, we believe that these types of studies pave the way to the direct integration of growth rates in biogeochemical models (Cullen et al., 1993; Baklouti et al., 2006) and should eventually provide a better assessment of the biogeochemical contribution of phytoplankton in oligotrophic ecosystems as well as better predict its evolution in the context of global change.

Appendix A: Identification of the phytoplankton functional groups by flow cytometry

Up to nine groups of phytoplankton have been identified on the cytograms (Fig. A1) thanks to their light scatter (forward scatter, FWS, and sideward scatter, SWS) and fluorescence intensities (red fluorescence, FLR, and orange fluorescence, FLO). These groups have been called the conventional names used by flow cytometrists, i.e., some groups relate to taxonomy (Synechococcus, cryptophytes), while others relate to a range of sizes (picoeukaryotes, nanoeukaryotes) as described by Sieburth et al. (1978). Synechococcus (Syn in Fig. A1c) is a prokaryotic picophytoplankton that can be distinguished from the other picophytoplankton owing to its high FLO intensity, induced by phycoerythrin pigment content. Cryptophytes (Crypto in Fig. A1c) were also discriminated from the other groups, as they also produce a characteristic orange fluorescence induced by phycoerythrin. Concerning the other phytoplankton groups, four eukaryotic picophytoplankton groups were put in evidence: Pico1 (on Fig. A1c), characterized by lower FLR and FLO intensities than Synechococcus; Pico2 and Pico3 (on Fig. A1d) with higher FWS, SWS and FLR intensities than Pico1; and PicoHFLR (on Fig. A1a), which has a high FLR signal induced by Chl a. We defined two distinct nanophytoplankton groups (SNano and RNano) according to their high FLR and FLO intensities. SNano exhibits a higher $SWS/FWS$ ratio and SWS intensity than RNano (Fig. A1b and a). Finally, microphytoplankton (Micro) is characterized by the highest FLR and FWS intensities (Fig. A1c).

Figure A1 Cytograms obtained with the CytoSense automated flow cytometer. Synechococcus are in dark blue (Syn), the picophytoplankton with the lowest FLO in orange (Pico1), the picophytoplankton with intermediate FWS in light blue (Pico2), the picophytoplankton with the highest FWS in purple (Pico3), the picophytoplankton with a high red fluorescence in pink (PicoHFLR), the nanophytoplankton with high $SWS/FWS$ ratio in yellow (SNano) and higher SWS intensities than the other nanophytoplankton (RNano) in green, the cryptophytes in grey (Crypto), and the microphytoplankton in red (Micro). The flow cytometry units for both fluorescence and light scatter are arbitrary (a.u.). Figure extracted from Tzortzis et al. (2021).

Figure A2 Reconstruction of irradiance during 24 h in the older AW (a) and the younger AW (b). Computation of trapezoidal integration of irradiance in the older AW, E1 = 286 µE m⁻², and in the younger AW, E2 = 299 µE m⁻².

Figure A3 Vertical velocities at 25 m (a) and 85 m (b), calculated with the omega equation. Figure extracted from Tzortzis et al. (2021).