Théo Garcia
email: theo.garcia@mio.osupytheas.fr

Daniela Banaru

Loïc Guilloux

Véronique Cornet

Gérald Gregori

François Carlotti

Daniela Bănaru

¶ Loïc

Véronique CORNET 1& Guilloux 1&

Gérald Gregori 1&

Temporal changes in zooplankton indicators highlight a bottom-up process in the Bay of Marseille (NW

1-Introduction

Zooplankton are short-lived organisms, highly sensitive to environmental changes, which makes their long-term monitoring an excellent approach to provide responses to management issues [START_REF] Beaugrand | Monitoring pelagic ecosystems using plankton indicators[END_REF][START_REF] Pitois | A first approach to build and test the Copepod Mean Size and Total Abundance (CMSTA) ecological indicator using in-situ size measurements from the Plankton Imager (PI)[END_REF]. Previous time series analyses have shown that zooplankton shifts can be observed at various scales: at large scale through climate-related processes [START_REF] Planque | Long-term changes in zooplankton and the climate of the North Atlantic[END_REF][START_REF] García-Comas | Zooplankton long-term changes in the NW Mediterranean Sea: Decadal periodicity forced by winter hydrographic conditions related to large-scale atmospheric changes[END_REF], at regional scale in relation with regional meteorological phenomena, and at local scale (e.g. bays) most often in response to human pressure [START_REF] Berline | Intercomparison of six Mediterranean zooplankton time series[END_REF][START_REF] Fanjul | Zooplankton variability at four monitoring sites of the Northeast Atlantic Shelves differing in latitude and trophic status[END_REF]. These multiple stressors alter the physiology, behavior and phenology of organisms, and consequently population distributions [START_REF] Hays | Climate change and marine plankton[END_REF][START_REF] Mackas | Changing zooplankton seasonality in a changing ocean: Comparing time series of zooplankton phenology[END_REF]. Nonetheless the multiple scales of forcing make it difficult to disentangle causation.

The increase in zooplankton observations has made it possible to validate several indicators such as those based on individual taxa, species assemblages, diversity indices or size spectra, which are useful for characterizing the state of pelagic ecosystems [START_REF] Beaugrand | Monitoring pelagic ecosystems using plankton indicators[END_REF]. The development of multi-indicator approaches has made it possible to better compare them, and thus to identify their complementary contributions. Therefore nowadays quantification of the roles of zooplankton in biogeochemical fluxes and trophic fluxes towards higher levels requires, in addition to information on biomass and size spectra, knowledge of taxonomic distributions and associated ecological traits.

In the Gulf of Lion (NW Mediterranean), earlier studies dedicated to zooplankton were focused on quantification of either biogeochemical flux or pelagic trophic flows, and were all carried out during oceanographic cruises over short periods (less than a month) [START_REF] Pagano | Vertical Migrations and Feeding Rhythms of Mesozooplanktonic Organisms in the Rhône River Plume Area (North-west Mediterranean Sea)[END_REF][START_REF] Razouls | Spatial distribution and seasonal variation of mesozooplankton biomass in the Gulf of Lions (northwestern Mediterranean)[END_REF][START_REF] Champalbert | Characteristics of zooplankton standing stock and communities in the Western Mediterranean Sea: Relations to hydrology[END_REF][START_REF] Gaudy | Space and time variations in zooplankton distribution south of Marseilles[END_REF][START_REF] Plounevez | Diet, feeding behaviour and trophic activity of the anchovy (Engraulis encrasicolus L.) in the Gulf of Lions (Mediterranean Sea)[END_REF][START_REF] Gaudy | Biomass, metabolism and nutrition of zooplankton in the Gulf of Lions (NW Mediterranean)[END_REF][START_REF] Raimbault | Research activities in the Gulf of Lion (NW Mediterranean) within the 1997-2001 PNEC project[END_REF][START_REF] Palomera | Small pelagic fish in the NW Mediterranean Sea: An ecological review[END_REF][START_REF] Espinasse | Defining zooplankton habitats in the Gulf of Lion (NW Mediterranean Sea) using size structure and environmental conditions[END_REF][START_REF] Saiz | The Zooplankton[END_REF][START_REF] Bourg | Trophic niche overlap of sprat and commercial small pelagic teleosts in the Gulf of Lions (NW Mediterranean Sea)[END_REF]. Recent increased interest in the role of zooplankton as trophic resources for planktivorous fish has been prompted by the reduction in size of two planktivorous fishes (European pilchard, Sardina pilchardus, and European anchovy, Engraulis encrasicolus) in the Gulf of Lion in 2008 (NW Mediterranean) [START_REF] Van Beveren | Rapid changes in growth, condition, size and age of small pelagic fish in the Mediterranean[END_REF][START_REF] Bănaru | Trophic structure in the Gulf of Lions marine ecosystem (north-western Mediterranean Sea) and fishing impacts[END_REF][START_REF] Van Beveren | Predator-prey interactions in the face of management regulations: changes in Mediterranean small pelagic species are not due to increased tuna predation[END_REF][START_REF] Gfcm | Scientific Advisory Committee on Fisheries (SAC) Working Group on Stock Assessment of Small Pelagic species[END_REF]. Various studies have led to rejection of the hypotheses of the impact of overfishing, increased predation, disease and a decline in recruitment [START_REF] Saraux | Small pelagic fish dynamics: A review of mechanisms in the Gulf of Lions[END_REF]. The observation of a drastic reduction in the size of individuals [START_REF] Van Beveren | Rapid changes in growth, condition, size and age of small pelagic fish in the Mediterranean[END_REF][START_REF] Saraux | Small pelagic fish dynamics: A review of mechanisms in the Gulf of Lions[END_REF][START_REF] Brosset | Linking small pelagic dietary shifts with ecosystem changes in the Gulf of Lions[END_REF] has been linked to a deficit in individual growth attributed to potential changes in plankton resources, although this has not yet been proven. To shed light on this potential planktonic change, the data from oceanographic cruises over the last few decades in the Gulf of Lion are too few and limited in time, and the closest locations of zooplankton monitoring time series in the Ligurian Sea [START_REF] García-Comas | Zooplankton long-term changes in the NW Mediterranean Sea: Decadal periodicity forced by winter hydrographic conditions related to large-scale atmospheric changes[END_REF][START_REF] Serranito | Decadal shifts of coastal microphytoplankton communities in a semi-enclosed bay of NW Mediterranean Sea subjected to multiple stresses[END_REF][START_REF] Fullgrabe | Zooplankton dynamics in a changing environment: A 13-year survey in the northwestern Mediterranean Sea[END_REF] or the Balearic Sea [START_REF] De Puelles | Zooplankton time-series in the Balearic Sea (Western Mediterranean): Variability during the decade 1994-2003[END_REF] are far apart.

A mesozooplankton monitoring survey that has been carried out since 2005 in the Bay of Marseille (BoM) provides a zooplankton time series in the Gulf of Lion which encompasses the critical shift period. Zooplankton sampling was added to an existing monitoring programme started in 1994 including hydrological and nutrient measurements, as well as sampling for a few biotic variables (pigment, bacteria and phytoplankton). The sampling station at the center of the BoM is subject to multiple regional environmental influences such as winter deep convection phenomena from the Provencal basin off the Gulf of Lion, north-westerly winds (known as Mistral), which induce coastal upwelling events [START_REF] Millot | Wind induced upwellings in the gulf of lions[END_REF], intrusions from the Rhone River under certain wind conditions [START_REF] Pinazo | MASSILIA: Modélisation de la baie de MArSeILLe : Influence des apports Anthropiques de la métropole sur l'écosystème marin[END_REF], intrusions from the northerly current from the Lingurian Sea [START_REF] Millot | Circulation in the Western Mediterranean Sea[END_REF], and local anthropogenic influences due to the proximity of the Marseille agglomeration (e.g. sewage plants) [START_REF] Millet | Unexpected spatial impact of treatment plant discharges induced by episodic hydrodynamic events: Modelling Lagrangian transport of fine particles by Northern Current intrusions in the bays of Marseille (France)[END_REF].

In this paper, we aimed to test the hypothesis of a bottom-up relationship between environmental conditions, phytoplankton, and zooplankton in the BoM. To achieve this goal we will first describe a number of zooplankton indices and examine their interannual and seasonal variations in order to identify major changes over the last two decades. We will then relate zooplankton changes to variations in other monitored abiotic and biotic parameters, and discuss how the observed changes might be the result of bottom-up forcing. Finally, we will discuss how this could affect planktivorous fish.

2-Material and methods

2.1-Environmental data

Sampling was carried out in the Bay of Marseille (N-W Mediterranean Sea) at the Frioul monitoring station (60m depth, 43.24°N; 5.29°E) of the long-term national littoral observation program SOMLIT (www.somlit.fr), in the eastern part of the Gulf of Lion (Fig 1). This monitoring survey has generated an extensive collection of physical, chemical, and biological data (see summary of the variables in Table 1). Twice a month since 1994, a vertical CTDoxygen-fluorometer cast of the whole water column (0-55m) and Niskin bottle sampling were performed at three depths (subsurface -1 m-, bottom -55 m-and fluorescence maximumvariable). The values of each measurement obtained at the three depths were averaged. Climatic and meteorological indices completed the environmental database (Table 1).

2.2.1-Zooplankton size-fractioned dry weights

The aliquot of cod-end content maintained in marine water was fractioned on a sieve column in 6 size fractions: 80-200 µm, 200-300 µm, 300-500 µm, 500-1000 µm, 1000-2000 µm and >2000 µm. The material from each sieve was recovered on pre-combusted and pre-weighed Whatman GF/F filters and dried for 48h at 60°C for biomass measurement. Because the WP2 200 µm mesh size plankton net does not efficiently sample the smallest and largest size fractions, the 80-200 µm and >2000 µm fractions were not considered in our study.

2.2.2-Zooplankton digitalization

The other aliquot conserved in the formaldehyde solution was used for plankton digitalization following the procedure in [START_REF] Gorsky | Digital zooplankton image analysis using the ZooScan integrated system[END_REF] with a Zooscan. During the image processing, a proxy of individual size has been obtained by means of Equivalent Spherical Diameter (ESD, see [START_REF] Jennings | Particle size measurement: the equivalent spherical diameter[END_REF]).

A Convolutional Neural Network (CNN) has been trained for vignette classification [39]. The CNN predicted zooplankton vignettes in the following 13 classes: appendicularians, bivalves, calanoid copepods, ergasilida copepods, harpacticoid copepods, oithonoid copepods, nauplii copepods, crustaceans (which includes holoplankton, e.g. amphipod, cladocera, and meroplankton, e.g. decapod), chaetognaths, cnidarians, fish eggs, pteropods and salps. See the S2 file for more details concerning CNN classification, training, and performance.

2.3-Data analysis

Data analyses were conducted on R software version 4.1.2 [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF].

2.3.1-Description of zooplankton community structure indices

Size spectra were constructed using images with ESDs between 200 µm (limit of particle detection by Zooscan) and 1660 µm (i.e. outliers were excluded according to the Tukey method [START_REF] Tukey | Exploratory data analysis[END_REF]). Following the method described by Nerini and Ghattas [START_REF] Nerini | Classifying densities using functional regression trees: Applications in oceanology[END_REF], the shape of the zooplankton size spectra was studied by means of Functional Principal Component Analysis (FPCA). The FPCA is a statistical tool that ordinates sample units according to the characteristics (i.e. the shape) of their functional entity (i.e. here a curve/size spectrum). This methodology was applied as follows: (i) zooplankton size spectrum densities for all samples were estimated using Kernel density estimator; (ii) the density functions were converted into a continuous functional object, a curve, using 250 B-splines functions; (iii) the curves obtained from all samples were ordinated by means of a FPCA [START_REF] Ramsay | Introduction to Functional Data Analysis[END_REF]. We applied the FPCA on the coefficients of the B-spline expansion model as explained by Pauthenet et al. in [START_REF] Pauthenet | A Linear Decomposition of the Southern Ocean Thermohaline Structure[END_REF]. These new variables which concentrate the variance of the system are the principal components (PCs) and represent the most significant modes of the data variation.

The sample scores on the first two PCs were used as time series to describe the variations of the size spectrum shapes and by extension the structure of the zooplankton community size.

In addition, multidimensional scaling (MDS) was applied on the species abundance table to disentangle the diversity structure with Manhattan distance. The sample scores of the three MDS axes were used as time series to describe the variations of the diversity structure.

The R packages 'fda' [START_REF] Ramsay | Introduction to Functional Data Analysis[END_REF] and 'vegan' [START_REF] Dixon | a package of R functions for community ecology[END_REF] were respectively used to perform the FPCA and MDS analysis.

2.3.2-Exploring interannual changes in zooplankton data

Time series data were averaged on a monthly basis and missing values were imputed using the Multiple Imputation with Principal Component Analysis (MIPCA) procedure using the R package 'MissMDA' [START_REF] Audigier | Multiple imputation for continuous variables using a Bayesian principal component analysis[END_REF]. Trends in monthly time series were highlighted using Local Polynomial Regressions (LOESS, using 75% of the neighborhood data). Breakpoint detection procedure ('strucchange' package in R [START_REF] Zeileis | strucchange: An R Package for Testing for Structural Change in Linear Regression Models[END_REF]) was applied to highlight the structural change on every single time-series. Until three breakpoints were investigated, the optimum number of breakpoints was defined according to the Bayesian Information Criterion (BIC).

2.3.3-Zooplankton seasonality

Level plots were used to describe seasonal patterns in the various zooplankton series: biomass, abundance (see 2.2-Zooplankton data) and structure descriptors (see 3.1-Community structure indices). In order to obtain a smoothed picture of the level plot, LOESS functions were applied over 30 equally spaced points on the year/month grids. Seasonal patterns were assessed by fitting sinusoidal models. The seasonal pattern was considered significant when the sinusoidal effect was significant (p-value<0.05). The seasonal onset timing of each zooplankton group was defined as the date when 20% of the yearly cumulative value was reached [START_REF] Mackas | Changing zooplankton seasonality in a changing ocean: Comparing time series of zooplankton phenology[END_REF][START_REF] Greve | Predicting the seasonality of North Sea zooplankton[END_REF].

2.3.4-Investigating relationships between environmental and zooplankton time series

Dynamic Factor Analysis (DFA, [49]) is a statistical method that aims to find common trends among observational data. In our case DFA was applied in order to find common trends between environmental and zooplankton data. DFA belongs to the family of the state-space models and allows determination of M hidden common trend(s) among N standardized (to reduce the impact of the variable scale on the analysis) data time series, when M<<N. DFA models relate observational time series (y) to the sum of the M hidden trends linear combination and noise.

This can be represented by the following equation (1):

𝑦 𝑡 = 𝛤𝛼 𝑡 + 𝜀 𝑡 ; 𝑤ℎ𝑒𝑟𝑒 𝜀 𝑡 ~𝑁(0, 𝑅) (1)
Γ is a matrix of dimensions N x M and contains the factor loadings of the N factors for M detected hidden trends, αt is a vector of length M containing the M common trends at time t. εt is the noise component which follows a normal distribution of mean 0 and variance-covariance matrix R.

DFA models require user specification for certain parameters:

-Firstly, different assumptions on the noise, R, variance-covariance matrix among factors can be made:

-same variance and no covariance ('diagonal and equal')

-different variances and no covariances ('diagonal and unequal')

-same variances and same covariances ('equalvarcov')

-different variances and covariances ('unconstrained')

-Secondly, the number of M hidden trends must be specified as low as possible to facilitate interpretation (M<<N).

Correction of Akaike Information Criterion (AICc) was applied to select the most likely data model.

Multiple DFA models were computed on:

(1) -all environmental data for the 2005-2020 period. From one to three common trends were investigated among the 34 environmental times series with the four different assumptions on the R matrix. A total of 12 models were computed on environmental time series.

(2) -zooplankton data for the 2005-2020 period. From one to three common trends were investigated among 11 zooplankton times series with the four different assumptions on the R matrix. A total of 12 models were computed on biological data.

(3) -environmental winter conditions (average values between January and March) and dates of zooplankton seasonal onset. From one to three common trends were investigated among the 34 environmental and 19 zooplankton time series with the four different assumptions on the R matrix. A total of 12 models were computed.

As DFA aimed at finding common trends among short time series, in (1) and ( 2), time series were averaged at a biannual time step. This allowed to reduce the computing time for the models. The series were de-seasonalized (i.e. when seasonal patterns were detected, by means of partial autocorrelation functions) to focus on interannual variations. DFA can be performed with missing data. For model 1 and 3 some time series were shorter (i.e. micro-, pico-, nanoplankton data started in 2009), the results of the models fitting computed for those variables were therefore not discussed before 2009.

Dynamic factor analysis models were computed with the MARSS package [START_REF] Holmes | MARSS: Multivariate Autoregressive State-space Models for Analyzing Timeseries Data[END_REF].

3-Results

3.1-Community structure indices

The two first PCs of the FPCA performed on the size spectra curves explained 66. bivalves; cala: calanoids copepods; chae: chaetognaths; cnid: cnidarians; crust: crustaceans (without copepods); eggs: fish eggs; erga: ergasilida copepods; harp: harpacticoids copepods; naup: copepods nauplii; oith: oithonids copepods; pter: pteropods; salp: salp.

3.2-Zooplankton interannual variations and relation with the environmental trend

The

3.3-Zooplankton seasonality and relationship between springsummer production and winter environment

Most of the zooplankton series presented a clear seasonal pattern (Fig 7) but with different shapes according to groups and size fractions. Seasonality assessments were significant (p-values<0.05) for every variable except for the second axis of FPCA (see S4). The total zooplankton abundance and biomass, the four-biomass size-fractions between 200 and 2000 µm, and the sample scores on the first axis of the MDS followed a similar seasonal pattern to calanoids, oithonoids, ergasilida, bivalves, copepod nauplii and harpacticoids with their highest abundance values in spring (April-May). The remaining zooplankton groups (appendicularians, chaetognaths, cnidarians, crustaceans, pteropods, salps, fish eggs) and the two remaining MDS axes seemed to reach their highest values during summer.

4-Discussion

4.1-Seasonal and interannual variations of multiple zooplankton indicators

The seasonal patterns of the main zooplankton groups in the Bay of Marseille during the time series are consistent with zooplankton successions observed in Mediterranean coastal areas [START_REF] Berline | Intercomparison of six Mediterranean zooplankton time series[END_REF][START_REF] Fullgrabe | Zooplankton dynamics in a changing environment: A 13-year survey in the northwestern Mediterranean Sea[END_REF][START_REF] De Puelles | Seasonal spatial pattern and community structure of zooplankton in waters off the Baleares archipelago (Central Western Mediterranean)[END_REF][START_REF] Romagnan | Comprehensive Model of Annual Plankton Succession Based on the Whole-Plankton Time Series Approach[END_REF][START_REF] Del | Spatial and Temporal Long-Term Patterns of Phyto and Zooplankton in the W-Mediterranean: RADMED Project[END_REF][START_REF] Feuilloley | Temporal fluctuations in zooplankton size, abundance, and taxonomic composition since 1995 in the North Western Mediterranean Sea[END_REF], with a zooplankton community dominated by copepods in response to the spring phytoplankton blooms and a second late summer -fall peak of abundance. In the Gulf of Lion, phytoplankton-zooplankton peaks are associated for both periods with strong grazing activity [START_REF] Gaudy | Biomass, metabolism and nutrition of zooplankton in the Gulf of Lions (NW Mediterranean)[END_REF]. The second peak is associated with higher diversity indices (as shown in Fig. 7 and8) linked to a relative diminution of the copepod contribution (mainly calanoids and oithonoids) and the appearance of other taxonomic groups (either herbivorous taxa feeding on small size phytoplankton, such as salps and pteropods, or carnivorous taxa, such as chaetognaths, crustaceans and cnidarians), a general pattern already described in the NW Mediterranean Sea [START_REF] Champalbert | Characteristics of zooplankton standing stock and communities in the Western Mediterranean Sea: Relations to hydrology[END_REF].

Our results showed different interannual patterns among the size-fraction biomass dynamics (see Fig. 4). The sudden decline in 2008 for the largest size-fraction (1000-2000 µm) was not simultaneously observed for the three other fractions (between 200 and 1000 µm) which decreased between 2 and 4 years later. No major change in the community diversity was observed in 2008, while the diminution of the larger size fraction biomass was accompanied by a shift in the size structure of the zooplankton community. Previous studies in the Bay of Marseille have shown that the taxonomic assemblages are not the same in the four size fractions [START_REF] Bănaru | Seasonal variation of stable isotope ratios of size-fractionated zooplankton in the Bay of Marseille (NW Mediterranean Sea)[END_REF][START_REF] Chen | Seasonal variation in biochemical and energy content of size-fractionated zooplankton in the Bay of Marseille (North-Western Mediterranean Sea)[END_REF], impacting differently the biomass dynamics of each size fraction. The size-fractioned approach enabled us to highlight these differences between the size fraction biomass dynamics that would have gone unnoticed if considering the whole biomass. This revealed the strong link between the biomass size-fraction dynamic and community size structure.

Our results are also consistent with the zooplankton dynamics of the Bay of Calvi [START_REF] Fullgrabe | Zooplankton dynamics in a changing environment: A 13-year survey in the northwestern Mediterranean Sea[END_REF] and the Bay of Villefranche-sur-Mer [START_REF] Feuilloley | Temporal fluctuations in zooplankton size, abundance, and taxonomic composition since 1995 in the North Western Mediterranean Sea[END_REF] during the same decade. The main changes in the zooplankton community in the Bay of Marseille observed between 2012 and 2014 (see Fig. 4 & 6) appear to have occurred earlier than similar variations in the whole mesozooplankton community observed after 2015 at Villefranche-sur-Mer [START_REF] Feuilloley | Temporal fluctuations in zooplankton size, abundance, and taxonomic composition since 1995 in the North Western Mediterranean Sea[END_REF], where crustaceans (including copepods in their study) showed a diminution of abundance and a lower size spectrum slope (mainly due to the loss of small copepods). In their analysis of multiple time series in the Mediterranean Sea, Berline et al. [START_REF] Berline | Intercomparison of six Mediterranean zooplankton time series[END_REF] did not find any common zooplankton interannual patterns at basin scale.

Nevertheless, our results and those at Villefranche-sur-mer and Calvi suggest that the coastal regions of the Ligurian Sea might be similarly impacted by the environmental changes and by their connections through the Ligurian current [START_REF] Berline | A Connectivity-Based Eco-Regionalization Method of the Mediterranean Sea[END_REF]. As the re-analysis of time series with recent data proved that some environmental-biological processes relationships may break down [START_REF] Myers | When Do Environment-recruitment Correlations Work?[END_REF].

Reanalyses of zooplankton time series in the Mediterranean Sea with longer time series are required to assess some large-scale processes.

The different community level indicators (i.e. size-fraction biomasses, total zooplankton abundance, size-and diversity-structure) used in our study highlighted shifts occurring in different years (Fig. 4 &6) denoting structural changes that would not have been perceived by more aggregated variables (e.g. total abundance or total biomass). In 2014 the assemblage composition from image analysis was sufficient to explain the shift. However, a finer taxonomic reanalysis would certainly be necessary to understand the 2008 shift potentially due to changes within the calanoid group, a group with high diversity in species composition, size distribution and trophic behavior [START_REF] Kleppel | On the diets of calanoid copepods[END_REF][START_REF] Beaugrand | Biodiversity of North Atlantic and North Sea calanoid copepods[END_REF]. In addition, such a more detailed taxonomic identification, up to separation of adult and copepodite stages, could help to better interpret the fine changes in size spectrum, highly sensitive to the underlying assemblage compositions and their associated size distribution.

Overall the multiple indicator approach used in our study helped us to better characterize and understand the zooplankton community dynamic. Pitois et al. [START_REF] Pitois | A first approach to build and test the Copepod Mean Size and Total Abundance (CMSTA) ecological indicator using in-situ size measurements from the Plankton Imager (PI)[END_REF] also advocated a multiindicator approach as necessary for describing zooplankton structural changes.

4.2-Variations of the environmental conditions and link with zooplankton

In the present study, we aimed to understand the observed zooplankton variations within an ecosystemic framework at the study site (as shown in -As autotrophic planktonic groups show preferences in nutrient uptake [START_REF] Probyn | Nitrogen uptake by size-fractionated phytoplankton populations in Antarctic surface waters1[END_REF], the variations in the stoichiometry of nutrient sources and their concentrations may have . This supports the hypothesis of a regime shift in the whole pelagic food web through bottom-up cascading effects [START_REF] Sommer | Pelagic food web configurations at different levels of nutrient richness and their implications for the ratio fish production:primary production[END_REF][START_REF] Stibor | Copepods act as a switch between alternative trophic cascades in marine pelagic food webs[END_REF],

similarly to those described in Rietkerk et al. [START_REF] Rietkerk | Self-Organized Patchiness and Catastrophic Shifts in Ecosystems[END_REF]. The concomitant decline of calanoid copepods and increase in appendicularians and salps represented an indicator of such change in the zooplankton community. Because appendicularians and salps can feed on smaller particles, the trophic competition with calanoid copepods (feeding on larger particles) seemed to be reduced [START_REF] Flood | Filtration of colloidal melanin from sea water by planktonic tunicates[END_REF][START_REF] Katechakis | Feeding selectivities and food niche separation of Acartia clausi, Penilia avirostris (Crustacea) and Doliolum denticulatum (Thaliacea) in Blanes Bay (Catalan Sea, NW Mediterranean)[END_REF]. When copepods are unable to feed on large microphytoplankton, they often switch to microzooplankton [START_REF] Stoecker | Predation on Protozoa: its importance to zooplankton[END_REF][START_REF] Irigoien | Feeding selectivity and egg production of Calanus helgolandicus in the English Channel[END_REF] resulting in a higher trophic level of copepods [START_REF] Sommer | Pelagic food web configurations at different levels of nutrient richness and their implications for the ratio fish production:primary production[END_REF]. The hypothesis of such a relationship between the phytoplankton size and zooplankton assemblages has been suggested in the NW Mediterranean Sea [START_REF] Feuilloley | Temporal fluctuations in zooplankton size, abundance, and taxonomic composition since 1995 in the North Western Mediterranean Sea[END_REF][START_REF] Calbet | Annual Zooplankton Succession in Coastal NW Mediterranean Waters: The Importance of the Smaller Size Fractions[END_REF]. These changes in the Bay of Marseille were accompanied with the appearance of predators such as chaetognaths or cnidarians. In the Catalan Sea, the predation impact of chaetognaths on the standing stock of copepod was low, but chaetognaths can exert a high level of predation on specific copepod species and stages under food-limited conditions [START_REF] Duró | Distribution and trophic ecology of chaetognaths in the western Mediterranean in relation to an inshore-offshore gradient[END_REF]. A species level taxonomy would help to quantify the top-down pressure of predators on copepods species dynamics.

In parallel, by analyzing the winter environmental conditions, we observed a relationship between NAO, precipitation, temperature, and zooplankton phenology (Fig 9). Although at Mediterranean scale no relationship was found between NAO and zooplankton by Berline et al. [START_REF] Berline | Intercomparison of six Mediterranean zooplankton time series[END_REF], several studies suggested that winter conditions at NW Mediterranean Sea regional scale (inducing strong variation in the NAO index) impacted the following zooplankton spring peak [START_REF] García-Comas | Zooplankton long-term changes in the NW Mediterranean Sea: Decadal periodicity forced by winter hydrographic conditions related to large-scale atmospheric changes[END_REF][START_REF] Fullgrabe | Zooplankton dynamics in a changing environment: A 13-year survey in the northwestern Mediterranean Sea[END_REF][START_REF] De Puelles | Zooplankton time-series in the Balearic Sea (Western Mediterranean): Variability during the decade 1994-2003[END_REF], as shown by our results (and discussed in the previous section). Since 2014, winter deep convection (associated with climatic forcing NAO) declined in the NW Mediterranean Sea causing a stratification and water warming [START_REF] Margirier | Abrupt warming and salinification of intermediate waters interplays with decline of deep convection in the Northwestern Mediterranean Sea[END_REF]. This was concomitant with the increase of winter temperature, NAO and precipitation in the Bay of Marseille, and (i) a later seasonal onset of the 200-1000 µm zooplankton biomass, copepods (calanoids and ergasilida), appendicularians, total zooplankton (ii) an earlier seasonal onset of salps. Therefore, we cannot exclude the hypothesis that winter large-scale forcing impacted the coastal zooplankton dynamics in the Bay of Marseille in the 2010s.

4.3-Implications of the variations in the pelagic ecosystem

The conceptual schema, The results of the Bay of Marseille time series contribute to the documentation of the alterations in the Mediterranean pelagic ecosystem [START_REF] De Madron | Marine ecosystems' responses to climatic and anthropogenic forcings in the Mediterranean[END_REF], particularly changes occurring in the Gulf of Lion.

Espinasse et al. [START_REF] Espinasse | Defining zooplankton habitats in the Gulf of Lion (NW Mediterranean Sea) using size structure and environmental conditions[END_REF], on the basis of a large-scale survey in winter and spring, defined three zooplankton habitats in the Gulf of Lion: the continental shelf, the Rhone influence zone and the Occitan littoral zone. They found that the zooplankton habitat in the BoM presented the same environmental and biological characteristics as most of the Gulf of Lion continental shelf.

Therefore, the Bay of Marseille time series can be considered as a sentinel zone for monitoring environmental effects from the open sea as well as the coastal area. [20] and with observed changes in the diets of small pelagic fishes [START_REF] Bourg | Trophic niche overlap of sprat and commercial small pelagic teleosts in the Gulf of Lions (NW Mediterranean Sea)[END_REF]. This supports the hypothesis of a bottom-up control of the pelagic food chain up to the planktivorous fishes in the Gulf of Lion [START_REF] Saraux | Small pelagic fish dynamics: A review of mechanisms in the Gulf of Lions[END_REF][START_REF] Feuilloley | Concomitant changes in the environment and small pelagic fish community of the Gulf of Lions[END_REF], as shown in other areas in the Mediterranean Sea [START_REF] Macias | Biogeochemical control of marine productivity in the Mediterranean Sea during the last 50 years[END_REF][START_REF] Piroddi | Historical changes of the Mediterranean Sea ecosystem: modelling the role and impact of primary productivity and fisheries changes over time[END_REF].

By taking into account indicators related to the zooplankton populations seasonal onset, we

show that the process of match/mismatch with small pelagic fish [START_REF] Cushing | Plankton Production and Year-class Strength in Fish Populations: an Update of the Match/Mismatch Hypothesis[END_REF] may be a process at work in the Gulf of Lion. An increased mismatch between the spring peak of zooplankton biomass (mainly calanoids) and the growth phase of small pelagic fish [START_REF] Costalago | Feeding of European pilchard (Sardina pilchardus) in the northwestern Mediterranean: from late larvae to adults[END_REF][START_REF] Nikolioudakis | Trophodynamics of anchovy in a non-upwelling system: direct comparison with sardine[END_REF][START_REF] Chen | Temporal variation in prey selection by adult European sardine (Sardina pilchardus) in the NW Mediterranean Sea[END_REF] could certainly explain their body condition. In their experimental study, Queiros et al. [START_REF] Queiros | Is starvation a cause of overmortality of the Mediterranean sardine?[END_REF], showed that sardine can display adaptative phenotypic plasticity to food condition changes. Under lower food quantity and quality conditions, the smallest phenotypes experienced lower mortality by starvation than the larger ones during the critical post-reproductive period (i.e. at the end of winter). In this context, a delay in the availability of these preferred fish prey might favor smaller fish phenotypes. In addition, the zooplankton biochemical composition and energy content may vary seasonally [START_REF] Chen | Seasonal variation in biochemical and energy content of size-fractionated zooplankton in the Bay of Marseille (North-Western Mediterranean Sea)[END_REF][START_REF] Barroeta | Energy density of zooplankton and fish larvae in the southern Catalan Sea (NW Mediterranean)[END_REF] due to the reorganization of different zooplankton assemblages. Further investigations of the interannual variations of the mesozooplankton quality, both in species composition and biochemical content, are needed to improve our understanding of changes in the trophic environment of the small pelagic fishes.

Fig 1 .

 1 Fig 1. Map of the study site. In the left panel, localization of the Bay of Marseille (red) at

  7% of the data variability (Figs 2A and 2B). In each figure, positive and negative influences of the factorial axes on the size spectra shape are displayed. The variation modes of the size spectra shape on the first axis explained 45.9% of the data variability. Positives values on the first factorial axis (Fig 2 A, red curve) the size spectra were associated with a displacement of size spectra mode towards higher ESD. Negatives values on first factorial axis (Fig 2 A, blue curve) was associated with a displacement of size spectra mode toward lower values. The second axis of the FPCA explained 20.8% of the data variability. Positives values on the second PC (Fig 2B, red curve) were associated with the displacement of the size spectra peak toward lower ESD and the increase of the proportion of larger zooplankton organisms (>500 µm ESD). Negative values on the second PC (Fig 2B, blue curve) were associated with the displacement of the size spectra peak toward higher ESD and the decrease of the proportion of larger zooplankton organisms (>500 µm ESD). The MDS (Figs 2C and 2D) computed on the abundances of the thirteen taxa best recognized by the CNN model summarized the structure in species of the zooplankton community (stress <0.2). The ordination method separated in the first axis the centroid of two copepod taxa (calanoids and oithonoids, positive values) to chaetognaths and salps (negative values). The second axis (Fig 2C) separated bivalve, pteropods and harpacticoid centroids (positive values) from all other species centroids. The last axis (Fig 2D) discriminated crustacean centroid (positive values) from copepods nauplii centroid (negative values).

Fig 2 .

 2 Fig 2. Results of the FPCA (A & B) and MDS (C & D) analysis to derive community

  zooplankton community in the Bay of Marseille is dominated by calanoids(Fig 3A) which represent more than 50% of the individuals sampled. The second most represented taxa were appendicularians, oithonoids and crustaceans which represent with calanoids, on average, more than 80% of the zooplankton community identified in this study. Obvious changes in relative abundance of some species occurred over the observation period. The dominance of calanoids seemed to have declined from more than 60% before 2012 to 50% after 2017, in contrast to the relative abundance of crustaceans which almost doubled from 5 to 10% during the same period.The changes in the relative abundances were mainly attributed by the diminution of copepod calanoid and oithonoid absolute abundances (Fig3B). More details concerning the patterns of variation in the zooplankton community structure in species were observable on the interannual series depicted at monthly scale (Fig4). The abundances of bivalves, chaetognaths, harpacticoids, crustaceans, pteropods, salps, appendicularians, fish eggs, and cnidarians seemed to have increased during the study and breakpoints were detected for all these taxa (except for cnidarians) in the last five years of the time series. The series of the two main copepod taxa (calanoids and oithonoids), copepod nauplii and total zooplankton abundance showed a diminution, and a breakpoint was detected for calanoids and oithonoids series at the end of 2012 (in November and September). Abundance estimations of ergasilida presented a breakpoint in April 2007, resulting in an increase in abundance. Interannual variations of the sample scores on the first MDS axis diminished (with a breakpoint detected in April 2014); this observation was based mainly on the diminution of the copepod dominance while abundances of other taxonomic groups such as chaetognaths and salps increased. The second axis highlighted the increase in abundances of bivalves, harpacticoids and pteropods (breakpoint in October 2012). The third axis represented the diminution of nauplii and increase of crustaceans abundance with the breakpoint in April 2017. The smoothing of the first axis of the FPCA revealed that after a displacement toward lower ESD values (between October 2009 -March 2012), the size spectra shifted toward higher ESD after March 2012. No outstanding interannual variations of the second PC were observable and no interannual structural change of the series was detected. The biomass of the 1000 -2000 µm size-fraction presented a breakpoint associated with a decrease in September 2008. The three other zooplankton biomass size fractions decreased later, but no breakpoints were detected.

Fig 3 .

 3 Fig 3. Year-to-year variations in the A) relative and B) absolute abundance of the 13 best

