Loïck Kléparski
email: loick.kleparski@hotmail.fr

Grégory Beaugrand

Martin Edwards

François G Schmitt

Richard R Kirby

Elsa Breton

François Gevaert

Emeline Maniez

Morphological traits, niche-environment interaction and temporal changes in diatoms

Keywords: annual diatom succession, cell elongation, ecological niche, morphological traits, phenology, North Sea, CPR survey

come   L'archive ouverte pluridisciplinaire

Introduction

Diatoms are a major phytoplankton clade that have evolved during the Middle Triassic when sea levels were rising and continental margins were flooding [START_REF] Falkowski | The Evolution of Modern Eukaryotic Phytoplankton[END_REF]. Diatoms' subsequent increase in diversity during the mid-Cenozoic Era strongly influenced the global carbon cycle [START_REF] Katz | Biological overprint of the geological carbon cycle[END_REF]. Today, diatoms account for about 40% of total marine primary production and almost 40% of total particulate organic carbon exported to the deep ocean annually [START_REF] Jin | Diagnosing the contribution of phytoplankton functional groups to the production and export of particulate organic carbon, CaCO3 , and opal from global nutrient and alkalinity distributions[END_REF][START_REF] Nelson | Production and dissolution of biogenic silica in the ocean: Revised global estimates, comparison with regional data and relationship to biogenic sedimentation[END_REF][START_REF] Tréguer | Influence of diatom diversity on the ocean biological carbon pump[END_REF]. Annually, the abundance of diatoms rises and falls on a seasonal basis following a predictable pattern named "Annual Diatom Succession" (called hereafter ADS).

In the North Sea, ADS is characterised by the succession of various species throughout the year, with species such as Skeletonema costatum and Thalassiothrix longissima dominating in spring, Guinardia striata and G. flaccida in summer and Poroboscia indica and Bidulphia alternans in autumn [START_REF] Caracciolo | Annual phytoplankton succession results from niche-environment interaction[END_REF]. ADS therefore determines the pulses of energy that influence the dynamics of the whole marine ecosystem, the life cycle of many zooplankton and fish being coupled with the peaks in primary production or species dominance [START_REF] Cushing | Plankton Production and Year-class Strength in Fish Populations: an Update of the Match/Mismatch Hypothesis[END_REF][START_REF] Platt | Spring algal bloom and larval fish survival[END_REF]. Recently, anthropogenic climate change has begun to alter diatom phenology and biogeography, and this is likely to have strong consequences on trophic interaction and ecosystem functioning [START_REF] Chivers | Mismatch between marine plankton range movements and the velocity of climate change[END_REF][START_REF] Edwards | Impact of climate change on marine pelagic phenology and trophic mismatch[END_REF].

To predict the future consequences of climate change on marine ecosystems, we need to i) better understand how diatom communities naturally form and reorganise at multiple time-scales and ii) to identify the key elements and processes that control phenology. Morphological traits (i.e. cell size and shape) are known to play an important role in diatom succession because they influence resource uptake, buoyancy and predation, which are all properties that may explain seasonal changes in species dominance [START_REF] Behrenfeld | Abandoning Sverdrup's Critical Depth Hypothesis on phytoplankton blooms[END_REF][START_REF] Caracciolo | Annual phytoplankton succession results from niche-environment interaction[END_REF][START_REF] Karp-Boss | The Elongated, the Squat and the Spherical: Selective Pressures for Phytoplankton Shape[END_REF][START_REF] Margalef | Life-forms of phytoplankton as survival alternatives in an unstable environment[END_REF][START_REF] Naselli-Flores | Life in suspension and its impact on phytoplankton morphology: an homage to Colin S. Reynolds[END_REF][START_REF] Naselli-Flores | Invited Review -Fight on Plankton! or, Phytoplankton Shape and Size as Adaptive Tools to Get Ahead in the Struggle for Life[END_REF]. The niche of a species (i.e. the set of environmental parameters that enable a species to grow, maintain and reproduce) is also of paramount importance as it controls the range of environmental conditions that a species can tolerate [START_REF] Hutchinson | Concluding remarks[END_REF]. Here, we investigate the relationships between species' morphological traits (size and shape), ecological niches, and phenology with the aim of identifying the key elements and processes responsible for ADS.

Materials and Methods

Biological data

