Coralie Picoche

Laure Carassou

Isabelle Auby

F Barraquand
email: frederic.barraquand@u-bordeaux.fr

D Maurer

A political sociology of ECOsystem sciences: theories, narratives, interactions and GOVernance View project Gironde estuary Exotic species View project

Keywords: competition, diatoms, dinoflagellates, hydrodynamics, interactions, MAR, mutualism, network, population dynamics, seasonality, time series

Phytoplanktonic communities maintain a high diversity in a seemingly homogeneous environment, competing for the same set of resources. Many theories have been proposed to explain this coexistence despite likely competition, such as contrasted responses to temporal environmental variation. However, theory has developed at a faster pace than empirical evaluation using field data, which requires inferring drivers of community dynamics from observational time series. Here, we combine autoregressive models with a data set spanning more than 20 years of phytoplankton counts every two weeks, together with nutrients and physical variables. By comparing models dominated by nutrients or physical variables (hydrodynamics and climate), we first explore which abiotic factors contribute more to phytoplankton growth and decline. We find that physical drivers -such as irradiance, wind, and salinity -explain some of the variability in abundances unexplained by biotic interactions. In contrast, responses to nutrients explain less of the phytoplankton variability. Concerning biotic drivers of community dynamics, multivariate autoregressive models reveal that competition between different groups (genera) has a much weaker effect on population growth rates than competition within a group. In fact, the few biotic interactions between genera that are detected are frequently positive. Hence, our system is unlikely to be best represented as a set of competitors whose differing responses to fluctuating environments allow coexistence, as in "paradox of the plankton" models with a storage effect or a relative nonlinearity of competition. Coexistence is more likely to result from high intragroup density-dependence. Competition between phytoplanktonic groups and nutrient limitation are often invoked as drivers of phytoplankton dynamics; our findings suggest instead that more attention should be given to the physical structure of the environment and natural enemies such as grazers and pathogens, that could explain the high intraspecific density dependence found here.

Introduction

What maintains the diversity of species-rich communities is an old yet still challenging question for ecological theory [START_REF] Hutchinson | The paradox of the plankton[END_REF][START_REF] Jewson | Coexistence of two Cyclotella diatom species in the plankton of Lake Baikal[END_REF][START_REF] Li | The effects of dynamical rates on species coexistence in a variable environment: the paradox of the plankton revisited[END_REF]. The continued coexistence of phytoplanktonic taxa is amongst the most puzzling diversity patterns: species, genera and even classes of plankton are sometimes competing for the same limited resources in a seemingly homogeneous environment [START_REF] Titman | Ecological competition between algae: experimental confirmation of resource-based competition theory[END_REF][START_REF] Tilman | Phytoplankton community ecology: the role of limiting nutrients[END_REF]. This state of affairs led [START_REF] Hutchinson | The paradox of the plankton[END_REF] to refer to the "paradox of the plankton", a paradox for which many theoretical answers have been proposed [START_REF] Record | The paradox of the "paradox of the plankton[END_REF]. Classical mechanisms for biodiversity maintenance, such as niche separation or spatial variation in the environment, are not seen as obvious in the case of plankton, because species apparently compete for the same set of resources in environments that can be rather well-mixed in the absence of strong stratification of the water column (Huisman et al. 1999b).

Neutrality, i.e., per capita equivalence of birth and death rates, is another potential explanation [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32)[END_REF]. There are varied life-histories in phytoplankton (Litchman and Klausmeier 2008), but neutrality does not require organisms to be equal, just that trade-offs in life history traits equalize their fitnesses [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32)[END_REF][START_REF] Doncaster | Ecological equivalence: a realistic assumption for niche theory as a testable alternative to neutral theory[END_REF]). However, experimental work does suggest a difference in net reproduction rates among plankton species under varying nutrient concentrations [START_REF] Tilman | Phytoplankton community ecology: the role of limiting nutrients[END_REF], which means that birth and death rates are in fact likely to vary a lot between species, depending on environmental conditions.

Another leading hypothesis relates to the temporal variation in the environment [START_REF] Hutchinson | The paradox of the plankton[END_REF][START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF][START_REF] Litchman | Competition of phytoplankton under fluctuating light[END_REF][START_REF] Li | The effects of dynamical rates on species coexistence in a variable environment: the paradox of the plankton revisited[END_REF]. There is often a rather strong seasonal variation in the environment experienced by plankton (e.g., temperature and irradiance). A number of abiotic variables have some extra temporal variation as well, often due to perturbations from nearby ecosystems, especially in coastal systems like estuaries and lagoons (e.g., nitrogen (N), phosphorus (P), silicon (Si), due to terrestrial inflows). Theory posits that diversity maintenance can be precisely due to this temporal variation in the environment [START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF][START_REF] Li | The effects of dynamical rates on species coexistence in a variable environment: the paradox of the plankton revisited[END_REF], as in models of the storage effect [START_REF] Chesson | The roles of harsh and fluctuating conditions in the dynamics of ecological communities[END_REF]. Such models are different from neutral models, as species have differentiated responses to environmental variables. Some near-neutral models do include temporal responses to the environment [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF], but near-neutral models usually assume that species compete in a zero-sum game, where no species can increase without a decrease in some other species (but see [START_REF] Jabot | Non-random correlation of species dynamics in tropical tree communities[END_REF]. However, total biomass across species/genera often fluctuates across several orders of magnitude throughout the year, blooms being the most obvious defining feature of planktonic dynamics. This makes phytoplanktonic communities fundamentally different from other species-rich communities, such as forests or coral reefs, where the fluctuations in biomass are clearly milder and near-neutral dynamics may be more plausible [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32)[END_REF]; but see [START_REF] Ser-Giacomi | Ubiquitous abundance distribution of non-dominant plankton across the global ocean[END_REF] for rare plankton species).