Fig 4 .

 4 Fig 4. Representations of the monthly interannual times series of the zooplankton

Fig 5 .

 5 Fig 5. The two trends depicted by the best DFA model (diagonal and equal R matrix and

Fig 6 .

 6 Fig 6. The two trends depicted by the best DFA model (diagonal and unequal R matrix)

Fig 7 .

 7 Fig 7. Boxplot of the seasonal (monthly) pattern of the zooplankton time series. When a

Fig 8 .

 8 Fig 8. Level plots of the monthly interannual times series of the zooplankton variables.

Fig 9 .

 9 Fig 9. Results of the best DFA model (diagonal and equal R matrix, one trend) performed

Fig 5 )

 5 including changes in the abiotic conditions and in the lower trophic levels (micro-nano-pico-plankton). We described multiple environmental trends:-Higher values of salinity and oxygen between 2005 and 2007 (Fig 5, trend 1). This anomaly can be interpreted as the signature of stronger offshore deep-water inflows during this period. By analyzing satellite images in those years, Mayot et al.[START_REF] Mayot | Influence of the Phytoplankton Community Structure on the Spring and Annual Primary Production in the Northwestern Mediterranean Sea[END_REF] highlighted particularly high phytoplankton productivity induced by strong deep water convection events occurring over a broad area in the Provencal basin with effects extending to surrounding coastal areas, including the Bay of Marseille.-The reduction in nutrients and particulate matter inputs observed in our time series (Fig 5, trend 1) is consistent with observations already made in the Bay of Marseille, and appears to be a general trend in French Mediterranean coastal waters[START_REF] Lheureux | Bi-decadal variability in physicobiogeochemical characteristics of temperate coastal ecosystems: from large-scale to local drivers[END_REF][START_REF] Goberville | Climate-driven changes in coastal marine systems of western Europe[END_REF]. The physical and chemical characteristics of the Bay of Marseille might have been impacted by the diminution of the Rhône nutrient inputs in the Gulf of Lion[START_REF] Feuilloley | Concomitant changes in the environment and small pelagic fish community of the Gulf of Lions[END_REF] and/or by the proximity to the Marseille sewage treatment plant where the treatment efficiency has improved since 2008[START_REF] Millet | Unexpected spatial impact of treatment plant discharges induced by episodic hydrodynamic events: Modelling Lagrangian transport of fine particles by Northern Current intrusions in the bays of Marseille (France)[END_REF][START_REF] Raimbault | Chapitre 7 : Le milieu marin autour de Marseille[END_REF]. While ammonium and phosphate decreased, nitrate concentration was stabilized (and even presented a slight increase in the second part of the study, see S3 file).

  affected their distribution. Our observations evidenced changes in size and abundance of the micro-, nano-, pico-plankton in 2013 after the strong decrease in nutrient concentrations which were concomitant with changes in the mesozooplankton community structure (Fig 5B, Fig 6B and Fig 6C)