We used abundance data of 45 diatom species/taxa commonly sampled by the Continuous Plankton Recorder (CPR) survey in the North Sea (51°N, 60°N, -3°E, 9°E; Supplementary Figure S1) between 1958 and 2017 [START_REF] Beaugrand | Plankton effect on cod recruitment in the North Sea[END_REF]. This long-term plankton monitoring programme has collected plankton on a monthly basis in the extratropical regions of the North Atlantic Ocean and its adjacent seas since 1946. The sampling machine is a high-speed plankton recorder towed behind voluntary merchant ships at a depth of approximately 7 m (Beaugrand et al., 2003;[START_REF] Beaugrand | Plankton effect on cod recruitment in the North Sea[END_REF] (see Supplementary Text S1 for an overview of the CPR limitations and Supplementary Text S2 for CPR data availability) (Beaugrand et al., 2003;[START_REF] Beaugrand | Plankton effect on cod recruitment in the North Sea[END_REF]. Each value of abundance corresponds to a number of cells per CPR sample, which corresponds on average to 3 m 3 of seawater filtered [START_REF] Jonas | The volume of water filtered by a Continuous Plankton Recorder sample: the effect of ship speed[END_REF]. We calculated a daily climatology of the abundance of 45 diatoms sampled by the CPR survey. To account for all CPR samples collected between 1958 and 2017, we based our climatology on 366 days. We applied twice a 6-order simple moving average on the matrix (366 days x 45 species/taxa) and standardized the abundance data between 0 and 1 for each species, following the procedure described in Caracciolo and colleagues [START_REF] Caracciolo | Annual phytoplankton succession results from niche-environment interaction[END_REF].

Environmental data

Supplementary Text S2 provides key references and internet sites for each environmental dataset we used. Sea Surface Temperature (SST; °C) originated from the NOAA OI SST V2 high resolution dataset provided by the Earth System Research Laboratory, Physical Science Division, Boulder Colorado, USA.

Mixed Layer Depth (MLD; m) was provided by the Copernicus Marine Environment Monitoring service (CMEMS) and originated from the Global Ocean Physical Reanalysis product (GLOBAL_REANALYSIS_ PHY_001_030). We used the MLD as a proxy of turbulence level (the deeper the MLD is, the more active the wind-induced turbulence is likely to be) and stratification (a lower MLD is in general associated with a higher stratification, which prevents nutrients recycling to the surface).

Phosphate, silicate and nitrate concentrations (mmol.m -3 ) were provided by the CMEMS and came from the Global Ocean Biogeochemistry Hindcast (GLOBAL_REANALYSIS_BIO_001_029).

Photosynthetically Active Radiation (PAR; mole /m²/day) was provided by Globcolour Hermes. PAR was estimated from satellites data. Euphotic depth of sea water (m) was provided by the NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group. Euphotic depth was calculated with Lee algorithm. Bathymetry (m) was provided by the British Oceanographic Data Centre (BODC). Distance to nearest coastlines (km) was provided by the Nasa's Ocean Biology Processing Group, Goddard Space Flight Center.

We calculated a daily climatology for each environmental parameter in the North Sea (Figure 2b-h and Supplementary Figures S1 andS2a-c), based on the time period 1997-2017. The climatology of the euphotic depth of sea water was based on the time period 2002-2017. Here also, we based our climatologies on 366 days to match the diatom climatologies. We assumed that the difference of time periods between biological (1958-2017) and environmental (1997-2017) datasets did not alter our conclusions because the seasonal changes observed here, are much more important than decadal changes (see Figure 4). The use of the full time period for CPR data was important to have a more reliable daily climatology and we know that daily variability is much higher than decadal variability, an assumption that has been checked for copepod biodiversity assessed from the CPR data (Beaugrand et al., 2003).

The kinematic viscosity (in m 2 .s -1 ) was assessed from SST daily climatology using equations from Pilson [START_REF] Pilson | An Introduction to the Chemistry of the Sea[END_REF] and assuming a salinity of 35 (i.e. mean salinity of the North Sea).

Data on monthly long-term changes in SST (°C), mean surface downward short-wave radiation flux (W.m -2 ), and northward and eastward wind components (m.s -1 ) originated from the ERA5 dataset. This monthly gridded dataset, which covered the period 1958-2017, were provided by the Copernicus climate change service. Data on mean surface downward short-wave radiation flux (solar radiations) were used as a proxy of PAR.

Morphological traits

Cell shape

We gathered together information on the geometrical shape of each species from the Global Diatoms Database [START_REF] Leblanc | A global diatom database -abundance, biovolume and biomass in the world ocean[END_REF]. Six different shapes were distinguished: cylinder, cylinder + 2 half spheres, rectangular box, prism on an elliptic base, prism on a parallelogram base and prism on a triangular base (Supplementary Table S1).

Oblates and prolates

We also retrieved information on minimum and maximum cell dimensions (diameter/height or length/width/height depending upon cell shapes). For cylindrical cells, we defined height as the dimension perpendicular to cell diameter. For other geometrical shapes, we first converted them to a cylinder of equivalent size: we calculated the area of the slice perpendicular to cell height (as defined in the database) and estimated the diameter of the equivalent disk (i.e. the disk of equivalent surface). Based on these cell dimensions, we divided diatoms into three groups: (i) oblates (i.e. minimum height < minimum diameter and maximum height < maximum diameter), (ii) prolates (i.e. minimum height > minimum diameter and maximum height > maximum diameter) and (iii) oblates/prolates for the remaining species/taxa, meaning that they can be both oblates and prolates (Supplementary Table S2).

Other morphological indices

Based on mean cell height (ℎ) and diameter (𝑑) (Supplementary Table S2), we estimated the area (𝑠), the volume (𝑣) and the area/volume ratio (𝑠𝑣 -1 ) of each species by assimilating their shape to a cylinder:

𝑠 𝑐𝑦𝑙𝑖𝑛𝑑𝑒𝑟 = 𝜋𝑑 ( 𝑑 2 + ℎ) [1] 𝑣 𝑐𝑦𝑙𝑖𝑛𝑑𝑒𝑟 = 𝜋 4 𝑑 2 ℎ [2]
We assessed an index of cell morphology, which had the advantage to be dimensionless:

𝛿 = 𝑚𝑠𝑣 -1 [3]
with 𝑚 the cell maximum linear dimension (i.e. the mean cell diameter or height for oblates and prolates, respectively). 𝛿 describes cell attenuation (i.e. reduction of one of the dimensions of the cell volume or cell flattening) and its departure from a sphere [START_REF] Naselli-Flores | Life in suspension and its impact on phytoplankton morphology: an homage to Colin S. Reynolds[END_REF].

We estimated the maximum area of light interception normalized by volume (β). For oblates, we calculated the ratio between disk area and volume as follows:

𝛽 𝑜𝑏𝑙𝑎𝑡𝑒𝑠 = 𝜋( 𝑑 2 ) 2 𝑣 [4]
For prolates, we calculated the ratio between rectangle area (projection of a cylinder on a surface) and volume:

𝛽 𝑝𝑟𝑜𝑙𝑎𝑡𝑒𝑠 = 𝑑ℎ 𝑣 [5]

Predation

We used abundance data of 74 small (≤ 2 mm) and 196 large (> 2 mm) copepods species/taxa collected by the CPR survey in the North Sea (Supplementary Figure S1). We summed the total abundance for small/large copepods in each CPR sample between 1958-2017 and then calculated a daily abundance climatology for small/large copepods. We therefore obtained two vectors of 366 days (one vector for small and another for large copepods) with the daily mean abundance for small/large copepods. As above, we applied twice a simple 6-order moving average and standardized the abundance data between 0 and 1 (Supplementary Figure S2c).

Phylogenetical classification

Information on the phylogenetical classification of the 45 diatoms (family, order, superorder, subclasses) was retrieved from the WOrld Register of Marine Species (WORMS) (http://www.marinespecies.org/index.php). (Supplementary Figure S3 and Table S3).

Analyses

Principal Component Analysis (PCA)

We characterized Annual Diatom Succession (ADS) by applying a standardized PCA on a matrix of daily abundance climatology (366 days) x diatom species (45 species) (Figure 1a). We examined the first three axes (principal components and eigenvectors) of the PCA because a test based on a broken-stick distribution [START_REF] Legendre | Numerical Ecology[END_REF] revealed that the first three axes of the PCA were significant, i.e. they were above average eigenvalues calculated by the broken-stick model (Supplementary Table S4). Seasonal changes in diatom abundance were examined by the first three principal components (PCs) and relationships between species were investigated by using the first three normalized eigenvectors, which represented the correlation of the 45 diatoms with the first three PCs. We tested the relationships between the first three principal components and the daily climatology for SST, euphotic depth, PAR, MLD, mean small/large copepods abundance, viscosity, silicate, phosphate and nitrate concentration by means of a Spearman rank correlation coefficient.

The significance of the correlation coefficients was tested by applying a Montecarlo test with 100,000 simulations (Figure 3 and Supplementary Table S5).

Clustering

We applied a cluster analysis based on the first three normalized eigenvectors of the PCA. First, we calculated a squared matrix of Euclidean distance using the first three eigenvectors. Then we applied a Ward linkage algorithm [START_REF] Legendre | Numerical Ecology[END_REF]. The resulting dendrogram was cut using a threshold of Euclidean distance of 8; this threshold was chosen because it enabled the identification of three species groups that closely matched the seasonal patterns exhibited by the first three PCs (Figure 1a, 3 and Supplementary Figures S2d-f, S4 and S5).

Relationships between phenology and morphological traits

We examined the relationships between the species groups and morphological traits in the space defined by the first three eigenvectors: (i) geometrical shapes (Figure 1b), (ii) oblates/prolates (Figure 1c), (iii) height-diameter difference (Figure 1d), (iv) area/volume ratio (Figure 1e), (v) δ (Figure 1f), (vi) mean cell diameter (Supplementary Figure S6a) and (vii) mean cell height (Supplementary Figure S6b).

We tested the differences in the mean of the morphological indices (height-diameter (Supplementary Table S6), mean cell diameter (Supplementary Table S7), mean cell height (Supplementary Table S8), area/volume ratio (Supplementary Table S9) and δ (Supplementary Table S10)) among species groups by means of a Kruskal-Wallis test and post-hoc tests.

Comparison between morphological traits based on cylindrical and spheroidal shapes

In our analyses, we assimilated diatom shapes to cylinders. To test whether our conclusions were affected by this assumption, we also compared our results with those obtained by assimilating diatom shapes to spheroids. Based on mean cell height (ℎ) and diameter (𝑑) (Supplementary Table S2), we re-estimated the area (𝑠), the volume (𝑣) and the area/volume ratio (𝑠𝑣 -1 ) of each species by assimilating their shape to spheroids as follows:

𝑠 𝑠𝑝ℎ𝑒𝑟𝑜𝑖𝑑 = 𝜋𝑑 2 (𝑑 + ℎ 2 √ℎ 2 -𝑑 2 sin -1 √ℎ 2 -𝑑 2 𝑑 ) [6] 𝑣 𝑠𝑝ℎ𝑒𝑟𝑜𝑖𝑑 = 𝜋 6 𝑑 2 ℎ [7]
We examined the relationships between morphological traits (area, volume, area/volume ratio and δ) when cells were based on cylinders (see equations 1-3) and spheroids by means of Pearson correlation coefficients (Supplementary Figure S7a-d).

We also investigated the relationships between species groups and morphological traits estimated with spheroidal shape in the space defined by the first three eigenvectors using (i) area/volume ratio and (ii) δ (Supplementary Figure S7e-h).

Finally, we tested the differences in the mean of the morphological indices among species groups, based on spheroidal shapes for area/volume ratio (Supplementary Table S11) and δ (Supplementary Table S12), by means of a Kruskal-Wallis test and post-hoc tests.

Relationships between morphological indices and the environment

We also examined the relationships between 𝑚 and area/volume ratio (Figure 1g-h), PAR and β (Supplementary Figure S8), MLD and 𝛿 (Supplementary Figure S9), kinematic viscosity and 𝛿 (Supplementary Figure S10) and between copepod abundance and height-diameter difference (Supplementary Figure S11).

We calculated the Spearman rank correlation coefficient between 𝑚 and area/volume ratio (Figure 1g), PAR and β (Supplementary Figure S8), MLD and 𝛿 (Supplementary Figure S9), kinematic viscosity and 𝛿 (Supplementary Figure S10) and between copepod abundance and height-diameter difference (Supplementary Figure S11). The correlations were tested by means of a Montecarlo test using 100,000 permutations [START_REF] Jackson | Are probability estimates from the permutation model of Mantel's test stable?[END_REF].

Phenology and niche

We estimated the optimum and amplitude (breadth) of the niche using nine ecological dimensions:

SST, euphotic depth, PAR, MLD, silicate phosphate and nitrate concentration, bathymetry and distance to nearest coastlines. To do so, we estimated the weighted mean as a proxy for optimum (𝑋 𝑜𝑝𝑡 ) and the weighted standard deviation as a proxy for the amplitude (𝜎) (Figure 2 and Supplementary Figure S12). These proxies were calculated for each ecological dimension:

𝑋 𝑜𝑝𝑡 = ∑ 𝑤 𝑖 𝑥 𝑖 𝑡 𝑖=1 ∑ 𝑤 𝑖 𝑡 𝑖=1 [8] 𝜎 = ∑ 𝑤 𝑖 (𝑥 𝑖 -𝑥 𝑜𝑝𝑡 ) 2 𝑡 𝑖=1 𝑡-1 𝑡 ∑ 𝑤 𝑖 𝑡 𝑖=1 [9]
where 𝑤 is the abundance of the diatom species, 𝑥 the value of the environmental parameter or copepod abundance and 𝑡 the number of days in the climatology (for SST, euphotic depth, PAR, MLD, silicate phosphate and nitrate concentration, and small/large copepod abundance) or the total number of CPR samples when bathymetry or distance to nearest coastlines was used.

We examined the relationships between diatom niche optimum/amplitude, their day of maximum abundance and the daily climatology for SST, euphotic depth, PAR, MLD, silicate phosphate and nitrate concentration (Figure 2b-h). We also investigated the relationships between diatom niche optimum/amplitude and height-diameter difference (Supplementary Figure S12).

We tested the niche separation between oblates and prolates by means of a Krukal-Wallis test performed on niche optima and amplitude (Supplementary Tables S13-S14). The relationships between species niche optima/amplitudes and species height-diameter difference was examined by applying a Spearman rank correlation coefficient and tested with a Montecarlo test using 100,000 simulations (Supplementary Figure S12) [START_REF] Jackson | Are probability estimates from the permutation model of Mantel's test stable?[END_REF].

Long-term changes in the abundance of oblates and prolates

We investigated long-term monthly changes in the mean abundance of oblates and prolates collected by the CPR survey between 1958 and 2017 (Figure 4a). Oblate abundance in April and October and prolate abundance in August were smoothed by means of a second-order simple moving average (Figure 4b-d). The relationships between species abundance (in April and October for oblates and August for prolates) and the long-term changes in SST, wind and solar radiations were examined by applying a Spearman rank correlation coefficient and tested with a Montecarlo test using 100,000 simulations (Supplementary Table S15) [START_REF] Jackson | Are probability estimates from the permutation model of Mantel's test stable?[END_REF].

Results and Discussion

We first characterised ADS by means of a Principal Component Analysis (PCA) performed on a matrix of daily abundance climatology (366 days) x diatom species (45 species; Methods). The abundance matrix originated from data collected by the Continuous Plankton Recorder (CPR) survey between 1958 and 2017 in the North Sea (Supplementary Figure S1). We retained the first three Principal Components (PCs), which were significant after applying a broken-stick test (Methods and We next studied relationships between species phenology and morphological traits. Cell shape and size are major traits affecting nutrient acquisition, light capture, grazing resistance, and vertical motion [START_REF] Clifton | Enhanced sedimentation of elongated plankton in simple flows[END_REF][START_REF] Karp-Boss | The Elongated, the Squat and the Spherical: Selective Pressures for Phytoplankton Shape[END_REF][START_REF] Litchman | Trait-Based Community Ecology of Phytoplankton[END_REF][START_REF] Naselli-Flores | Life in suspension and its impact on phytoplankton morphology: an homage to Colin S. Reynolds[END_REF][START_REF] Naselli-Flores | Invited Review -Fight on Plankton! or, Phytoplankton Shape and Size as Adaptive Tools to Get Ahead in the Struggle for Life[END_REF][START_REF] Padisak | Sinking properties of some phytoplankton shapes and the relation of form resistance to morphological diversity of plankton -an experimental study[END_REF]. To do so, we compiled data on cell size and geometrical shape from the Global Diatom Database [START_REF] Leblanc | A global diatom database -abundance, biovolume and biomass in the world ocean[END_REF] for 45 species/taxa (Methods). Spring-dominated and autumn-dominated groups comprised a mixture of species of different shapes (cylindrical, rectangular box and prism) whereas the LSSA group mainly comprised cylindrical diatoms (Figure 1b and Supplementary Table S1). By assuming all diatom shapes to be cylindrical (note that we also performed the same analyses by assimilating diatom shapes to spheroids, see below), we separated diatoms as oblates (i.e., cell diameter > cell height, flattened cell) and prolates (cell diameter < cell height, elongated cell). We identified 20 oblates and 21 prolates, the remaining being identified as oblates/prolates (i.e. species/taxa that can be both during their life cycle; Methods and Supplementary Table S2). The spring-and autumn-dominated groups mainly comprised oblates (62.5% and 58% in the spring-and autumn-dominated group, respectively; Figure 1c and Supplementary Table S1) and the LSSA group mostly comprised prolates (82%; Figure 1c and Supplementary Table S1).

Following this result, we investigated the variations of the mean cell height-diameter difference between groups more deeply (Figure 1d). Highest and lowest positive differences characterized the LSSA group and the spring-dominated group, respectively. Values were more variable from the autumn-dominated group. A Kruskal-Wallis test revealed significant differences and post-hoc tests showed that differences occurred between LSSA and spring/autumn-dominated groups, but not between spring-and autumn-dominated groups (Supplementary Table S6). Interestingly, negative differences (shown in blue on the colour bar) were smaller than positive ones (from green to red;

Figure 1d), suggesting that cell height rather than diameter is modified during ADS; this result is confirmed by insignificant differences in cell diameter and significant differences in cell height among groups (Supplementary Figure S6 and Supplementary Tables S7-S8). Previous studies have suggested that the size of phytoplankton is rarely greater than 300 µm (Naselli-Flores and Barone, 2011) because of fundamental structural constraints imposed by natural selection, i.e. limited minimum thickness and maximum diameter. For example, to avoid turbulence-induced damage and to benefit from the viscosity of water, planktonic cells have to remain smaller than the Kolmogorov scale (i.e.

micro eddies) [START_REF] Naselli-Flores | Life in suspension and its impact on phytoplankton morphology: an homage to Colin S. Reynolds[END_REF][START_REF] Reynolds | Plants in motion: Physical -biological interaction in the plankton[END_REF] and the area/volume ratio needs to remain relatively low to ensure chemical exchanges between the cell and its environment [START_REF] Raven | Biology[END_REF].

Phytoplankton exhibit different cell shapes that modify their area/volume ratio to optimize nutrient exchange and buoyancy with their surrounding environment (Karp-Boss and Boss, 2016; [START_REF] Pahlow | Impact of cell shape and chain formation on nutrient acquisition by marine diatoms[END_REF]. Energy flux occurs through the cell surface while metabolism is related to the cell volume, which explains why the ratio between surface area and volume should be high, although the building of a large surface has also a negative metabolic cost [START_REF] Brown | Toward a metabolic theory of ecology[END_REF]. We did not find any significant seasonal changes in the ratio between surface area and volume amongst our three groups (Figure 1e and Supplementary Table S9). Similar results have been documented in tropical lakes and coastal marine waters where the area/volume ratio was preserved through morphological changes [START_REF] Lewis | Surface/Volume Ratio: Implications for Phytoplankton Morphology[END_REF][START_REF] Stanca | Geometric shape as a trait to study phytoplankton distributions in aquatic ecosystems[END_REF]. We tested whether elongation (i.e. an increase in the maximum linear dimension or 𝑚) could be a more relevant parameter. To do so, we used an index of cell morphology 𝛿 (i.e. 𝑚 times the area/volume ratio; Methods) [START_REF] Naselli-Flores | Life in suspension and its impact on phytoplankton morphology: an homage to Colin S. Reynolds[END_REF]. In contrast to the area/volume ratio (Figure 1e), we found significant differences in cell morphology (δ) among the three groups (Figure 1f and Supplementary Table S10). Assimilating diatom shapes to spheroids led to similar conclusions (Supplementary Figure S7 and Tables S11-S12). A closer look on 𝑚 versus 𝑠𝑣 -1 allowed a better understanding of why there was no significant change in the area/volume ratio (Figure 1g-h). The two figures suggest that the dominance of prolates in the LSSA group during the stratified low-nutrients period does not affect the area/volume ratio in contrast to oblates where changes in size do influence the area/volume ratio. These results suggest that cell elongation is a way of altering shape without changing the area/volume ratio that would otherwise affect nutrient exchange [START_REF] Karp-Boss | The Elongated, the Squat and the Spherical: Selective Pressures for Phytoplankton Shape[END_REF][START_REF] Pahlow | Impact of cell shape and chain formation on nutrient acquisition by marine diatoms[END_REF].

Further, we examined whether phylogeny (from family to subclass level) was responsible for the observed differences in traits. We found no evidence for a phylogenetic effect at the family, order and superorder levels. This was not so at the subclass level; although the LSSA and autumndominated groups were mainly composed of Coscinodiscophycidae, the spring-dominated group was characterized by a diversity of subclasses (Supplementary Figure S3 and Table S3). This analysis, therefore, shows that differences in morphological traits are also influenced by phylogeny. Other studies have suggested that major morphological traits that drive strategies in nutrient uptake are conserved over evolutionary time at a high phylogenetic level, a phenomenon known as trait conservatism [START_REF] Litchman | The role of functional traits and trade-offs in structuring phytoplankton communities: scaling from cellular to ecosystem level[END_REF]. Trait conservatism results from stabilizing selection and niche conservatism, i.e. species keeping the shape of their niche unaltered over time because of their propensity to remain in places where the environment is highly suitable for them [START_REF] Ackerly | Community Assembly, Niche Conservatism, and Adaptive Evolution in Changing Environments[END_REF].

(Not all phytoplankton traits exhibit the same level of conservatism, some being more labile than others. [START_REF] Litchman | Linking traits to species diversity and community structure in phytoplankton[END_REF][START_REF] Schwaderer | Eco-evolutionary differences in light utilization traits and distributions of freshwater phytoplankton[END_REF]).

Three other hypotheses could also explain the morphological changes we observed: (i) light capture optimization (Naselli-Flores and Barone, 2011), (ii) buoyancy [START_REF] Naselli-Flores | Life in suspension and its impact on phytoplankton morphology: an homage to Colin S. Reynolds[END_REF][START_REF] Naselli-Flores | Invited Review -Fight on Plankton! or, Phytoplankton Shape and Size as Adaptive Tools to Get Ahead in the Struggle for Life[END_REF][START_REF] Padisak | Sinking properties of some phytoplankton shapes and the relation of form resistance to morphological diversity of plankton -an experimental study[END_REF] and (iii) response to predation [START_REF] Smetacek | A watery arms race[END_REF]. We investigated relationships between maximum area of light interception normalized by volume (β) and

Photosynthetically Active Radiation (PAR) (Methods, Supplementary Figure S8). We found no clear relationship between PAR and β. Noticeably, β diminished when maximum PAR increased, which is in agreement with physiology [START_REF] Reynolds | Physical determinants of phytoplankton succession[END_REF]. However, small values of β were observed for a large range of PARs, which suggests that changes in β are not modulated by PAR for most species considered here. We acknowledge that physiological processes by which the cell can adapt to different light intensity are complex, however (e.g. photosynthetic pigments, see Supplementary Text S3) [START_REF] Reynolds | Physical determinants of phytoplankton succession[END_REF].

Many other studies have provided evidence that variations in cell shape might result from an adaptation to turbulence through a modification in sinking velocity [START_REF] Clifton | Enhanced sedimentation of elongated plankton in simple flows[END_REF][START_REF] Margalef | Life-forms of phytoplankton as survival alternatives in an unstable environment[END_REF][START_REF] Naselli-Flores | Life in suspension and its impact on phytoplankton morphology: an homage to Colin S. Reynolds[END_REF][START_REF] Naselli-Flores | Invited Review -Fight on Plankton! or, Phytoplankton Shape and Size as Adaptive Tools to Get Ahead in the Struggle for Life[END_REF][START_REF] Padisak | Sinking properties of some phytoplankton shapes and the relation of form resistance to morphological diversity of plankton -an experimental study[END_REF]. We tested this hypothesis by investigating the relationship between Mixed Layer Depth (MLD, used as a proxy for wind-induced turbulence intensity) and 𝛿 (Supplementary Figure S9, Methods). This analysis revealed a negative relationship between maximum MLD values and 𝛿. While oblates were found for a large range of MLDs, prolates dominated low-turbulence waters. We therefore concluded that diatom elongation observed during ADS may result mainly from an adaptation to buoyancy, although resource acquisition and predation probably also play a role [START_REF] Karp-Boss | The Elongated, the Squat and the Spherical: Selective Pressures for Phytoplankton Shape[END_REF][START_REF] Pahlow | Impact of cell shape and chain formation on nutrient acquisition by marine diatoms[END_REF][START_REF] Sommer | Silicate and the functional geometry of marine phytoplankton[END_REF].

Studies have also found negative relationships between temperature and phytoplankton sinking velocity through an alteration in viscosity [START_REF] Smayda | The suspension and sinking of phytoplankton in the sea[END_REF][START_REF] Sournia | Form and function in marine phytoplankton[END_REF]. From Stokes' equation, and considering a cell of 20 µm of diameter, sinking velocity increases of 4% for each increase of 1°C because of a decrease in viscosity [START_REF] Smayda | The suspension and sinking of phytoplankton in the sea[END_REF]. However, sinking velocity depends on cell shapes [START_REF] Padisak | Sinking properties of some phytoplankton shapes and the relation of form resistance to morphological diversity of plankton -an experimental study[END_REF][START_REF] Smayda | The suspension and sinking of phytoplankton in the sea[END_REF][START_REF] Sournia | Form and function in marine phytoplankton[END_REF]. We therefore tested the relationship between kinematic viscosity and δ and found a weak but significant negative relationship (Supplementary Figure S10; Methods). This result suggests that prolates may also occur in stratified warm waters because their shapes slow down sinking velocity [START_REF] Padisak | Sinking properties of some phytoplankton shapes and the relation of form resistance to morphological diversity of plankton -an experimental study[END_REF]. We acknowledge that we have not considered all processes that can impact phytoplankton buoyancy. First, many species produce organic substances that increase the viscosity of the surrounding environment and that may then trigger cell aggregation, sometimes associated with minerals and particles, which may ultimately increase sinking velocity [START_REF] Smayda | The suspension and sinking of phytoplankton in the sea[END_REF][START_REF] Smetacek | Role of sinking in diatom life-history cycles: ecological, evolutionary and geological significance[END_REF]. Second, phytoplankton cells can also develop protuberances such as spines and setae that directly increase the frictional drag or help turn the longest cell dimension perpendicular to the settling axis [START_REF] Padisak | Sinking properties of some phytoplankton shapes and the relation of form resistance to morphological diversity of plankton -an experimental study[END_REF][START_REF] Smayda | The suspension and sinking of phytoplankton in the sea[END_REF].

Third, some species form chains that may affect sinking speed, although contrasting results have been observed for different species with either an increase or a decrease in buoyancy [START_REF] Padisak | Sinking properties of some phytoplankton shapes and the relation of form resistance to morphological diversity of plankton -an experimental study[END_REF][START_REF] Smayda | The suspension and sinking of phytoplankton in the sea[END_REF]. Lastly, we have not considered intra-specific variations in size and seasonal morphological plasticity, which can be important in many species [START_REF] Ligowski | Morphological adaptation of a planktonic diatom to growth in Antarctic sea ice[END_REF].

Lastly, we investigated whether there is a relationship between diatom shape and predation. We examined the relationship between the abundance of small (≤2mm)/large(>2mm) copepods and morphological shape (Methods, Supplementary Figure S11). Although correlation cannot imply causation, we found a nonlinear positive correlation between zooplankton abundance for small and large copepods and height-diameter difference. This result is in agreement with works that have suggested that changes in diatom shape may be influenced by zooplankton predation [START_REF] Karp-Boss | The Elongated, the Squat and the Spherical: Selective Pressures for Phytoplankton Shape[END_REF][START_REF] Naselli-Flores | Invited Review -Fight on Plankton! or, Phytoplankton Shape and Size as Adaptive Tools to Get Ahead in the Struggle for Life[END_REF][START_REF] Smetacek | A watery arms race[END_REF], although other works have shown that elongated cells can be ingested by copepods when turbulence is low [START_REF] Visser | On the reorientation of non-spherical prey particles in a feeding current[END_REF]. We caution that defence mechanisms in diatoms are diverse, ranging from physiological (allelopathy) to phenomena such as colony formation and escape response [START_REF] Pančić | Phytoplankton defence mechanisms: traits and trade-offs: Defensive traits and trade-offs[END_REF] and that all mechanisms involved could not be realistically tested in this study.

Recent works have suggested that annual plankton succession results from the niche-environment interaction [START_REF] Caracciolo | Annual phytoplankton succession results from niche-environment interaction[END_REF]. We therefore explored the relationships between cell shape and ecological niche for the 45 species/taxa. Here, the ecological niche was characterised by its optimum and amplitude (i.e. niche breadth) [START_REF] Hutchinson | Concluding remarks[END_REF]; these niche parameters were assessed for each species/taxa and nine environmental parameters (Methods; Figure 2a). Our results suggest that morphological traits influence species' niche of diatoms and that oblates and prolates had significant differences in ecological optima for all niche dimensions. Amplitudes were also distinct for four niche dimensions: nutrients (nitrate, phosphate and silicate) and MLD (Figure 2b-h, Supplementary Tables S13-S14). Although the relationship was less obvious for sea surface temperature and euphotic depth (Figure 2b-c), optima and amplitudes plotted against daily climatology of each parameter showed a clear dominance of prolates during low-nutrients, turbulence and PAR conditions (Figure 2d-h). In addition, oblates were significantly more euryoecious (i.e. larger niche amplitude) than prolates with respect to nutrients and MLD (Figure 2e-h, Supplementary Figure S12 and Table S14). Prolates are therefore K-strategists (i.e. equilibrium species), adapted to stable environment (i.e. predictable environment with low perturbations) with respect to nutrients, turbulence and stratification (low MLD) whereas oblates are more r-strategists (i.e. opportunist species), adapted to a more unstable environment (i.e. unpredictable) [START_REF] Southwood | Ecological strategies and population parameters[END_REF]. In other words, prolates (K-strategist)

thrive in stratified water with low turbulence and constant low-nutrient concentration whereas oblates (r-strategists) thrive in turbulent waters in the absence of stratification with more variable levels of nutrients concentration (Figure 2). Because oblates and prolates have distinct niches, this analysis demonstrates that morphological traits influence species niche and phenology. We summarize the relationships between diatom morphological traits, phenology and the environment in Figure 3 (see also Supplementary Table S5). S2). Oblates are in blue, prolates in red and oblates/prolates in green. Species names are those recorded in the CPR dataset. Blue area delineates the space occupied by the springdominated group, green area the space occupied by the autumn-dominated group and red area the space occupied by the LSSA group. Correlation between the first three principal components and the environmental parameters are displayed in Supplementary Table S5. More information on the principal components can be found in Figure 1a.

Finally, we investigated long-term monthly changes in the abundance of oblates and prolates between 1958 and 2017. Our results show that both groups exhibit distinct patterns of temporal changes (Figure 4a), with oblates having a constant and high level of abundance in April (Figure 4b) but a more variable abundance in October with blooms of decreasing intensity during the 1960s and 1970s, followed by a strong increase after circa 1980 (Figure 4c). Similar patterns are also observed for prolates in August but with an abundance level below the maximum reached in 1960 (Figure 4d).

The period of low abundance found for oblates and prolates between 1978 and 1982 was probably the consequence of a cold hydro-meteorological anomaly observed in the North Sea [START_REF] Edwards | Ocean climate anomalies and the ecology of the North Sea[END_REF]. We subsequently investigated whether these long-term changes in abundance were correlated with long-term changes in Sea Surface Temperature (SST), wind and solar radiations (Supplementary Table S15). We found a nonlinear positive correlation between oblate abundance in October and SST and a negative correlation with eastward wind as well as wind intensity. A negative correlation was also found between oblate abundance in April and eastward wind. No significant correlations with northward wind and solar radiations were found. In the Northeast Atlantic region, long-term changes in SST/wind variations have led to a marked increase in diatoms, associated with a diminution in dinoflagellates [START_REF] Hinder | Changes in marine dinoflagellate and diatom abundance under climate change[END_REF]. Because oblates and prolates exhibit specific long-term changes, our results suggest they may respond to different environmental forcing. However, it is widely assumed that diatom contribution to primary production and carbon export will decline in the mid and low latitudes because ocean warming will increase stratification and decrease nutrients supply to the surface ocean [START_REF] Bopp | Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models[END_REF][START_REF] Bopp | Response of diatoms distribution to global warming and potential implications: A global model study[END_REF][START_REF] Dutkiewicz | Winners and losers: Ecological and biogeochemical changes in a warming ocean[END_REF][START_REF] Fu | Climate change impacts on net primary production (NPP) and export production (EP) regulated by increasing stratification and phytoplankton communitystructure in the CMIP5 models[END_REF]. Such predictions are made because of an oversimplification of Margalef's mandala [START_REF] Margalef | Life-forms of phytoplankton as survival alternatives in an unstable environment[END_REF] assimilating diatoms to a single Plankton Functional Type (PFT) that thrives when nutrients and turbulence are high [START_REF] Fu | Climate change impacts on net primary production (NPP) and export production (EP) regulated by increasing stratification and phytoplankton communitystructure in the CMIP5 models[END_REF][START_REF] Kemp | The case of the diatoms and the muddled mandalas: Time to recognize diatom adaptations to stratified waters[END_REF][START_REF] Tréguer | Influence of diatom diversity on the ocean biological carbon pump[END_REF]. Therefore, current state-of-the-art biogeochemical models only consider oblates in their simulations and not prolates, but our results clearly show that diatoms, through morphological adaptations, are able to exploit a wider range of environment, which is also confirmed by (i) a Tara Oceans expedition study that has found comparable diatom diversity in oligotrophic and coastal regions [START_REF] Malviya | Insights into global diatom distribution and diversity in the world's ocean[END_REF] and by (ii) fossil records of the Mediterranean sea, where sapropels (one of the most carbon-rich sediments) formation is associated with highly stratified conditions and high diatom production and export [START_REF] Kemp | High diatom production and export in stratified waters -A potential negative feedback to global warming[END_REF]. Hence, we suggest that biogeochemical and earth-system models should consider diatom morphology by defining two PFTs (i.e. oblates and prolates) in their settings to improve their predictions on the consequences of climate change on primary production and carbon export.

Conclusions

By characterizing ADS in the North Sea with a trait-based approach, we have shown that this phenomenon is controlled by the niche-environment interaction that results from seasonal changes in morphological traits. Although four different hypotheses have been proposed to explain changes in cell shapes [START_REF] Karp-Boss | The Elongated, the Squat and the Spherical: Selective Pressures for Phytoplankton Shape[END_REF][START_REF] Litchman | Trait-Based Community Ecology of Phytoplankton[END_REF][START_REF] Naselli-Flores | Life in suspension and its impact on phytoplankton morphology: an homage to Colin S. Reynolds[END_REF][START_REF] Naselli-Flores | Invited Review -Fight on Plankton! or, Phytoplankton Shape and Size as Adaptive Tools to Get Ahead in the Struggle for Life[END_REF][START_REF] Smetacek | A watery arms race[END_REF], our results suggest that cell adaptation to turbulence and viscosity is of primary importance. Cell elongation gives a competitive advantage to diatoms in low-turbulence and low-viscous waters because it diminishes sinking velocity [START_REF] Padisak | Sinking properties of some phytoplankton shapes and the relation of form resistance to morphological diversity of plankton -an experimental study[END_REF]. Our results have also shown that elongation does not affect surface/volume ratio and therefore the ability of cells to absorb nutrients (Figure 1g-h). We did not find any evidence that elongation might enhance light capture despite some suggestions that it might (Naselli-Flores and Barone, 2011). Finally, although some studies have suggested that zooplankton can feed on large chains and elongated cells [START_REF] Djeghri | High prey-predator size ratios and unselective feeding in copepods: A seasonal comparison of five species with contrasting feeding modes[END_REF][START_REF] Visser | On the reorientation of non-spherical prey particles in a feeding current[END_REF], the positive correlation between copepod abundance and elongation that we found suggests predation may also be an important selective agent. Our study suggests that prolates and oblates are key adaptations to specific environments, conferring to species a specific niche. The resulting niche-environment interaction explains phenology and therefore the precise position of a species in the sequence of ADS. Understanding the relationships between morphological traits and species niche can help us to better anticipate community shifts and their functional implications for the ecosystems in the context of global climate change [START_REF] Litchman | Phytoplankton niches, traits and eco-evolutionary responses to global environmental change[END_REF]. To conclude, as both oblates and prolates exhibit different long-term changes in abundance, suggesting that they might also have different responses to climate forcing, we encourage biogeochemical and earth-system modelers to implement these two diatom groups to improve the assessment of marine primary production and carbon export.

viable. The authors also want to thanks Eric Lecuyer, Dimitra-Ioli Skouroliakou, Monica Michel-Rodriguez and Lucie Courcot for their helpful comments during the design of this study.

Figure 1 .

 1 Figure 1. Annual Diatom Succession (ADS) and morphological traits. (a) ADS characterized by a PCA applied on a species abundance matrix 366 days x 45 species. The first three principal components (PCs) are shown. (b-f) Relationships between the first three eigenvectors and morphological traits. Relationships between species groups and (b) morphological traits, (c) oblates and prolates, (d) height-diameter differences, (e) area/volume ratio and (f) δ. Blue area delineates the space occupied by the spring-dominated group, green area the space occupied by the autumn-dominated group and red area the space occupied by the LSSA group. (d-f) Size and colour of the bullets are scaled with (d) mean cell height-diameter difference, (e) mean area/volume ratio and (f) δ. (g-h) Relationships

Figure 2 .

 2 Figure 2. Relationships between day of maximum abundance and ecological optima/amplitude of each diatom and environmental parameter. (a) Theoretical ecological niche with Xopt representing species niche optimum and σ denoting niche amplitude assessed by the weighted standard deviation.The Y axis corresponds to the species abundance and the X axis to an hypothetical environmental parameter. Relationships between day of maximum abundance and niche optimum/amplitude for (b) SST, (c) euphotic depth, (d) PAR, (e) mixed layer depth (MLD), (f) silicate, (g) phosphate and (h) nitrate concentration. Blue lines correspond to the North Sea daily climatology of the parameter considered in the panel. The niche optimum can be read on the ordinate axis. The size of the bullet is a function of species' rank of the weighted standard deviation; larger bullets are associated with larger niche amplitude and therefore ranks. Oblates are in blue, prolates in red and oblates/prolates in green.

Figure 3 .

 3 Figure 3. Summary of the relationships between the first three eigenvectors, diatom shapes and the environment. Cell shapes are represented by cylinders scaled with species height and diameter (Supplementary TableS2). Oblates are in blue, prolates in red and oblates/prolates in green. Species names are those recorded in the CPR dataset. Blue area delineates the space occupied by the springdominated group, green area the space occupied by the autumn-dominated group and red area the space occupied by the LSSA group. Correlation between the first three principal components and the environmental parameters are displayed in Supplementary TableS5. More information on the principal components can be found in Figure1a.

Figure 4 .

 4 Figure 4. Long-term changes in oblate and prolate abundance (i.e. number per CPR sample) between 1958 and 2017. (a) Long-term monthly changes in total oblate and prolate abundance. Dashed lines highlight the months that were considered in the examination of long-term changes in oblates or prolates on panels (b), (c) and (d). (b-d) Long-term changes in the abundance of oblates in (b) April and (c) October and (d) of prolates in August. Dashed lines represent the long-term changes in diatom abundance from panel (a) and black bold lines show the smoothed long-term changes after applying a second-order simple moving average. Note that the exceptional nil abundance value observed for oblates and prolates in 1978 is below the y axis in the panels c and d.

Acknowledgements:

We would like to thank the owners and crews that have towed the CPRs on a voluntary basis for over 80 years contributing to one of the world's largest and longest ongoing ecological experiments. Without these early pioneers of citizen science and broad-scale volunteer monitoring projects this unique ecological dataset would never have been financially or logistically

Funding information: This work is supported by the Centre National de la Recherche Scientifique (CNRS), the Region Hauts-de-France and the Marine Biological Association of the UK.

Data sharing statement: Data used in the production of this manuscript are already freely available (see Materials and Methods section and Supplementary Text S2).

Supplementary Table S4). The first PC (26% of the total variance) showed a strong peak in spring at a time of weak stratification and low temperature, reflecting the spring bloom, and a smaller peak in late autumn at a time of decreasing temperature and stratification breakdown [START_REF] Caracciolo | Annual phytoplankton succession results from niche-environment interaction[END_REF] (Figure 1a). The second PC (25%) exhibited a similar pattern with two peaks, the lowest in spring and the highest in autumn, again reflecting both spring and autumn blooms. The third PC (15%) showed high positive values from late spring to the end of summer, identifying diatoms that were dominant during the stratified low-nutrients period at a time of higher zooplankton abundance (Supplementary Figure S2a-c). We therefore identify three key periods in ADS, confirming a pattern seen in temperate seas [START_REF] Caracciolo | Annual phytoplankton succession results from niche-environment interaction[END_REF][START_REF] Colebrook | Environmental influences on long-term variability in marine plankton[END_REF][START_REF] Colebrook | Continuous plankton records: relationships between species of phytoplankton and zooplankton in the seasonal cycle[END_REF].