Temporal variation is therefore a likely cause of coexistence in plankton -as initially proposed by [START_REF] Hutchinson | The paradox of the plankton[END_REF]. In this mechanism, all species are assumed to outcompete others at some point in time, with alternation in the ranking of competitors preventing competitive exclusion [START_REF] Descamps-Julien | Stable coexistence in a fluctuating environment: an experimental demonstration[END_REF]. This mechanism is likely to work if species have differing nonlinear responses to the environment (i.e., relative nonlinearity of competition), or a covariance between environmental conditions and competitive strength, as in the storage effect [START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF][START_REF] Fox | The intermediate disturbance hypothesis should be abandoned[END_REF][START_REF] Li | The effects of dynamical rates on species coexistence in a variable environment: the paradox of the plankton revisited[END_REF]. The plausibility of a fluctuation-driven coexistence mechanism has been empirically verified on diatoms that belong to different genera (e.g., [START_REF] Descamps-Julien | Stable coexistence in a fluctuating environment: an experimental demonstration[END_REF] though some experiments and models also suggest that planktonic communities might be in a chaotic or seasonally-driven chaotic state [START_REF] Huisman | Biological conditions for oscillations and chaos generated by multispecies competition[END_REF][START_REF] Benincà | Chaos in a long-term experiment with a plankton community[END_REF][START_REF] Dakos | Interannual variability in species composition explained as seasonally entrained chaos[END_REF]. Such a chaotic state invokes other fluctuation-driven coexistence mechanisms, with nonlinearities in functional forms promoting endogeneously-driven fluctuations. Of course, these mechanisms are not exclusive, and coexistence is generally due to the joint influence of equalizing factors (i.e., species have the same fitness when averaged over the variation in the environmental variables) and stabilizing mechanisms (i.e., species increases when rare) that make coexistence all the more likely [START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF]. Temporal partitioning of resources, highlighted above, is not the only way to produce coexistence through niche-based mechanisms here, as indirect interaction pathways can increase intraspecific interaction strength. These include most notably natural enemies such as predators [START_REF] Haydon | Pivotal assumptions determining the relationship between stability and complexity: an analytical synthesis of the stability-complexity debate[END_REF][START_REF] Vallina | Maximal feeding with active prey-switching: A killthe-winner functional response and its effect on global diversity and biogeography[END_REF][START_REF] Barabás | Emergent neutrality or hidden niches?[END_REF] or parasites, as well as toxin production [START_REF] Felpeto | Allelopathy prevents competitive exclusion and promotes phytoplankton biodiversity[END_REF]. In all of these mechanisms, some degree of temporal fluctuation in phytoplankton abundances is expected, even though not necessarily the main force driving the coexistence as in fluctuating-environment models.

Therefore, any framework examining competition among planktonic species should account for both the broad variations in environmental conditions and in planktonic abundances over time. Multivariate autoregressive (MAR) modeling is one such dynamic framework that has been increasingly used to examine interactions between planktonic groups in a dynamic environment [START_REF] Klug | Compensatory dynamics in planktonic community responses to pH perturbations[END_REF][START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF]Hampton and Schindler 2006;[START_REF] Huber | The role of predation for seasonal variability patterns among phytoplankton and ciliates[END_REF][START_REF] Scheef | Inferring plankton community structure from marine and freshwater long-term data using multivariate autoregressive models[END_REF][START_REF] Griffiths | Phytoplankton community interactions and environmental sensitivity in coastal and offshore habitats[END_REF][START_REF] Gsell | Quantifying change in pelagic plankton network stability and topology based on empirical long-term data[END_REF]. MAR models enable the estimation of interaction strengths between taxa [START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF][START_REF] Mutshinda | What drives community dynamics?[END_REF] as well as the dependence of population growth rates on abiotic variables [START_REF] Hampton | Quantifying effects of abiotic and biotic drivers on community dynamics with multivariate autoregressive (MAR) models[END_REF], both necessary to model planktonic dynamics. These models are linear on a log-scale, hence they represent population growth processes as multiplicative, power-law functions of densities that can approximate more complex nonlinear functional forms [START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF]. For example, for 3 species or groups, the equation for species 1 dynamics is N t+1,1 = e a 1 N b 11 t,1 N b 12 t,2 N b 13 t,3 e 1t with N t,i the species counts and b ij the interaction coefficients (see Material and methods). One can also consider an increasing degree of nonlinearity using phase-based MAR models (also called "Threshold AR" models, these are piecewise log-linear, [START_REF] Stenseth | Testing for effects of climate change on competitive relationships and coexistence between two bird species[END_REF], in order to examine in detail the plausibility of nonlinear coexistence mechanisms.

Most plankton-based MAR analyses so far (but see [START_REF] Huber | The role of predation for seasonal variability patterns among phytoplankton and ciliates[END_REF] have aggregated data at the class level, therefore preventing any attempt to examine the strengths of competitive interactions among species or genera. Recent papers suggest that a finer taxonomic resolution would allow for a better understanding of community dynamics [START_REF] Griffiths | Phytoplankton community interactions and environmental sensitivity in coastal and offshore habitats[END_REF][START_REF] D'alelio | Ecological-network models link diversity, structure and function in the plankton food-web[END_REF]. Here, we take advantage of a long-term dataset of coastal phytoplankton monitoring (>20 years of counts every two weeks) with fine taxonomic resolution -to the genus level -to investigate:

(1) What are the abiotic drivers of population dynamics? Nutrients are classically assumed to be key drivers of phytoplankton dynamics. How do their influences compare to other physical drivers such as light and hydrodynamics?

(2) What are the strengths of interactions between and within phytoplankton groups? And what can interaction strengths tell us about biotic interactions?

(3) What are the implications of our results, obtained on a rich diatom and dinoflagellate assemblage, for theory on phytoplankton diversity maintenance?

Material and methods

Study area and sampling details

Arcachon Bay (AB) is a 155 km 2 coastal lagoon located in the south-west of France (Fig. 1). AB is connected to the Atlantic Ocean by two main channels (41 km 2 ), and freshwater inflow comes from the Eyre river (83%) and Porge canal (11%), as well as rainfall and groundwater (6%, [START_REF] Rimmelin | Study of annual and seasonal dissolved inorganic nitrogen inputs into the Arcachon Lagoon, Atlantic Coast (France). -Estuarine[END_REF]. Using a 2D-hydrodynamics model, [START_REF] Plus | Hydrodynamic characterization of the Arcachon Bay, using model-derived descriptors[END_REF] estimated AB's tidal prism, i.e. the volume of water leaving the bay at ebb tide, around 384 million m 3 . This volume can be compared to the much lower maximum daily water discharge of 12 million m 3 from the Eyre, the main river entering the bay. [START_REF] Plus | Hydrodynamic characterization of the Arcachon Bay, using model-derived descriptors[END_REF] also computed a flushing time for the lagoon between 13 and 30 days and highlighted a high return flow factor (fraction between 0.94 and 0.95, 1 being the maximum value), meaning that a compound carried out of the lagoon is likely to return to it. The environment is therefore well-mixed at large spatial scales, and mostly marine-influenced. In addition, the wind regime has a strong influence on the hydrodynamics of the bay: westerly winds represent about 77% of total wind and can increase water residence time in the bay up to 10 days, leading to surges up to 50 cm [START_REF] Rimmelin | Study of annual and seasonal dissolved inorganic nitrogen inputs into the Arcachon Lagoon, Atlantic Coast (France). -Estuarine[END_REF][START_REF] Plus | Hydrodynamic characterization of the Arcachon Bay, using model-derived descriptors[END_REF][START_REF] Plus | Phytoplankton versus macrophyte contribution to primary production and biogeochemical cycles of a coastal mesotidal system. A modelling approach[END_REF]. The oceanic temperate climate is characterized by an amount of precipitation averaging 785 mm year -1 , a total irradiance around 475 kJ cm -2 year -1 and an average water temperature around 15.6°C, with a strong seasonality for all parameters (rainfall is 40% less abundant in summer than in winter, irradiance is four times higher, and temperature varies between 9°C and 23°C). More details on AB (e.g., absence of stratification of the water column) are given in Supplementary material Appendix 1: Section A1.1.1.

Water samples were collected at 2 sites: Teychan (44°40'25" N / 1°09'31" W, water depth = 17.30m) and Buoy 7 (B7, 44°32'32" N / 1°15'49" W, water depth = 14.5m), from 1987 to 2015 (709 dates) and 2003 to 2015 (311 dates), respectively. For more information on water sampling, see Supplementary material Appendix 1: Section A1.1.3. Teychan is located near the limit drawn by [START_REF] Bouchet | Stratifications, fronts halins dans une lagune mésotidale (Bassin d'Arcachon-France). -In: 3e Colloque International d'Océanographie du Golfe de Gascogne[END_REF] between water masses influenced by oceanic waters or by continental inputs, in the middle of AB (Fig. 1), while Buoy 7 is at the oceanic entrance of the channel. Teychan is therefore expected to be more sensitive to continental inputs while Buoy 7 should be mainly subject to marine influences [START_REF] Bouchet | Stratifications, fronts halins dans une lagune mésotidale (Bassin d'Arcachon-France). -In: 3e Colloque International d'Océanographie du Golfe de Gascogne[END_REF][START_REF] Glé | Typology of environmental conditions at the onset of winter phytoplankton blooms in a shallow macrotidal coastal ecosystem, Arcachon Bay (France)[END_REF].

Phytoplankton data

The National Phytoplankton and Phycotoxin Monitoring Network (REPHY http://envlit. ifremer.fr/infos/rephy_info_toxines) involves phytoplankton collection every two weeks, within two hours of high tide (for more details, see Supplementary material Appendix 1: Section A1.1.3). Data collection started in 1987 in AB (Teychan), so the full dataset consists of 709 sampling dates. The full dataset is used in all analyses except for multivariate autoregressive models including all phytoplanktonic groups (see "MAR(1) models"). As cryptophytes were not recorded properly before 1996, and this could affect the counts of other taxa, we focused on observations post-1996 for MAR(1) models. Taxonomic units were aggregated at the genus level following previous work on this and other plankton datasets [START_REF] Hernandez | Evaluation of the taxonomic sufficiency approach for ichthyoplankton community analysis[END_REF][START_REF] Hernández Fariñas | Assessing phytoplankton realized niches using a French national phytoplankton monitoring network. -Estuarine[END_REF], according to the plankton experts' knowledge (see Acknowledgements). We made only two exceptions to this "genus rule": cryptophytes and euglenophytes could not be consistently identified below the class level, but their abundance throughout the time series justified retaining them in the analyses.

Our analyses therefore focused on the best-resolved and most abundant planktonic groups over time (Table 1), in the microphytoplankton fraction (above 20 µm, [START_REF] Reynolds | The ecology of phytoplankton[END_REF]. Hereafter, we use the word "population" for groups described in Table 1. The term "community" encompasses all studied planktonic taxa potentially interacting in AB.

Phytoplankton abundance values (cell L -1 ) were linearly interpolated over a regular time sequence with a 14-day time step to correct for some sampling irregularity (following Hampton et al. 2006). We evaluated several methods for reconstructing time series without gaps and they are detailed in Supplementary material Appendix 1: Section A1.2.1; we followed Hampton and Schindler (2006).

Abiotic data

Environmental variables collected at Teychan and Buoy 7 included water temperature (TEMP in °C), salinity (SAL, in g kg -1 ), nutrients (ammonia NH + 4 , silicates Si(OH) 4 , phosphates PO 3- 4 and nitrates-nitrites NO x =NO 2 + NO - 3 , all in µmol L -1 ) and suspended particulate matter (SPM in mg L -1 ), with an organic part assessment (SPOM, in mg L -1 ). They were measured following Aminot andKérouel (2004, 2007) (see Supplementary material Appendix 1: Section A1.1.3 for more details). Some abiotic variables (nutrients and SPM) could not be recorded at Teychan until 2007.

Table 1: Name and composition of the phytoplanktonic groups used in the paper, with their average proportion in the phytoplanktonic community, calculated as the ratio of their summed abundance over the summed abundance of all identified phytoplanktonic organisms, and the frequency of detection, calculated as the ratio of samples in which the group was present over the total number of samples Daily meteorological data -rainfall (mm), irradiance (J cm -2 ), wind direction and velocity (m s -1 ) -were provided by Météo France for the nearby Cap Ferret station (44°37' N / 1°14' W) and were used as inputs for both sampling stations. Wind energy was extracted from these data as the squared velocity (m 2 s -2 ).

Daily North Atlantic Oscillation (NAO) index was downloaded from the National Oceanic and Atmospheric Administration (NOAA) Weather Service for Climate Prediction Centre website. Monthly Atlantic Multidecadal Oscillation (AMO) index was downloaded from the NOAA Physical Science Division website.

All environmental data were linearly interpolated on the same dates as phytoplankton sampling dates, removing missing values (see Supplementary material Appendix 1: Fig. A1.3).

We also transformed abiotic variables to represent meaningful biological and physical processes. We summed NH + 4 and NO x to obtain the total inorganic nitrogen availability (Ntot). We integrated inflow, rainfall and irradiance between t -1 and t to represent the growth conditions of phytoplankton groups between two sampling dates [START_REF] Glé | Typology of environmental conditions at the onset of winter phytoplankton blooms in a shallow macrotidal coastal ecosystem, Arcachon Bay (France)[END_REF]. NAO was averaged over the same period and a variable summarizing the effect of wind was constructed using the mean wind energy over 3 days (F. Ganthy, Ifremer, personal communication).

Environmental drivers of population dynamics

The high number of possible explanatory variables with some degree of auto-correlation and crosscorrelation (see Supplementary material Appendix 1, Figs. A1.4 -A1.7) led to a high collinearity indicated by a very high condition index (CI=150, well above the thresholds of 10 or 30 indicated by [START_REF] Belsley | A guide to using the collinearity diagnostics[END_REF]. To select relevant variables, our approach was three-fold: (a) spectral analyses enabled us to detect significant coherency (i.e., correlation in the spectral domain) between biotic and abiotic cycles to select variables that seemed coupled to phytoplankton abundance, (b) nutrients being highly correlated (Supplementary material Appendix 1: Fig. A1.4-Fig. A1.7), we selected the most relevant ones based on literature values on phytoplanktonic requirements, which then led us to (c) compare a nutrient-limitation model using N and P as covariates (including the possibility for saturating functions of nutrients) to a physics-driven model using irradiance, salinity and wind energy as covariates. This comparison was based on log-linear, single species autoregressive modeling. A "full model" including both types of drivers was also compared to the aforementioned more parcimonious models. Climate indices (NAO and AMO) were included in a monthly model but explained little of the observed variation, in line with previous results for AB [START_REF] David | Spatial and long-term changes in the functional and structural phytoplankton communities along the French Atlantic coast. -Estuarine[END_REF]. As they were also more difficult to interpret ecologically, they were discarded from further analyses.

Spectral analyses

Phytoplankton groups show a periodicity that can be related to seasonal cycles. Correlating biotic and abiotic time series at different frequencies therefore helps to extract planktonic responses to environmental variability. For this, we compared each genus periodogram, using a modified Daniell kernel with a smoothing window of 2 time steps (1 month) and no tapering, with abiotic variables. Coherence significance was determined with a corrected (Bonferroni) 5% threshold.

Nutrient saturation

Nutrient uptake is a highly variable process both within and between phytoplankton genera [START_REF] Paasche | Silicon and the ecology of marine plankton diatoms. II. Silicate-uptake kinetics in five diatom species[END_REF][START_REF] Fisher | Phytoplankton, nutrients and turbidity in the Cheesapeake, Delaware, and Hudson Estuaries. -Estuarine[END_REF][START_REF] Sarthou | Growth physiology and fate of diatoms in the ocean: a review[END_REF][START_REF] Reynolds | The ecology of phytoplankton[END_REF][START_REF] Litchman | The role of functional traits and trade-offs in structuring phytoplankton communities: scaling from cellular to ecosystem level[END_REF]). To assess potential limitations to growth and insert ecologically meaningful predictors into autoregressive models (Fig. 2), we evaluated whether nutrient uptake was possibly saturated; details can be found in Supplementary material Appendix 1 (section A1.1.8). Nutrient uptake functions were then used as covariates in some of the linear autoregressive models. Single-species autoregressive models Before using multivariate autoregressive (MAR) models, each phytoplanktonic group was studied on its own. The modeling philosophy adopted here classically assumes that a multidimensional dynamical system might be represented as a delayed unidimensional dynamical system [START_REF] Takens | Detecting strange attractors in turbulence[END_REF][START_REF] Turchin | Complex Population Dynamics: A Theoretical/Empirical Synthesis[END_REF], even in a stochastic context [START_REF] Abbott | Environmental variation in ecological communities and inferences from single-species data[END_REF]. Accordingly, each population can be described by an autoregressive (AR) model where the number of time lags depends on the dimensionality of the system. The goal was to identify potential abiotic drivers of population dynamics before the multivariate analyses, in order to fit later on more parsimonious MAR models (Fig. 2; bear in mind that a full-matrix MAR(1) model with covariates can require to estimate up to 204 parameters here).

Lag order selection For log-linear AR modeling, focusing on one species at a time, we assumed that phytoplankton ecology did not vary greatly between sites and concatenated data from both stations to increase the effective sample size in the statistical analyses. We fitted a log-linear AR(p) model, which allowed the determination of the maximum lag order p. The shape of growthabundances curves (see Supplementary material Appendix 2: Fig. A2.1) and analyses of residuals (Supplementary material Appendix 2: Section A2.1.4) confirmed that log-linearity on a log-scale -Gompertz rather than Ricker growth [START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF] -was an appropriate approximation. We used the arima function from the stats package [START_REF] Venables | the R Core Team. -An Introduction to R. Notes on R: A Programming Environment for Data Analysis and Graphics[END_REF] in R (version 3.3.2) to estimate coefficients and the corresponding Akaike Information Criterion (AIC) for each AR(p) model. The number of time lags considered varied between 1 and 30, the maximum reasonable value for the length of our dataset (780 dates for taxa monitored after 1996). The likelihood was estimated using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization method, an approximation of Newton's method, with a maximum of 100 iterations.

Adding environmental variables AR(p) models with abiotic predictors affecting growth rates -sometimes called ARX(p) -were then used to study the environmental effects on phytoplanktonic growth, according to equation (1).

µ s,t = α s + p-1 i=0 β s,i n s,t-i + V j=1 γ s,j u j,t + s,t (1)
where µ s,t = n s,t+1 -n s,t is the population growth rate on a log scale for genus s, α s is the genus-specific productivity when population size is zero on a log scale, β s,i is an AR coefficient that can represent intragroup density-dependence at lag i, u j,t is environmental variable j for time t, with effect γ s,j on the population growth rate, and s,t is a random variable modeling both environmental and observation error with mean 0 and variance σ 2 [START_REF] Dennis | Estimating density dependence, process noise, and observation error[END_REF].

All abiotic variables were standardized in order to make their effects comparable. As Si and N were highly correlated (Supplementary material Appendix 1: Fig. A1.5, A1.6; a strong laggedcorrelation also appeared in Fig. A1.7(d)) and Si seemed less limiting than N (Table 2), P and N were the only nutrients retained in linear ARX(3) models (see Supplementary material Appendix 2 for the lag chosen, section A2.1.3). Two models were used to assess N availability: the first one used raw concentrations while the other one used a saturation function (N t = 1 -e -λNtot , with λ determined by N half-saturation constant) to represent limits to nutrient uptake. We also considered models using nutrient ratios (Si/N and P/N) instead of raw concentrations. All abiotic variables were then transformed to take into account seasonality. First, a seasonal component was computed from a linear regression of temperature data against trigonometric functions with an annual frequency (Supplementary material Appendix 1: Section A1.1.9). Abiotic variables were then regressed against this seasonal component and residuals from this regression were added to the model as abiotic variables (see Eq. 2).

T t = a T + S 1 cos(ωt) + S 2 sin(ωt) + T,t X v,t = a v + b v (S 1 cos(ωt) + S 2 sin(ωt)) + η v,t (2)
where T t is the temperature at time t, S 1 cos(ωt) + S 2 sin(ωt) is the seasonal component with a normal noise T,t and η v,t is the residual component of covariable X v,t when season has been taken into account, a v and b v being associated regression coefficients.

Model evaluation was performed in two steps. We first used physics-based variables (salinity, irradiance and wind energy) and nutrients (nitrogen and phosphorus) together to define a "full" model including all covariables. Three representations of nutrient variables were considered and compared: raw concentrations, saturating functions of concentrations -see above "Nutrient saturation" -or nutrient concentration ratios. We considered, for each model, versions with and without a dedicated seasonal component. This led to six different full-model formulations. The minimum AICc for a majority of phytoplankton groups determined the best full model representation. The full model was then compared to its "nutrients-only" and "physics-only" counterparts (i.e., models with only nutrient or physical variables as predictors, respectively), in order to choose the minimum amount of abiotic predictors for the following MAR(1) analysis.

MAR(1) models

Multivariate autoregressive (MAR) models are multivariate extensions of AR models. Instead of focusing on one phytoplanktonic group at a time, these models describe the change in abundance of all groups at the same time, conditioned by past values of both the V environmental variables and the abundances of the S other phytoplanktonic groups. We used [START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF] formulation (Eq. 3), which explains growth between times t and t + 1 by the abiotic variables at time t + 1 (because of the rapid growth of phytoplankton, using variables both at t and t + 1 could be valid choices, see [START_REF] Hampton | Quantifying effects of abiotic and biotic drivers on community dynamics with multivariate autoregressive (MAR) models[END_REF] for an alternative). This formulation was chosen because it (1) led to smaller AICc and Bayesian Information Criterion (BIC) for otherwise identical models and (2) resulted in better consistency between estimates at Teychan and Buoy 7 sites (77% of covariate effects are qualitatively similar). This choice does not impact qualitatively the results (see Supplementary material Appendix 3). The MAR(1) model is described by Eq. ( 3)

n t+1 = Bn t + Cu t+1 + e t , e t ∼ N S (0, Q) (3)
where n t is a S × 1 log abundance vector of phytoplanktonic genera at time t, B is a S × S community (interaction) matrix. When using such MAR(1) model, we need to substract one from autoregressive, intragroup coefficients in order to make them comparable to intergroup interaction values: no intragroup regulation leads to b ii = 1 whereas a strong intragroup effect is characterized by b ii close to or below zero (see [START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF] for details). For i = j, b ij describes the effect of genus j on genus i. C is a S × V environment matrix describing the effects of variables u on phytoplanktonic growth. Therefore c ij is the effect of variable j on genus i. e t is a noise vector which covers both process and observation error (Supplementary material Appendix 1: Section A1.2.1, describes the handling of missing or incomplete observations), following a multivariate normal distribution with a variance-covariance matrix Q. In the following, Q was chosen diagonal, which was more parsimonious and did not change qualitatively the results (see Supplementary material Appendix 4). We used the MARSS package [START_REF] Holmes | MARSS: Multivariate autoregressive state-space models for analyzing time-series data[END_REF][START_REF] Holmes | MARSS: Multivariate autoregressive state-space models for analyzing time-series data[END_REF] to fit MAR(1) models. The likelihood optimization procedure is fully described in Supplementary material Appendix 2, section A2.2.2.

All time series were centered and standardized. We used a physics-only model with a dedicated seasonal component to represent phytoplankton dynamics, based on previous single-species analyses.

Five interaction matrices were considered. Except for the null model which allowed no net interactions between groups, all interaction matrices included the only documented effect in the studied phytoplanktonic community: possible predation between dinoflagellates and cryptophytes [START_REF] Moeller | Acquired phototrophy stabilises coexistence and shapes intrinsic dynamics of an intraguild predator and its prey[END_REF]. Because all other interactions were unknown, we considered several contrasted interaction scenarios for model fitting. Eq. (4) defines for example a scenario differentiating between pennate and centric diatoms. Three other interaction matrices are presented in Supplementary material Appendix 2: Section A2.2.3.

B =                                AST b 1,1 NIT b 1,2 PSE b 1,3 SKE 0 CHA 0 GUI 0 LEP 0 RHI 0 GYM 0 PRP 0 CRY 0 EUG 0 b 2,1 b 2,2 b 2,3 0 0 0 0 0 0 0 0 0 b 3,1 b 3,2 b 3,3 0 0 0 0 0 0 0 0 0 0 0 0 b 4,
0 0 0 0 0 0 0 0 0 0 0 b 12,12                                AST NIT PSE        Pennate diatom SKE CHA GUI LEP RHI                    Centric diatom GYM PRP    Dinoflagellate CRY EUG (4) 0
In our first interaction scenario, "the full model", the interactions were totally unconstrained. The second and third scenarios only took into account intra-phylum competition or facilitation (i.e., within diatoms or dinoflagellates). The third interaction scenario added a further difference between pennate and centric diatoms (eq. 4). Finally, the fourth model only took into account interphylum interactions, which ecologically represents a scenario where, within a phylum, organisms are too close to each other to compete. All models include intra-group interactions that are needed to maintain model stability. For each interaction scenario and model fit, 95% confidence intervals (CIs) were calculated for all parameters by bootstrapping. AICc was used to select models, and this model selection was compared to results obtained by BIC and cross-validation.

We thoroughly tested our model selection approach using 10 simulated datasets resembling our empirical dataset, for 3 contrasted ecological settings: (1) effects of environment only, (2) effects of interactions only, and (3) both environment and interactions (see Supplementary material Appendix 2: Section A2.2.1). The results on simulated data showed that AICc selected the appropriate model in each ecological setting, thus there was no need for more complex model selection using AICb or refined procedures as did [START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF] (results are provided in Supplementary material Appendix 2: Section A2.2.1).

Non-linear dynamics: bivariate phase-dependent (TAR) models

The log-linear MAR(1) framework essentially assumes power-law functional forms on the untransformed scale. For example, for three groups, the equation for group 1 is N t+1,1 = e a 1 N b 11 t,1 N b 12 t,2 N b 13 t,3 e 1t . Hence MAR(1) already encompasses some of the nonlinearities in ecological dynamics. However, complex nonlinear dynamics (sensu [START_REF] May | Biological populations with nonoverlapping generations: stable points, stable cycles, and chaos[END_REF][START_REF] Turchin | Complex Population Dynamics: A Theoretical/Empirical Synthesis[END_REF] may not be completely captured by linear models, even in the logarithmic scale (Stenseth et al. 1998a[START_REF] Stenseth | Modelling non-additive and nonlinear signals from climatic noise in ecological time series: Soay sheep as an example[END_REF]. Modeling accurately differential responses to environmental variations like the storage effect [START_REF] Chesson | The roles of harsh and fluctuating conditions in the dynamics of ecological communities[END_REF] or relative nonlinearity of functional forms warrants models with a finer temporal structure. Fortunately, threshold autoregressive (TAR) models (Stenseth et al. 1998b[START_REF] Stenseth | Modelling non-additive and nonlinear signals from climatic noise in ecological time series: Soay sheep as an example[END_REF] can account for different temporal regimes by introducing different phases separated by threshold values -they are piecewise log-linear. TAR models have been successfully applied to diverse ecological studies (Stenseth et al. 1998a[START_REF] Stenseth | Modelling non-additive and nonlinear signals from climatic noise in ecological time series: Soay sheep as an example[END_REF]. They can be described by Eq. 5, where in our case two phases are defined:

n t+1 = a + B (1) n t + U (1) u t + (1) t if x ≤ θ B (2) n t + U (2) u t + (2) t if x > θ (5)
where n t is the log abundance at time t. The parameters B (1/2) , U (1/2) and

(1/2) t

have the same meaning as in the usual MAR model but can have different values in different phases. The phase is defined by x which can be both an intrinsic state variable (population abundance or growth rate in Stenseth et al. 1998a;[START_REF] Barraquand | Overcompensation and phase effects in a cyclic common vole population: between first and second-order cycles[END_REF] or an extrinsic variable [START_REF] Stenseth | Testing for effects of climate change on competitive relationships and coexistence between two bird species[END_REF] which is compared to a threshold θ. We describe below and in Supplementary material Appendix 5 how the definition of phases and thresholds map to ecological mechanisms like the storage effect and relative nonlinearity of competition.

In order to keep the analyses tractable and the number of parameters reasonable, given the length of the time series, we focused on the two main genera in Arcachon Bay: AST (Asterionellopsis) and CHA (Chaetoceros, see Table 1 for details). According to our findings using log-linear MAR models (see above), they are not supposed to interact, being pennate and centric diatoms respectively. We tested their interactions using the phase-based models to make sure that the absence of reported interaction was not due to hidden nonlinearities.

The first phase-based models, with phase defined through population densities (x is defined by AST and CHA abundance in Eq. 5), model possibly complex functional forms of competition between AST and CHA. Competition and response to environmental variables have then separate coefficients for high and low densities of the two plankton groups. This corresponds well to a "relative nonlinearity of competition" (see e.g., [START_REF] Fox | The intermediate disturbance hypothesis should be abandoned[END_REF]. The second, environment-based phase models (x and θ are defined by clusters of environmental conditions in Eq. 5) allow for coefficients to differ depending on environmental conditions. The latter models are particularly relevant to see if a so-called "storage effect" could happen, in which the effect of abiotic variables on population growth covaries with the values of the environmental variables themselves [START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF][START_REF] Ellner | How to quantify the temporal storage effect using simulations instead of math[END_REF].

Results

Environmental drivers of population dynamics

Here, we describe the correlation patterns of environmental drivers and phytoplankton densities in the frequency domain, that helped us, together with ecological knowledge, to select the key environmental drivers used later in ARX(p) and MAR(1) models (Fig. 2).

Spectral correlation analyses

Salinity was related to both freshwater inflow and precipitation, and having all three variables in our final models would not be desirable given their high correlation. As salinity cycles were more often correlated with plankton dynamics (Supplementary material Appendix 1: Fig. A1.18), and this variable also makes biological sense for the plankton niche (see Discussion), we retained only salinity. Irradiance and temperature were also correlated, which was of course logical given that a higher irradiance locally increases the water temperature. Both could have been considered, but we chose to use irradiance as the leading variable summarizing solar energy. Furthermore, on purely statistical grounds, CRY, the only group with differential responses to irradiance and temperature, were more sensitive to irradiance cycles (Supplementary material Appendix 1: Fig. A1.18(c)). Salinity, irradiance and wind energy were therefore the drivers in our physics-based models. We note that nutrients, especially N and Si, were also dependent upon freshwater inflow, which explained their high correlation (Supplementary material Appendix 1: Fig. A1.7).

Most plankton genera showed a strong seasonal response, except for CRY, GUI at Buoy 7 and SKE. The CRY group differed from other algae in its phenology by exhibiting much weaker (if any) periodicity. GUI and SKE may have a stronger long-term (supra-annual) component than other planktonic species (Supplementary material 6: Fig. A6.2(c) and Fig. A6.2(j)). Finally, shorter cycles (6 months) were apparent for AST, GUI and NIT; these cycles indicate two blooms per year.

Nutrient limitations

In a coastal lagoon with continental inputs such as AB, the odds that concentrations in nutrients would be high enough for plankton uptake to saturate are high. However, using half-saturation values from the literature, we could not demonstrate that nutrient concentrations were saturating for all genera at both sites (Table 2). Si seemed to be non-limiting (despite its general ecological importance for diatoms), as indicated by a saturating uptake function. The high correlation between Si and N led us therefore to select N rather than Si as the main nutrient driver, since it was more likely to be limiting. P seemed to be lacking in comparison to other nutrients, and was therefore considered potentially limiting as well.

Single-species autoregressive models

When taking all abiotic variables into account, the best model -in the sense of lowest AICcfor most planktonic groups used raw values of nutrient concentrations, as opposed to saturating functions of concentrations or nutrient ratios. This "absolute concentrations" model was consistent across genera: 7 genera out of 12 were within ±2 of the lowest AICc value across all considered models. The best model also included an explicit seasonal component (Table 3).

To decipher which of the abovementioned factors were most important to plankton dynamics, models with subsets of the environmental data were considered. Three models were evaluated (see Methods for details): the full model (physical drivers + nutrients), a model with physical variables only, and a model with nutrients only (Fig. 3). The physics-only model, which included irradiance, salinity, and wind energy, had the lowest AICc, followed by the full model and finally the nutrient-only model. For the two abundant and widely blooming genera (AST, CHA), the full and physics-only models were nearly equivalent in terms of ∆ AICc but clearly better than the nutrient-only model.

Even though the physics-only model did not account for all the variability in population growth rates (R 2 ranged between 0.12 and 0.27) and density-dependence was, in general, more influential to explain the overall dynamics than environmental factors, some environmental effects such as those of salinity and wind appeared consistent and common to most planktonic groups (Fig. 4). Seasonality had a strong impact on growth rates -the sign of which depends on the life history of the different genera. Residuals from integrated irradiance had a positive effect, except for AST which is known to be less sensitive to irradiance at Arcachon's latitude (early bloomer), and for PSE and RHI, whose growth rates were driven mostly by variations in wind energy and salinity. Except for a negligible effect on CRY and NIT, wind energy had strong negative impacts on population growth (Fig. 4).

The physics-only model was therefore used for MAR analyses in the following.

MAR(1) models

Model selection with AICc revealed that both the null model (no interaction) and the third scenario enabling interactions within centric diatoms, within pennate diatoms, and within dinoflagellates (eq. 4), were the best models to describe data at Teychan and Buoy 7 sites (Table 4).

The absence of interactions between pennate and centric diatoms was also confirmed by a more detailed analysis of interactions between the main pennate diatom AST and the main centric diatom CHA (see below). B and C matrices are presented for Teychan and Buoy 7 in Fig. 5, for the Physics-only ∆AICc cannot be seen for GUI, GYM and LEP as they are 0. ∆AICc for nutrients-only models were above 40 for CHA. Composition of planktonic groups on the x-axis is described in Table 1.

q -0.4 -0.2 0.0 0.2 0.4
Covariate effect q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q AST NIT PSE Figure 4: Mean values and standard errors of coefficients from linear models of plankton growth rates against irradiance, wind energy, salinity and a dedicated season component. Composition of planktonic groups on the x-axis is described in Table 1.

model differentiating pennate and centric diatoms (eq 4). Results based on other interaction matri- ) and comparisons between real data and simulations from the fitted model (Fig. 6, see also Supplementary material Appendix 2 and 5).

We found generally that intragroup interactions were at least two times higher than intergroup interactions (often four times higher or more). The interaction matrix structure was found to be similar between sites. Among diatoms, the strongest intragroup regulation was observed for CHA, the second most abundant genus in AB. AST, the most abundant genus, had an average intragroup interaction coefficient when compared to other phytoplanktonic groups. Intergroup coefficients were mostly low and positive when significantly different from zero (Fig. 5). 11 intergroup coefficients were positive at both Teychan and Buoy 7 sites, while only 5 of them were negative and consistent between sites. However, low positive and negative coefficient values might be artefacts, as the model with diagonal B was equally supported (Table 4).

Covariates played similar roles here as in single-genus models: salinity and wind energy had mostly negative effects while integrated irradiance had positive effects for the majority of planktonic genera. Our simulations of three contrasted ecological settings (Supplementary material Appendix 2: Section A2.2.1) followed by MAR model estimation revealed that the sign and rank of covariate effects were correctly estimated in all ecological settings (there could be a slight underestimation bias, c. 10%). and season) at Teychan (black) and Buoy 7 (blue), using an interaction matrix with no interaction between pennate (AST, NIT, PSE) and centric (SKE, GUI, CHA, GUI, LEP, RHI) diatoms. 95% significance of coefficients was determined by bootstrapping and is marked by asterisks (*). The figure should be read as element i having effect e ji on plankton group j. The identity matrix was subtracted to the interaction matrix (B-I) in order to make effects on growth rates comparable. The scale of coefficients is indicated by the blue bar of unit 1 in the upper right corner of the B-I matrix. Composition of planktonic groups is described in Table 1.

AST q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (a)

-3

-1 1 3 NIT q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (b) PSE q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (c)

SKE

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (d)

-3 -1 1 3
Scaled log abundance q Observation Estimation CHA q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (e) GUI q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (f) LEP q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (g) -3

-1 1 3 RHI q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (h) GYM q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (i) PRP q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (j)

2000 2005 2010 2015

-3 -1 1 3 CRY q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (k) 2000 2005 2010 2015 EUG q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (l) 2000 2005 2010 2015

Figure 6: Smoothed states obtained with the Kalman filter used in MARSS package (black line) compared to observations (red dots) for Teychan site, using an interaction matrix with no interaction between pennate (AST, NIT, PSE) and centric (SKE, GUI, CHA, GUI, LEP, RHI) diatoms. Composition of planktonic groups is described in Table 1.

Bivariate phase-dependent (TAR) models

The phases were defined either through the species densities (relative nonlinearity scenario) or the general environmental conditions for Teychan (storage effect scenario, Table 5). Phase definition is detailed in Supplementary material Appendix 5, and similar results for Buoy 7 can be found in always bloomed first and necessitated lower irradiance, being adapted to cold environments [START_REF] Kaczmarska | Cryptic diversity in a cosmopolitan diatom known as Asterionellopsis glacialis (Fragilariaceae): Implications for ecology, biogeography, and taxonomy[END_REF]. The effects of nutrients were not selected in autoregressive models, even though we considered absolute concentrations, saturating functions of these concentrations, and ratios of nutrient concentrations. There was, however, a possible impact of nitrogen but it was neither high nor consistent enough across groups to be considered a key driver, compared to other abiotic factors like irradiance or wind. Our results were similar for the two study sites, Teychan and Buoy 7, which were both largely influenced by the ocean, being respectively in the main channel and outer part of Arcachon Bay. Tidal movements may also connect Buoy 7 and Teychan populations; this could explain some of the similarity in the responses to environmental variables for the two stations.

Our models relied on the assumption of log-linearity in population growth rates (PGRs), which amounts to assume a multiplicative, power-law equation for the untransformed densities. Both PGR-log(density) curves (Supplementary material Appendix 2: Fig. A2.1) and examination of residuals (Supplementary material Appendix 2: Section A2.2.5) did not reject this hypothesis. Because examination of residuals of time series models can sometimes be ambiguous with noisy data, to be sure that alternative modeling choices could not alter our results, we fitted supplementary nonlinear competition models on the two main genera corresponding to pennate and centric diatoms (AST and CHA respectively; see Methods and Supplementary material Appendix 5). These included Ricker-based Lotka-Volterra (LV) competition (Ives et al. 1999b) (Supplementary material Appendix 5) and two phase-dependent models (i.e., piecewise log-linear Stenseth et al. 1998a[START_REF] Stenseth | Testing for effects of climate change on competitive relationships and coexistence between two bird species[END_REF]. The Ricker-LV model fitted the data to some extent but did not produce realistic dynamics (Supplementary material Appendix 5: Fig. A5.4) compared to MAR(1) models. The phase-dependent models mimicked relative nonlinearity of competition, by allowing the phase to depend on density values, and the storage effect mechanism [START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF], by allowing the phase to depend on environmental conditions (Supplementary material Appendix 5). The phase-dependent models produced very similar results to the basic MAR(1) log-linear model: no significant interactions were found in any phase, and the phases did not differ, which rejects the storage effect and relative nonlinearity. We are therefore confident that (1) the log-linear assumption at the heart of the MAR(1) model is appropriate for our data and (2) there is an absence of interactions between the two most abundant genera (AST and CHA), so that our MAR(1) model including phylogenetic information in the interaction matrix (i.e., weak interactions within pennate and centric diatom groups, as well as dinoflagellates) is quite robust.

We now discuss the likely reasons behind this absence of intergroup competition, as well as the effects of physical environmental drivers on phytoplankton population dynamics that we found.

Phytoplankton responses to environmental variation

There is a strong influence of the season on phytoplankton growth processes. We have taken into account this variation through an explicit seasonality variable to avoid collinear regressors in our models, as did most other MAR plankton studies (Ives et al. 1999a;[START_REF] Klug | Interactions among environmental drivers: Community responses to changing nutrients and dissolved organic carbon[END_REF][START_REF] Beisner | The effects of an exotic fish invasion on the prey communities of two lakes[END_REF]Hampton et al. 2006;Hampton and Schindler 2006;[START_REF] Hampton | Sixty years of environmental change in the world's largest freshwater lake -Lake Baikal, Siberia[END_REF][START_REF] Griffiths | Phytoplankton community interactions and environmental sensitivity in coastal and offshore habitats[END_REF][START_REF] Gsell | Quantifying change in pelagic plankton network stability and topology based on empirical long-term data[END_REF]. The season variable aggregates, of course, the effect of several covariates that have seasonal variation. Therefore, all the other variables, that are deseasonalized, have to be interpreted as deviations from the overall seasonal pattern.

With these limitations in mind, it is remarkable that both wind energy and salinity had mostly negative effects, meaning that more wind or more salinity than expected at a particular time of year has detrimental effects on planktonic population growth. Why wind has a consistently negative effect is not clear. Wind has different effects at different spatial scales. At very large spatial scales (>100 km 2 ), wind can have disruptive effects on oceanic currents and change the stratification outside coastal areas. In turn, a change in large-scale currents can create favorable conditions for blooms, as shown for the Bay of Biscay where AB is located [START_REF] Díaz | Climate variability and oceanographic settings associated with interannual variability in the initiation of Dinophysis acuminata blooms[END_REF]. At slightly smaller scales (m 2 to km 2 ), wind is generally expected to create large-scale turbulence, which is then cascaded to small-scale eddies and dissipated as heat at the micrometer scale. Despite the development of multiple modeling studies (Huisman et al. 1999a;[START_REF] Guasto | Fluid Mechanics of Planktonic Microorganisms[END_REF][START_REF] Nguyen | Hydrodynamics of diatom chains and semiflexible fibres[END_REF], field observations and experiments in the laboratory are not conclusive on the effect of small-scale turbulence (providing mixing and heat) generated by larger-scale phenomena [START_REF] Peters | Effects of turbulence on plankton: an overview of experimental evidence and some theoretical considerations[END_REF]. Microscale turbulence can lower the sinking rate of diatoms thanks to increased aggregation [START_REF] Reynolds | The ecology of phytoplankton[END_REF]) and improve nutrient assimilation [START_REF] Margalef | Life-forms of phytoplankton as survival alternatives in an unstable environment[END_REF][START_REF] Thornton | Diatom aggregation in the sea: mechanisms and ecological implications[END_REF]. However, other studies have shown that an increase in flow velocity results in decreased growth rates and increased mortality of phytoplankton [START_REF] Li | Effect of flow velocity on phytoplankton biomass and composition in a freshwater lake[END_REF][START_REF] Garrison | Effects of episodic turbulence on diatom mortality and physiology, with a protocol for the use of Evans Blue stain for live-dead determinations[END_REF], especially for dinoflagellates [START_REF] Peters | Effects of turbulence on plankton: an overview of experimental evidence and some theoretical considerations[END_REF][START_REF] Llaveria | Cell cycle and cell mortality of Alexandrium minutum (Dinophyceae) under small-scale turbulence conditions[END_REF]. Increased turbulence has also been shown to increase the encounter rate of prey and predators [START_REF] Kiørboe | Planktivorous feeding in calm and turbulent environments, with emphasis on copepods[END_REF], although the relationship between clearance rates and turbulence is nonmonotonic [START_REF] Saiz | Predatory and suspension feeding of the copepod Acartia tonsa in turbulent environments[END_REF]. Moreover, an increase in turbulence may lead to resuspension of sediments and therefore more turbid waters and lower light availability. The buffering of incoming solar energy could then explain the observed negative effect of wind energy on phytoplankton growth [START_REF] Gervais | Influence of small-scale turbulence and large-scale mixing on phytoplankton primary production[END_REF][START_REF] Glé | Typology of environmental conditions at the onset of winter phytoplankton blooms in a shallow macrotidal coastal ecosystem, Arcachon Bay (France)[END_REF].

Salinity is another dominant factor in other coastal areas, that can prove at times even more influential than nutrient loads [START_REF] Irwin | Phytoplankton niches estimated from field data[END_REF][START_REF] Gasiūnaitė | Seasonality of coastal phytoplankton in the Baltic Sea: Influence of salinity and eutrophication. -Estuarine[END_REF]. [START_REF] Scheef | Inferring plankton community structure from marine and freshwater long-term data using multivariate autoregressive models[END_REF] called for the use of salinity as a discriminating factor for estuarine environments in MAR analyses. Some diatoms have indeed lower growth rates when salinity increases [START_REF] Balzano | Effects of salinity on the growth rate and morphology of ten Skeletonema strains[END_REF], which could explain our mostly negative effects on plankton growth. An alternative hypothesis is that salinity is inversely related to freshwater inflow and therefore nutrient loads. As we found few direct effects of nutrients, it seems more logical to assume that salinity might have a negative effect for other reasons than its inverse relationship to nutrient load. Salinity is therefore most likely to influence phytoplankton growth directly through physiological adaptations (e.g., tolerance range due to osmotic processes, [START_REF] Kirst | Salinity tolerance of eukaryotic marine algae[END_REF] or indirectly through concentration effects other than nutrients (e.g., a higher freshwater inflow can dilute chemicals harmful to algae).

No (or weak) interspecific competition The absence of competition in our system was surprising, because we expected to find results similar to Descamps-Julien and Gonzalez (2005), i.e., a coexistence in spite of marked competition, mediated by differential responses to environmental variation. Our detailed comparison of nonlinear community models for AST versus CHA (Table 5, Supplementary material Appendix 5) is similar to their comparison of Fragilaria (pennate diatom) vs Cyclotella (centric diatom); except ours was field-based rather than based on experimental data. But we did not find such competition.

That said, nutrient depletion could occur on a fast timescale and stay undetected at the two weeks timescale. However, this fast resource sharing would likely imply nearly equal intra-and inter-group competition [START_REF] Scranton | Coexistence and emergent neutrality generate synchrony among competitors in fluctuating environments[END_REF], which we did not observe. Our observational results therefore contrast widely with the classic experiments of [START_REF] Titman | Ecological competition between algae: experimental confirmation of resource-based competition theory[END_REF] and [START_REF] Tilman | Phytoplankton community ecology: the role of limiting nutrients[END_REF] that showed competition for resources, in a context of well-defined niches related to nutrients. It is possible that under the differing circumstances considered here, entirely new drivers, e.g. predators or parasites, take over the regulation of phytoplankton growth. We discuss this possibility further below.

Among the few intergroup interactions that were consistent across sites, we observed more positive than negative interactions (Fig. 5). This can be due to hidden players that generate positive interactions (see below). Although we concede that there is a possibility for positive interactions to emerge if the seasonal variable does not fully correct for shared seasonal trends, positive interactions can also be genuine ecological phenomena, due for instance to predation (de Ruiter and Gaedke 2017). A number of plankton MAR studies have explicitely forbidden positive interactions between phytoplankton groups from their models (e.g., [START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF][START_REF] Huber | The role of predation for seasonal variability patterns among phytoplankton and ciliates[END_REF][START_REF] Hampton | Sixty years of environmental change in the world's largest freshwater lake -Lake Baikal, Siberia[END_REF], which in turn suggests that positive interactions could be more common than what previous multivariate time series analyses concluded.

What about neutral dynamics? Neutral dynamics sensu [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32)[END_REF], i.e., ecological drift in a zero-sum game, is likely not a main force here because of the fluctuations in numbers across 6 orders of magnitude, as we mentioned in our introduction. However, there are other models of neutral community dynamics that may apply here, such as [START_REF] Mutshinda | Integrating the niche and neutral perspectives on community structure and dynamics[END_REF] or [START_REF] Loreau | Species synchrony and its drivers: neutral and nonneutral community dynamics in fluctuating environments[END_REF]. [START_REF] Loreau | Species synchrony and its drivers: neutral and nonneutral community dynamics in fluctuating environments[END_REF] present a variant of the discrete-time Lotka-Volterra model (Ives et al. 1999b) with demographic and environmental stochasticity. In the neutral case, the species compete with each other by an equal amount (including intraspecifically), and yet vary synchronously under the influence of a joint environment. Stabilizing niche differences are progressively introduced through a parameter α that weights the importance of inter-relative to intra-specific/intra-group interactions. [START_REF] Mutshinda | Integrating the niche and neutral perspectives on community structure and dynamics[END_REF] do the same with Gompertz density-dependence. Our results may be relevant in two ways for theoretical community dynamics studies investigating similarly niche and neutral processes.

First, in our study, species not only vary in partial synchrony under a joint environment, but intragroup interactions also dwarf intergroup interactions (i.e., α small in the second model of [START_REF] Loreau | Species synchrony and its drivers: neutral and nonneutral community dynamics in fluctuating environments[END_REF]. In this case, the Loreau and de Mazancourt (2008) model is not neutral anymore and it contains stabilizing mechanisms for coexistence (sensu [START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF]. In our case, the reasons for such stabilizing niche differences are yet to be discovered. It would be of general interest to know whether other studies using long time series of phytoplankton can find results close to ours, or instead manage to get near-neutral dynamics. If coexistence occurs through stabilizing niche differences, observing intergroup interactions that are very small may be in fact logical: the criteria for coexistence in theoretical Lotka-Volterra communities require intergroup competition to be much lower than intragroup competition in large communities [START_REF] Barabás | The effect of intra-and interspecific competition on coexistence in multispecies communities[END_REF]. [START_REF] Mutshinda | Ecological equivalence of species within phytoplankton functional groups[END_REF] recently reported on ecological equivalence -neutrality -between species (or morphotypes) using a long-term phytoplankton dataset similar to ours (the L4 dataset in the English Channel). Ecological drift was modeled by assuming that the fraction of the phylum biomass attributed to each species was on average equal to the fraction realized in the previous week. Though we commend the extensive hierarchical modeling of [START_REF] Mutshinda | Ecological equivalence of species within phytoplankton functional groups[END_REF], the fraction realized is an aggregate community property, not directly related to neutral dynamics, unlike weights based on MAR coefficients developed by [START_REF] Mutshinda | Integrating the niche and neutral perspectives on community structure and dynamics[END_REF]. Better understanding the relationship between "neutrality" sensu [START_REF] Loreau | Species synchrony and its drivers: neutral and nonneutral community dynamics in fluctuating environments[END_REF]; Mutshinda andO'Hara (2011) vs. Mutshinda et al. (2016) might help to clarify whether the L4 phytoplankton community is truly neutral within phyla (diatoms vs dinoflagellates, see [START_REF] Mutshinda | Ecological equivalence of species within phytoplankton functional groups[END_REF]. If this was the case, the L4 phytoplankton community would indeed have widely different community dynamics than in AB where stabilizing mechanisms dominate.

Our second insight is that we found a clear log-linearity of the relationship between population growth rates and population densities (Supplementary material Appendix 2: Fig. A2.1, Appendix 2: Section A2.2.5, and Appendix 5), i.e.

r i,t = ln(N i,t+1 ) -ln(N i,t ) ∝ S j=1 b ij ln(N j,t ) (6)
This is the MAR formulation with Gompertz density-dependence, rather than Ricker densitydependence (used for discrete-time equivalents to the Lotka-Volterra model) as in Loreau and de Mazancourt ( 2008)

r i,t = ln(N i,t+1 ) -ln(N i,t ) ∝ S j=1 b ij N j,t (7)
where b ij = -r m α/K when i = j using Loreau and de Mazancourt ( 2008)'s notations. Gompertz community dynamics (Ives et al. 1999b;[START_REF] Mutshinda | Integrating the niche and neutral perspectives on community structure and dynamics[END_REF] might therefore provide a better niche or neutral null model than Ricker dynamics, despite the clear Lotka-Volterra connection of the latter [START_REF] Loreau | Species synchrony and its drivers: neutral and nonneutral community dynamics in fluctuating environments[END_REF]. Although some modifications of the Gompertz model might be useful around zero densities in a theoretical setting, it is clearly the better model for our data. An approximately Gompertz density-dependence might arise when large fluctuations in densities are combined with saturating functional relationships between process rates and densities, that weaken the effect of an added competitor or facilitator once high densities are reached.

Comparison to biotic drivers in similar ecosystems, and taxonomic issues We found no interaction between diatoms and our mixotrophic/heterotrophic dinoflagellates, which contrasts with the findings of [START_REF] Gsell | Quantifying change in pelagic plankton network stability and topology based on empirical long-term data[END_REF]. Such interactions between classes are highly variable in other studies: [START_REF] Scheef | Assessing marine plankton community structure from long-term monitoring data with multivariate autoregressive (MAR) models: a comparison of fixed station versus spatially distributed sampling data: Point versus spatial data MAR models[END_REF] only found a positive effect of dinoflagellates on diatoms while [START_REF] Griffiths | Phytoplankton community interactions and environmental sensitivity in coastal and offshore habitats[END_REF] estimated a negative effect of autotrophic and mixotrophic dinoflagellates on diatoms in their coastal site and a positive effect of diatoms on heterotrophic dinoflagellates in their offshore site (c. 40km away). Surprisingly, in our study, cryptophytes do not affect and are not affected by other phytoplankton groups at both sites; we expected this abundant group with known interactions to other groups to exhibit some significant interactions. In another plankton study, [START_REF] Klug | Compensatory dynamics in planktonic community responses to pH perturbations[END_REF] found a negative effect of cryptophytes on dinoflagellates in one location and no effect in another site [START_REF] Klug | Interactions among environmental drivers: Community responses to changing nutrients and dissolved organic carbon[END_REF]. Basic ecological information suggests that dinoflagellates can have a negative effects on cryptophyte population dynamics (through predation, [START_REF] Moeller | Acquired phototrophy stabilises coexistence and shapes intrinsic dynamics of an intraguild predator and its prey[END_REF], but these seem to be quite difficult to recover.

The above summary of plankton MAR studies suggests that patterns of interactions between broad plankton groups, at the class or phylum level (diatoms / dinoflagellates / cryptophytes etc.) are highly inconsistent between places and periods of study. These idiosyncratic results might be due -among other factors -to the aggregation of heterogenous genera and species with different traits and dynamics. For this reason, [START_REF] Griffiths | Phytoplankton community interactions and environmental sensitivity in coastal and offshore habitats[END_REF] suggested that improving taxonomic resolution, as we did, would help to estimate more precisely competitive or predatory interactions between plankton groups. However, we did not estimate stronger interactions despite having quite long time series, which suggests that interspecific competition at least may in fact not exist in many (phyto)plankton datasets. To our knowledge, our inference of interactions using MAR models on field data, at a fine taxonomic scale, can only be compared to the detailed freshwater study of [START_REF] Huber | The role of predation for seasonal variability patterns among phytoplankton and ciliates[END_REF]. Interspecific interactions within diatoms were only found between one centric species and one pennate species in [START_REF] Huber | The role of predation for seasonal variability patterns among phytoplankton and ciliates[END_REF], but were not a common pattern: their results tend therefore to concur with ours.

New techniques to structure interaction matrices applied to metabarcoding data may hold promise to infer interactions at the species level [START_REF] Ovaskainen | How are species interactions structured in species-rich communities? A new method for analysing time-series data[END_REF] in the future, although metabarcoding currently presents its unique set of challenges, such as data compositionality [START_REF] Cao | Inferring human microbial dynamics from temporal metagenomics data: Pitfalls and lessons[END_REF]. In any case, we recommend that future studies use simulations mimicking the study design as we (Supplementary material Appendix 2: Section A2.2.1) and [START_REF] Ovaskainen | How are species interactions structured in species-rich communities? A new method for analysing time-series data[END_REF] did. Simulations are essential to verify that statistical models are able to infer interactions for a given observational study design and a range of ecological settings.

Hidden (biotic) drivers Effects of phytoplankton natural enemies can be quite important [START_REF] Klug | Compensatory dynamics in planktonic community responses to pH perturbations[END_REF][START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF][START_REF] Huber | The role of predation for seasonal variability patterns among phytoplankton and ciliates[END_REF][START_REF] Gsell | Quantifying change in pelagic plankton network stability and topology based on empirical long-term data[END_REF]) but were unknown for AB.

There is only one study on the effect of viruses in AB, which shows no viral control of plankton but a high potential for viral lysis during blooms [START_REF] Ory | Pelagic food web patterns: do they modulate virus and nanoflagellate effects on picoplankton during the phytoplankton spring bloom?: Effects of viruses and nanoflagellates on picoplankton[END_REF]). Zooplankton is another potential compartment that could be involved in regulating phytoplankton, notably by creating the high intragroup density-dependence observed. Zooplankton populations show spatial heterogeneity in species composition [START_REF] Castel | Composition and differential distribution of zooplankton in Arcachon Bay[END_REF][START_REF] Tortajada | Network analysis of the planktonic food web during the spring bloom in a semi enclosed lagoon (Arcachon, SW France)[END_REF]) over AB, but most of the temporal biomass variation is synchronous (D'Elbee and [START_REF] D'elbee | Evolution spatio-temporelle du zooplancton dans le bassin d'Arcachon: le rôle des variables de milieu[END_REF]. If zooplankton plays a dynamical role here, this may explain why there are so little differences in coefficient estimates between Buoy 7 and Teychan.

Predators can induce strong intragroup competition [START_REF] Haydon | Pivotal assumptions determining the relationship between stability and complexity: an analytical synthesis of the stability-complexity debate[END_REF][START_REF] Vallina | Maximal feeding with active prey-switching: A killthe-winner functional response and its effect on global diversity and biogeography[END_REF][START_REF] Barabás | Emergent neutrality or hidden niches?[END_REF], which helps coexistence, but they can also contribute to coexistence in other ways. They can promote apparent facilitation between species, which could explain the positive interactions found in our study [START_REF] Abrams | Apparent competition or apparent mutualism? Shared predation when populations cycle[END_REF][START_REF] Barraquand | Indirect effects of primary prey population dynamics on alternative prey[END_REF][START_REF] De Ruiter | Emergent facilitation promotes biological diversity in pelagic food webs[END_REF]. However, some caution is needed; in a well-mixed environment, sharing predators can have the same dynamical consequences as sharing resources, and the species that resist better predation can end up excluding the others. This depends largely on the shape of predator numerical responses [START_REF] Holt | Predation, apparent competition, and the structure of prey communities[END_REF][START_REF] Abrams | Apparent competition or apparent mutualism? Shared predation when populations cycle[END_REF], that are clearly unknown in our case. Thus natural enemies can maintain diversity, but their effects vary depending on the context. However, they seem like a likely cause for explaining the high intraspecific density-dependence observed.

Parasitism, predation and competition with other algae (from pico-or nano-phytoplankton to macroalgae) are potential regulating factors [START_REF] Huber | The role of predation for seasonal variability patterns among phytoplankton and ciliates[END_REF], but other variables may emerge as drivers of plankton dynamics and coexistence. These range from toxic compounds and allelopathy [START_REF] Felpeto | Allelopathy prevents competitive exclusion and promotes phytoplankton biodiversity[END_REF] to auxotrophy [START_REF] Gobler | Effect of B-vitamins (B1, B12) and inorganic nutrients on algal bloom dynamics in a coastal ecosystem[END_REF][START_REF] Tang | Most harmful algal bloom species are vitamin B1 and B12 auxotrophs[END_REF] to hydrodynamic features.

Towards a more mechanistic study of phytoplankton community dynamics In land plants, coexistence can be relatively high-dimensional and depend on a great variety of traits, which are themselves not directly related to abiotic resources such as nitrogen or phosphorus [START_REF] Kraft | Plant functional traits and the multidimensional nature of species coexistence[END_REF]. This is also quite likely to happen for phytoplankton. The phytoplankton groups that we studied, diatoms in particular, often live in colonies of differing shapes and sizes: these factors are likely to impact community dynamics and coexistence. The colonial nature and shapes of diatoms are known to interact with patterns of very-small scale turbulence [START_REF] Margalef | Life-forms of phytoplankton as survival alternatives in an unstable environment[END_REF][START_REF] Reynolds | The ecology of phytoplankton[END_REF], and microscale hydrodynamics does affect plankton growth and coexistence (Huisman et al. 1999b). Much in the same way that small-scale spatial structure can be pivotal to land plant coexistence (e.g., through vertical shading, [START_REF] Onoda | Trade-off between light interception efficiency and light use efficiency: implications for species coexistence in one-sided light competition[END_REF], the interplay between microscale hydrodynamics and phytoplankton structure (cell shape, cell size, coloniality) is likely to be important for phytoplankton coexistence. Current models have started to investigate the joint effects of cell size and turbulence on persistence [START_REF] Portalier | Size-related effects of physical factors on phytoplankton communities[END_REF]), but much remains to be done with respect to cell shape and coloniality (see Nguyen and Fauci 2014 for a hydrodynamics study).

This study highlighted that intragroup competition is much stronger than intergroup competition in our coastal phytoplankton community. This contrasts widely with other types of models invoking temporal fluctuations in the environment (Descamps-Julien and Gonzalez 2005; Li and Chesson 2016) that counteract effects of competition for resources. However, this rejection of "paradox of the plankton models" does not imply that all interactions between species are unimportant to species coexistence. Other features of phytoplankton life histories may need to be explored to better comprehend the full effect of interactions on plankton interlinked population dynamics, e.g., sexual reproduction in phytoplankton ensures a rather rich demography [START_REF] D'alelio | The time for sex: A biennial life cycle in a marine planktonic diatom[END_REF]. Having several life-stages increases the possibility for interactions such as predation to affect community dynamics without these interactions being detectable using aggregated biomass or counts at the species or genus level [START_REF] Oken | How detectable is predation in stage-structured populations? Insights from a simulation-testing analysis[END_REF] -one of the limitations of count-based studies like ours.

To sum up, both demography and microscale hydrodynamics, together with natural enemies, might be needed to understand coexistence in phytoplankton. Although the sessile nature of land plants and the stability of forests contrast with the extremely motile, dynamic nature of the marine environment, some suggestions for further work may be transferable from plant to plankton ecology. A likely driver of stabilizing niche differences in tropical forests are Janzen-Connell effects where greater intraspecific than interspecific competition occurs at the seed and seedling stage, through combined effects of spatial structure of the seed rain and natural enemies [START_REF] Bagchi | Pathogens and insect herbivores drive rainforest plant diversity and composition[END_REF][START_REF] Comita | Testing predictions of the Janzen-Connell hypothesis: a meta-analysis of experimental evidence for distance-and density-dependent seed and seedling survival[END_REF]. Although the specific mechanisms at work will undoubtedly be different in phytoplankton, the idea that microscale spatial structure, stage structure and natural enemies may all interact to create stabilizing niche differences holds promise. A more demographic, small-scale perspective may help us to move beyond simple "paradox of the plankton" models and better explain the puzzling phytoplanktonic diversity. also like to acknowledge the contribution of Loïc Rigouin, who kindly provided Fig. 1, and Météo France, who contributed meteorological data. Grégoire Certain, Florian Hartig and Jarad Mellard provided constructive comments on drafts. We also thank three reviewers whose comments helped to improve the manuscript, as well as the associate editor. This study was supported by the French ANR through LabEx COTE (ANR-10-LABX-45).

Figure 1 :

 1 Figure 1: General view of Arcachon Bay and location of sampling points (Teychan and Buoy 7) within the study area. Exact coordinates and further information about the sampling sites are given in the main text.

  ,9 b 11,10 b 11,11

Figure 3 :

 3 Figure 3: Difference between the AICc of a full model, and the AICc of physics-only (salinity, irradiance, wind) and nutrients-only (nitrogen and phosphorus) models. All models use raw nutrient values and an explicit seasonal component. When ∆AICc is negative, the model performs better than the full model to describe the dynamics of each phytoplanktonic group.Physics-only ∆AICc cannot be seen for GUI, GYM and LEP as they are 0. ∆AICc for nutrients-only models were above 40 for CHA. Composition of planktonic groups on the x-axis is described in Table1.

Figure 5 :

 5 Figure5: Model coefficients for both variates and covariates (salinity, irradiance, wind energy and season) at Teychan (black) and Buoy 7 (blue), using an interaction matrix with no interaction between pennate (AST, NIT, PSE) and centric (SKE, GUI, CHA, GUI, LEP, RHI) diatoms. 95% significance of coefficients was determined by bootstrapping and is marked by asterisks (*). The figure should be read as element i having effect e ji on plankton group j. The identity matrix was subtracted to the interaction matrix (B-I) in order to make effects on growth rates comparable. The scale of coefficients is indicated by the blue bar of unit 1 in the upper right corner of the B-I matrix. Composition of planktonic groups is described in Table1.

Table 2 :

 2 Proportion of observations (%) at Teychan and Buoy 7 for which nutrient concentrations were above twice the half-saturation parameters found in the literature for different planktonic phyla (see the text for values and references). Dino. = Dinoflagellates.

				Phosphorus				Nitrogen		Silicon
		Diatoms Dino. Cryptophytes Euglenophytes Diatoms Dino. Cryptophytes Euglenophytes Diatoms
	Teychan	0	0	0	0	60	17	61	50	96
	Buoy 7	0	0	0	0	37	0.6	38	26	57

Table 3 :

 3 AICc for models explaining plankton dynamics in Arcachon Bay using physical variables (irradiance, wind energy, salinity) and nutrients (nitrogen and phosphorus), with and without a dedicated season variable. Nitrogen values are either used directly, or clipped with a saturating function. Nutrients were also input as ratios of P/N and Si/N. AICc values within 2 units of the minimum AICc are shown in bold letters for each group. Composition of planktonic groups is described in Table1

		No season	Season	No season	Season	No season Season
		No saturation No saturation Saturation Saturation	Ratio	Ratio
	AST	1430	1431	1426	1427	1431	1432
	CHA	1708	1703	1717	1706	1720	1710
	CRY	1187	1185	1186	1185	1185	1184
	EUG	1005	1003	1002	1002	996	995
	GUI	1268	1272	1269	1273	1271	1272
	GYM	1269	1260	1269	1260	1266	1264
	LEP	1115	1109	1113	1111	1113	1112
	NIT	1564	1545	1557	1546	1561	1546
	PRP	1128	1129	1129	1130	1121	1121
	PSE	1571	1566	1569	1569	1569	1569
	RHI	894	893	894	894	893	893
	SKE	804	811	806	806	808	810

Table 4 :

 4 Comparison of model selection criteria for different interaction matrices at Teychan and Buoy 7 sites (see eq. (4) and Supplementary material Appendix 2: Tables A2.5, A2.6, A2.7 for matrix definition). Values in bold correspond to the more parsimonious models.

		Null Unconstrained Pennate vs. Centric Diat vs. Dino Inter-phylum
	Teychan AICc 11970	12039	11964	11971	12025
	Teychan BIC 12271	12882	12398	12533	12631
	Teychan R 2 pred	0.52	0.55	0.50	0.50	0.54
	Buoy 7 AICc 8924	8980	8889	8920	8972
	Buoy 7 BIC 9196	9736	9280	9422	9517
	ces can be found in Supplementary material Appendix 2: Section A2.2.4. Model fit was checked by
	examination of residuals (Supplementary material Appendix 2: section A2.2.5

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 1 0 0.5 b ij

Acknowledgements

This study could not have been performed without the dedicated, long-term data collection of the REPHY programme undertaken by Ifremer, for which we are deeply grateful. This work owes much to Nadine Neaud-Masson, Myriam Rumèbe, and Claire Méteigner who performed phytoplankton counts, as well as to Florence d'Amico for hydrological measurements. We would

Supplementary material Appendix 5: Section A5.1.3.

The intergroup interaction coefficients were at least two times lower than the intragroup densitydependence and statistically non-significant. This was true in both bloom (high densities, or favorable environment) and non-bloom (low densities, or unfavorable environment) conditions at both sites. We therefore conclude that there are very few or no interactions between AST and CHA, even when considering more complex nonlinear models.

Discussion

Using time series data from both biotic and abiotic components of a coastal ecosystem, we have applied univariate and multivariate autoregressive models to pinpoint the factors driving the joint dynamics of a phytoplankton community. This phytoplankton community was dominated by several genera of pennate and centric diatoms as well as dinoflagellates (grouped at the genus level for most). Because most of these groups overlap in their resource requirements [START_REF] Reynolds | The ecology of phytoplankton[END_REF], we expected some degree of competition but relatively different, potentially nonlinear responses to environmental variables. This would have suggested a mechanism of coexistence through temporal variation of the environment [START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF], despite competition between genera. Our results painted in contrast a rather different picture, with high intragroup density-dependence but no or very weak competition between groups (genera or groups of genera). Although one of our best models included interactions between genera (within pennate diatoms, within centric diatoms and within dinoflagellates), the interaction coefficients that were significantly different from zero or consistent for both sites were usually positive. These emergent positive interactions may occur because the resulting net effects of all possible biotic interactions with unmeasured taxa are not necessarily negative (e.g., [START_REF] Stone | Conditions for a Species to Gain Advantage from the Presence of Competitors[END_REF], including those mediated by predators or parasites.

The responses to environmental variables were relatively similar among taxa, highlighting irradiance, wind energy and salinity, but there were some notable phenological differences. For example, the most abundant group, AST (mostly Asterionellopsis glacialis, see Table 1) almost