Fig 10 ,

 10 summarizes the interannual changes analyzed in the present paper concerning the environment of the Bay of Marseille and its mesozooplankton community.

Fig 10 .

 10 Fig 10. Summary of our main results by means of the revisited Ramon Margalef mandala

  Fig. 2 764

Fig. 4

 4 Fig. 4

Fig. 5

 5 Fig.5

Fig. 6

 6 Fig. 6

Fig. 7

 7 Fig. 7

Fig. 8

 8 Fig. 8

Fig. 9

 9 Fig. 9

Fig. 10

 10 Fig. 10

Table 1 : Summary of the variables used in this study. Variable Abbreviation Unit Category Time Source Wind stress (according to NW-SE axis)

 1

		Wind stress	-	Environmental	2005-2020	Calculation from Météo France
	Precipitation	Prec	mm	Environmental	2005-2020	Météo France
	Mixed Layer Depth	MLD	m	Environmental	2005-2020	Calculation from SOMLIT data
	Western Mediterranean Oscillation	WeMO	-	Environmental	2005-2020	[34]
	Northern Atlantic Oscillation	NAO	-	Environmental	2005-2020	NOAA
	Temperature	T	°C	Environmental	2005-2020	SOMLIT
	Salinity	S	-	Environmental	2005-2020	SOMLIT

100 2.2-Zooplankton data.

	Picophytoplankton counts (subsurface)	Pico	cells.mL -1	Environmental	2009-2020	SOMLIT
	Picophytoplankton size index (subsurface)	Pico size	-	Environmental	2009-2020	SOMLIT
	Nanoplankton counts (subsurface)	Nano	cells.mL -1	Environmental	2009-2020	SOMLIT
	Nanoplankton size index (subsurface)	Nano size	-	Environmental	2009-2020	SOMLIT
	Biomass for size fraction 1000-2000µm	Biom 1000-2000	mg of dry weight.m -3	Zooplankton	2005-2020	This study
	Biomass for size fraction 500-1000µm	Biom 500-1000	mg of dry weight.m -3	Zooplankton	2005-2020	This study
	Biomass for size fraction 300-500µm	Biom 300-500	mg of dry weight.m -3	Zooplankton	2005-2020	This study
	Biomass for size fraction 200-300µm	Biom 200-300	mg of dry weight.m -3	Zooplankton	2005-2020	This study
		Total Biomass 200-2000 µm	Total Biomass	mg of dry weight.m -3	Zooplankton	2005-2020	This study
	Zooplankton size structure (1 st PC)	Zoopk PC1	-	Zooplankton	2005-2020	This study
	Zooplankton size structure (2nd PC)	Zoopk PC2	-	Zooplankton	2005-2020	This study
	Community structure index MDS 1	MDS1	-	Zooplankton	2005-2020	This study
	Community structure index MDS 2	MDS2	-	Zooplankton	2005-2020	This study
	Community structure index MDS 3	MDS3	-	Zooplankton	2005-2020	This study
	Zooplankton abundances	Abundance	individuals.m -3	Zooplankton	2005-2020	This study
	101
	102	In parallel with the SOMLIT monitoring survey, mesozooplankton samples have been collected
	103	since 2005 through vertical hauls using a WP2 200 µm mesh size plankton net (water column
	104	sampled: 0-55 m). The samples were fractioned in two. The first half of the cod-end content
	105	was preserved in a 4% buffered formaldehyde solution for digitalizing; the other half was
	106	maintained in cold condition for zooplankton size-fractioned dry weight (biomass)
	107	measurements. Complementary information concerning the environmental and zooplankton
	108	monitoring is given in the S1 file.

Acknowledgements

The zooplankton time series was initiated in 2002 by F. Carlotti and was linked with the ongoing SOMLIT monitoring program. The authors acknowledge all participants in the Bay of Marseille SOMLIT monitoring program: Patrick Raimbault as PI of the program, Michel Lafont for the regular sampling of physical and biogeochemical parameters, and the crew of the R.V. Antedon II. The regular sampling of zooplankton was carried out before 2012 mainly by F.Carlotti's students as part of their PhD or MSc work, V. Riandey, E. Dieval, A. Nowaczyk, K. Morsly, N. Neffati, A. Dron, S. Martini, then after 2012, mainly by L. Guilloux (Research Engineer CNRS), sometimes supported by PhD students K. Donoso, C-T Chen and T Garcia.

We would like to thank all of them for their contribution to this sampling effort. We thank P. Raimbault, N. Garcia and M. Lafont for making available the physical and biogeochemical data, and V. Cornet, S. Boussabat and K. Leblanc for processing the phytoplankton samples whose data are used. Chemical and biological sample treatments benefited from the Mediterranean Institute of Oceanography platforms (PACEM for chemical analyses and MIM

Thanks are also addressed to Michael Paul, native speaker, for English proofreading.