Remote Sensing of Environment Trends in phytoplankton phenology in the Mediterranean Sea based on ocean-colour remote sensing

Seventeen years (1998 – 2014) of satellite-derived chlorophyll concentration (Chl) are used to analyse the sea- sonal and non-seasonal patterns of Chl variability and the long-term trends in phytoplankton phenology in the Mediterranean Sea. With marked regional variations, we observe that seasonality dominates variability re- presenting up to 80% of total Chl variance in oceanic areas, whereas in shelf-sea regions high frequency variations may be dominant representing up to 49% of total Chl variance. Seasonal variations are typically char- acterized by a phytoplankton growing period occurring in spring and spanning on average 170days in the western basin and 150days in the eastern basin. The variations in peak Chl concentrations are higher in the western basin (0.88 ± 1.01mgm − 3 ) compared to the eastern basin (0.35 ± 1.36mgm − 3 ). Di ﬀ erences in the seasonal cycle of Chl are also observed between open ocean and coastal waters where more than one phyto- plankton growing period are frequent (>0.8 probability). During the study period, on average in the western Mediterranean basin (based on signi ﬁ cant trends observed over ~95% of the basin), we show a positive trend in Chl of +0.015 ± 0.016mgm − 3 decade − 1 , and an increase in the amplitude and duration of the phytoplankton growing period by +0.27 ± 0.29mgm − 3 decade − 1 and +11 ± 7days decade − 1 respectively. Changes in Chl concentration in the eastern (and more oligotrophic) basin are generally low, with a trend of − 0.004 ± 0.024mgm − 3 decade − 1 on average (based on observed signi ﬁ cant trends over ~70% of the basin). In this basin, the Chl peak has declined by − 0.03 ± 0.08mgm − 3 decade − 1 and the growing period duration has decreased by − 12 ± 7days decade − 1 . The trends in phytoplankton Chl and phenology, estimated in this study over the period 1998 – 2014, do not reveal signi ﬁ cant overall decline/increase in Chl concentration or earlier/delayed timings of the seasonal peak on average over the entire Mediterranean Sea basin. However, we observed large regional variations, suggesting that the response of phytoplankton to environmental and climate forcing may be complex and regionally driven.


Introduction
Phytoplankton chlorophyll (Chl) and productivity exhibit variations in timing and magnitude in different regions of the oceans due to several processes affecting their growth and demise, such as incident solar irradiance, water column stratification, nutrient supply and grazing pressure, which may vary with latitude, as well as regional and local oceanographic conditions (Cushing, 1975). Temporal and spatial variations of phytoplankton Chl and production affect carbon fluxes and other elemental cycles in the oceans (Falkowski et al., 1998;Laufkötter et al., 2016) and can lead to upper trophic mismatch altering the function of marine ecosystems with important consequences for fisheries (Edwards and Richardson, 2004). Indeed, for many commercially-important fish species, their spawning season is closely tied with the phytoplankton seasonal cycle (Chambers and Trippel, 1997;Cushing, 1975;Kassi et al., 2018;Lo-Yat et al., 2011;Platt et al., 2003).
Seasonal dynamics of phytoplankton pigment concentrations in the oceans are indicative of variations in bio-physical processes such as primary production, grazing and/or structural changes driving phytoplankton diversity in response to environmental change in the ocean surface layer (Shevyrnogov and Vysotskaya, 2006). The seasonal bloom is one of the dominant features in the growth patterns of phytoplankton in the pelagic environments, particularly in temperate and cold oceans, and it responds primarily to changes in light and nutrient availability (Legendre, 1990;Sverdrup, 1953). Interannual and long-term variations in plankton phenology may be used as indicators to assess how changes in the marine environment propagate from primary producers to higher trophic levels. At high latitudes, the trophic transfer relies on a single seasonal peak of phytoplankton, providing the main energy source for zooplankton and fish production Racault et al., 2012). However, in temperate areas, a less intense secondary peak may occur in fall (Bosc et al., 2004;Mayot et al., 2017;Morel and André, 1991;Racault et al., 2016;Siokou-Frangou et al., 2010;Zingone et al., 1995).
The phase of the seasonal phytoplankton bloom can be characterized by a suite of phenological metrics: timings of initiation, termination, maximum, and duration of the phytoplankton growing period; and also amplitude of the phytoplankton Chl concentration (Platt and Sathyendranath, 2008;Racault et al., 2012). These metrics are valid for shelf areas where phytoplankton seasonal cycles are readily apparent in both satellite and in-situ time series (Corredor-Acosta et al., 2015;Fendereski et al., 2014;Friedland et al., 2015;Ji et al., 2010b;Jordi et al., 2009;Platt et al., 2009;Zhai et al., 2011). However, in these areas, local effects caused by the interaction of ocean processes with the seafloor and terrestrial sources can also play a fundamental role. Short spatial and temporal variability may account for a large fraction of total phytoplankton variance (Cloern and Jassby, 2008). Furthermore, the seasonal cycle may be masked by variability on longer timescales (Kahru and Mitchell, 2001). Non-seasonal variability in pigment concentrations have been associated with different physical drivers, such as synoptic fields of wind stress and sea level heights (Strub et al., 1990), frontal dynamics (Sokolov and Rintoul, 2007), or climatological variations (Lin et al., 2014;Thomas et al., 2012). In addition, in coastal shelf areas, factors regulating phytoplankton variability may be significantly different from the processes dominating in the open ocean. Signal variability in coastal areas can be driven by the effects of river plumes (Walker and Rabalais, 2006) or by topographically induced variability occurring in the vicinity of islands, headlands and straits (Arístegui et al., 1997;Gomez et al., 2017).
The Mediterranean is a marginal sea divided in two main sub-basins separated by the Sicily channel. It has been considered for many years as a 'small-scale ocean' presenting a wide range of physical processes including deep water formation and thermohaline circulation (Lacombe et al., 1981;Robinson and Golnaraghi, 1993;Volpe et al., 2012). External sources of water are the inflow of the Atlantic water jet in the surface layer through the Strait of Gibraltar, the Black and Marmara Seas through Bosporus Strait and freshwater rivers inputs mainly in the northern western side (The MerMex Group, 2011). Biogeochemically, the Mediterranean Sea presents an apparent west-east gradient of increasing oligotrophy. A deficiency of phosphorous relative to nitrogen increases eastwards (Béthoux et al., 1998). The dynamics of oceanic phytoplankton in the western basin are characterized by a winter phytoplankton growing period similar to that of temperate seas whereas in summer, low Chl concentration generally prevails (i.e. Marty et al., 2002). Phytoplankton dynamics in the eastern basin resemble that in subtropical waters, with low Chl variations throughout the year. D'Ortenzio and Ribera d'Alcalà (2009) and Lavigne et al. (2015) refer to this region as a 'no-bloom' region because changes in Chl values are limited. Due to the low nutrient availability, phytoplankton assemblages in the Mediterranean Sea are generally dominated by small picoplankton (< 2 μM) for most of the year although the contribution of larger nanoplankton and microplankton cells have been shown to increase during winter and in coastal areas (Ignatiades et al., 2009;Vidussi et al., 2001).
Several studies have addressed the seasonal and interannual variability of phytoplankton at regional and basin-wide scales. For instance, Mélin et al. (2011) analysed the phenology in the Adriatic Sea using monthly ocean-colour composites and signal decomposition applying a simplified X-11 algorithm for the period 1997-2004. Using empirical orthogonal functions analyses, Volpe et al. (2012) investigated the seasonal response of phytoplankton in practically the entire Mediterranean Sea during the period 1998-2006. Some studies have also examined the biogeography and changes in magnitude in Chl concentration in relation to changing environmental conditions during different time periods D'Ortenzio and Ribera d'Alcalà, 2009;Mayot et al., 2016) but, with some exceptions (e.g. Lavigne et al., 2013), less is known about the long-trends of phytoplankton Chl and phenology. Barale et al. (2008) reported symptoms of an increased nutrient-limitation, resulting from reduced vertical mixing due to a more stable stratification of the basin, in line with the general warming trend of the Mediterranean. While warming of Mediterranean Sea surface appears to be a consistent pattern (e.g. Belkin, 2009;Rixen et al., 2005;Skliris et al., 2012), some uncertainties remain on longterm phytoplankton response that it is the result of complex processes, including variations in nutrient sources and availability, as well as changes in food web interactions.
The present work aims to characterize patterns of variability and trends of phytoplankton Chl and phenology in the Mediterranean Sea. Our study uses a regionally-tuned Chl product for the Mediterranean Sea based on the OC-CCI dataset, which is the longest, continuous, climate quality-controlled, 17-year time series (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) of fine space-time resolution remotely sensed ocean-colour data. First, we assess the contribution of seasonality to overall variability and the longterm tendencies of deseasonalized Chl concentration. Then, using a threshold-based phenology approach, we identify and characterize spatial patterns in phytoplankton phenology over the entire Mediterranean Sea basin. Finally, we calculate the trends in the estimated phytoplankton phenological metrics for the period of study. The observed changes in phytoplankton Chl and phenology are key to understand functional and structural aspects of marine communities.

Ocean-colour remotely sensed data
The present study is based on sea surface Chl concentration (mg m −3 ) Level-4 reprocessed product, obtained from the EU Copernicus Marine Environment Monitoring Service (CMEMS) at 8-day and 1-km resolution, from the website http://marine.copernicus.eu/ (OCEANCOLOUR_MED_CHL_L4_ REP_OBSERVATIONS_009_078). The analysed series covers the period from January 1998 to December 2014 for the Mediterranean Sea (30-to 46°N and 6°W to 37°E, Fig. 1). This ocean-colour data record is a regionally-tuned reprocessing of the climate-quality, error-characterized, and bias-corrected merged product of multi satellite observations (SeaWiFS, MODIS-Aqua and MERIS sensors) initially developed by the European Space Agency Ocean-Colour Climate Change Initiative Program (ESA OC-CCI) Sathyendranath and Krasemann, 2014). The Chl product available from CMEMS has been tailored to the Mediterranean region by using the regional algorithm MedOC4 (Mediterranean Ocean-Colour 4 bands, Volpe et al., 2007) for Case-1 waters and the AD4 algorithm (ADriatic 4 band, Berthon and Zibordi, 2004; for Case-2 waters. In this product, the merging of Case-1 and Case-2 information was performed following D' Alimonte et al. (2003). In practice, the CMEMS processor ingests the OC-CCI remote sensing reflectance, which is the result of a merging procedure that accounts for the inter-sensor bias among different sensors and then applies the specific regional algorithm. The OC-CCI processing removes inter-sensor bias by band-shifting (Mélin and Sclep, 2015) and bias-correcting the MODIS and MERIS bands and values to the SeaWiFS reference ones and finally, by merging the individual sensors datasets with a weighted averaging (Mélin et al., 2016;Mélin et al., 2017). The overall correlation between measured and estimated Chl, reprocessing time series REP, L4 product in Mediterranean Sea yield a correlation coefficient of r 2 = 0.764 (Pitarch et al., 2016), and the specific correlation for Case-2 waters is r 2 = 0.56 (D' ALimonte and Zibordi, 2003).
In order to reduce missing data, Chl values were first re-gridded to a 4 km regular grid by averaging all available data points within each new, larger pixel. Then, remaining gaps were filled by applying linear interpolation scheme in the sequential steps: longitude, latitude, and time (Racault et al., 2014b). Specifically, the gaps were filled with the average value of the surrounding grid points along the indicated axis.
The averaging window had a width of three points and the surrounding points were weighted equally. Along the indicated axis, if one of the points bordering the gap was invalid, it was omitted from the calculation. If the two surrounding points were invalid, then the gap was not filled. Finally, the data were smoothed by applying moving average filter of previous and next time step. The study area with its bathymetry and main sea surface circulations patters is shown in Fig. 1.

Relative contribution of seasonal and non-seasonal variations in chlorophyll
Census X-13 seasonal adjustment methodology, an improved version of X-11, is used to remove the seasonal signal from the Chl time series and to investigate the temporal variation patterns (methodology applied following the X-13 seasonal adjustment equations from the U.S. Census Bureau, 2017). Briefly, the X-13 technique is based on an iterative bandpass-filtering method that relies on the application of successive filters. Each component is estimated with an iterative procedure consisting of an alternate computation of the trend component from the seasonally adjusted series, and of the seasonal component from the series corrected for the trend. A principal feature of the X-13 decomposition is that the seasonal term is defined locally in time allowing year-to-year variations, and therefore has a variable amplitude accounting for interannual variations in seasonality. The procedure assumes that an original time series X(t), here a 8-day Chl composites, can be decomposed into three components: as X(t) = S(t) + I(t) + T(t) where S, I and T, represent, respectively, a seasonal signal, and an irregular (or high frequency residual), and a low frequency (> 1 year) or interannual component (Shiskin, 1978).
Trends in Chl time series were calculated by the Sen-slope estimator to the T(t) component (see details in Section 2.4). The relative contribution of each of the three components (S,I,T) of Chl is calculated as in Vantrepotte and Mélin (2009). The total variance σ x 2 of the original time series can be written as σ x 2 = σ 2 S + σ 2 I + σ 2 T + 2cov(S, I, T) where σ 2 S , σ 2 I and σ 2 T represent the variance associated with the S, I and T components respectively, and cov(S, I, T) represents the covariance terms between these components. The covariance term 2cov(S, I, T) accounts only for few percentage of the total variance (typically, 5% in absolute value). Therefore, the relative contribution p z 2 of each component to the total variance (with z = S, T or I) is approximated as ρ z 2 = 100 σ z 2 (σ 2 S + σ 2 I + σ 2 T ) (expressed in %). Time series with > 12% of missing values were excluded from the analysis.

Estimation of phenological indices
We applied the algorithm of Racault et al. (2017) to calculate a series of phenological indices: timings of initiation (b i ), peak (b t ), termination (b e ), duration (b d ) and peak amplitude (b a ) of up to two phytoplankton Chl peaks during each year. The phenology algorithm uses Sea Surface Temperature (SST) seasonality to define two time intervals (i.e., a SST warming phase and a SST cooling phase) during which up to two phenological growing periods can be estimated. In the Mediterranean Sea, the main phytoplankton growing period is associated with the minimum winter SST. In winter, cool, dense surface waters deep and causes a vertical mixing convection process in the water column that brings nutrients to the surface photosynthetic layers (Marty et al., 2002). In spring and summer, there is a setting of stratification in the water column due to the increase in SST. A second growing period may be observed in fall when water-column becomes less stratified, generally associated with storms or wind events. The fall growing is generally shorter than the spring one, and Chl concentration remain lower (Morel and André, 1991). Similar approaches to the one we use in the present work have been successfully applied at global scale (including tropical and subtropical regions (e.g. Kostadinov et al., 2017;Racault et al., 2017), as well as at regional scale in tropical regions such as the Red Sea (e.g., Gittings et al., 2018;Racault et al., 2015).
In the present study, the method to estimate the phenological indices is based on a long-term Chl median plus 20% pixel-by-pixel threshold criterion and calculation of the derivative of the cumulative sum of Chl anomalies. The selection of this threshold depends on the magnitude of the seasonal Chl peak with respect to baseline variations. For example, Siegel et al. (2002) find little differences between using values ranging from 5 to 30% in the North Atlantic spring bloom. In our case, the threshold of 5% yields to the detection of multiple small Chl peaks and, therefore, we chose the median plus 20% threshold to avoid the detection of these small peaks, especially in coastal regions. The timing of initiation b i was obtained based on the first time step when the cumulative sum of Chl changes signs (i.e., from positive to  (Millot and Taupier-Letage, 2005;Pinardi et al., 2006;Robinson et al., 2001;Siokou-Frangou et al., 2010). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) negative), and the timing of termination b e is calculated as the previous time step before falling below the threshold criterion. To identify b i and b e , the Chl concentration needs to fulfil the latter conditions during at least two consecutive time steps. The duration b d is the time elapsed between b i and b e . The timing and amplitude of the Chl peak, b t and b a respectively, and correspond to the time and amplitude at which the Chl value reaches the highest signal. Lastly, the method of detection of the timing of events is centred around b t , hence b i and b e times can span the calendar year. In any particular year, if two phytoplankton growing periods are estimated, the main growing period is defined as the one with highest Chl peak whereas the second growing period is defined as the one with lower Chl concentration. Further details about the method implementation are provided in Racault et al. (2015) for Chl climatology and Racault et al. (2017) for interannual Chl time series. A schematic representation of the phenology method is shown in Fig. 2.
In addition to computing phenology maps, the seasonal climatology of ten selected regions was estimated from the average of 300 detrended time series representing a surface area of 4600 km 2 . We selected these regions based on previous studies looking at spatio-temporal features of phytoplankton variability in the Mediterranean Sea which showed that these areas exhibit distinct seasonal patterns D'Ortenzio and Ribera d'Alcalà, 2009;Lavigne et al., 2017;Macias et al., 2017;Mélin et al., 2017).

Trend calculation
Trends on the 8-day T component (de-seasonalized Chl time series obtained from X-13 decomposition) and on the yearly phenological indices time series were calculated pixel-by-pixel. The existence of linear trends and their significance was evaluated with the Sen slope estimator, which calculates the slope and the constant of the linear models as the median of those derived by linear interpolations of any pair of available data (Gilbert, 1987;Sen, 1968). A Mann-Kendall statistics test is applied to each pixel to identify the statistically significant trend values at a 95% confidence level (Salmi et al., 2002). The use of this non-parametrical test is suitable for non-normally distributed data and has been previously used in the trend examination of remote-sensing Chl time series Coppini et al., 2013;Erasmi et al., 2014;Kahru et al., 2011;Vantrepotte and Mélin, 2011). As missing data could mask the trend results, only de-seasonalized Chl time series with at least 90% of the pixels were used and linear interpolation was performed on the rest of the data. A minimum of nine years (~50% of the total number of years 1998-2014) with bloom indices was set as a requirement to perform the regression (i.e., if phenological indices were available for less than nine years then the trend was not calculated). Trends in phenological indices are given in days per decade, which have been calculated by multiplying the phenology timing index estimates by eight based on the 8-day resolution of the ocean-colour composites.
Trends in Chl concentration estimated from ocean-colour products that includes observations from single sensor, and from multiple sensors (i.e., merged products) have shown to give comparable results, and that changes in sampling rate have limited influence on the estimation of long-term trend in Chl . For the calculation of phytoplankton phenological metrics, the influence of the sampling rate (i.e., number of observations) has been estimated to introduce a bias < 2% and RMSE < 10% when the observation coverage is between 80 and 100% (Racault et al., 2014b). In the Mediterranean Sea, the individual sensor (SeaWiFS, MODIS, MERIS) can provide coverage within this range and allow us to estimate phenological metrics and trends with low uncertainty.
Considering the length of the ocean-colour record 1998-2014, we would like to note that our analysis is representative of the variations during the last two decades. Longer datasets (i.e. > 30 years) may be necessary to identify climate change driven trends . Phenological metrics are sensitive to gaps in the time series, which can result in error and bias in the estimates (Racault et al., 2014b). However, compared with other regions, the Mediterranean Sea tends to benefit from low cloudiness (Sanchez-Lorenzo et al., 2017). In addition, we also note that phenological changes may not be effectively detected when noise levels are larger than the signal amplitude (Verbesselt et al., 2010). Uncertainty in the estimation of phenological metrics can be improved by the use of preprocessing techniques such as gap-filling and smoothing (Ferreira et al., 2014;Racault et al., 2014b). We use 8-day composites of the OC-CCI dataset, which result in quasicontinuous time series (no gap) in the Mediterranean Sea. While averaging over 8 days reduces the resolution at which events in the phytoplankton growing period can be estimated, it does not significantly affect the spatial pattern of the estimated trends, as shown by Henson et al. (2018).

Seasonal and non-seasonal components of chlorophyll variability
The S, I and T components of Chl variability derived by the X-13 methodology are displayed in Fig. 3. Seasonality prevails in most of the oceanic areas, representing up to 80% of the total variance (Fig. 3a). Interestingly, most of the northern Mediterranean (i.e. above latitude 39 o N) presents higher high frequency (irregular) variability. Low frequency variations are particularly intense in the Ibiza and Mallorca channels, in the northern Aegean Sea and Otranto channel (Fig. 3c). Generally, Chl variations in coastal areas present a more complex scenario compared to oceanic regions. This is especially evident in the Northern Adriatic and along the Po River plume path, and also in the northern Alboran Sea, where the contribution of high frequency to total variance is > 60%. Regions with a cyclonic circulation and a blooming dynamic also present strong part of irregular variability (NWMS, South Adriatic and Rhodes gyre; Fig. 3b). Fig. 4 and Table 1 shows the seasonal, irregular and interannual Chl variations at the ten selected oceanic and coastal areas. The oceanic areas include the North western Mediterranean Sea (NWMS), the South Adriatic pit, the Levantine Sea, the Ionian Sea and the South Western Mediterranean Sea (SWMS) in the Algerian basin; whereas the coastal areas are represented by a subset of pixels of the Gulf of Lion, the North Adriatic coast, the Nile estuary, the Gulf of Gabes and the northern Alboran Sea. Variability in oligotrophic oceanic areas such as the Levantine and Ionian Seas (Chl < 0.06 mg m −3 ) and in the more productive SWMS (0.15 ± 0.02 mg m −3 ) is dominated by seasonality (ρ s 2 > 73% of the variance, see Table 1). Conversely, high-frequency . Black squares at each figure delimits the ten areas over which temporal decomposition is shown in Fig. 4 (zones, z, 1 to 10).
variations constitute an important part of total variance in the NWMS and the South Adriatic oceanic areas. Among the coastal areas, the Gulf of Gabes is the one that most resembles the oceanic areas, except for its higher Chl values (0.83 ± 0.54 mg m −3 ), this region is dominated by the seasonal component (ρ s 2 = 71% of the variance). Other coastal areas such as the Gulf of Lion and Nile Estuary present a lower seasonality, but still, it dominates total Chl variance (ρ s 2 = 58 and 57% respectively). Contrastingly, the contribution of irregular variations in the North Adriatic (ρ I 2 = 61%) and Alboran (ρ I 2 = 64%) largely exceeds seasonality. Interannual variations represent a small proportion of the observed variability (ρ T 2 generally between 3 and 7%, and exceptionally 12% in the Northern Adriatic Sea) even though the contribution of this component is higher in coastal areas. Fig. 5 shows the global distribution of the long-term Chl trend (T). Most regions in the western Mediterranean and Adriatic Sea present a positive trend. The average value in the western basin is +0.015 ± 0.016 mg m −3 decade −1 and it intensifies in the NWMS region to +0.107 ± 0.004 mg m −3 decade −1 and in the northern Alboran Sea (+0.138 ± 0.007 mg m −3 decade −1 , Table 1). Large positive trends are found in the Adriatic Sea where the mean increase is +0.047 ± 0.085 mg m −3 decade −1 , and maximum values are found along the north-western shore (+0.241 ± 0.022 mg m −3 decade −1 , Fig. 4. Mean seasonal (black), irregular or high frequency (blue) and interannual trend (red) components, of Chl variability over the ten selected oceanic and coastal areas. The location of each zone z1 to z10 is displayed in Fig. 2. Note that different scales are used in some subplots. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 1 Mean Chl; variances contribution explained by the seasonal (ρ s 2 ), irregular (ρ I 2 ), and interannual (ρ T 2 ) components; slope of the interannual Chl variability (T component) and its associated r 2 value, as derived from the X-13 procedure of the ten selected oceanic and coastal areas.

Basin-wide phytoplankton phenology
The spatial distribution of mean phytoplankton phenological indices in the Mediterranean Sea for the period 1998-2014 is presented in Fig. 6 (see Suppl. Fig. 1 for Standard Deviation and Coefficient of Variation values). The seasonal cycle in oceanic areas (> 200 m) is typically characterized by a single growing period initiating early, around the 8-day period −8 (i.e., corresponding to the beginning of November) and terminating between the 8-day periods +13 and + 20 (corresponding to the end March and mid-May respectively). Largest spatial differences in b i and b e are observed in the NWMS where the main growing period is delayed by approximately 60 days with respect to other surrounding areas (the growing period starts in mid-January and ends in late May; Fig. 6a-b). Chl peaks from late January to April, with clear regional differences (Fig. 6c). The growing period duration is shorter in the eastern basin (b d~1 50 days) compared to the western basin (b d~1 70 days) (Fig. 6d). In the same way, b a is lower in the Only 95% significant Theil-Sen trends (p < 0.05) are shown. Black squares at each figure delimits the ten areas over which trend estimation is shown in Fig. 4 and Table 1 (zones, z1 to z10). Grey contour line delimits coastal waters (200 m isobath). eastern basin, with Chl values expanding from 0.1 to 0.5 mg m −3 (mean b a~0 .35 ± 1.36 mg m −3 ) compared to the western basin, where mean values range from 0.5 to 2.5 mg m −3 (mean b a~0 .88 ± 1.01 mg m −3 ; Fig. 6e). A spatial variation is seen in the NWMS where the main seasonal peak is more intense but shorter (b a values of up to 3 mg m −3 and b d spanning~120-150 days) compared to the rest of the western oceanic region (Fig. 6d-e).
Shelf regions (< 200 m) present similar phenological patterns to the oceanic areas, but with some particularities. For example, in the Gulf of Gabes, the Adriatic coast of Italy, and in the eastern Mediterranean coast expanding form Egypt to Lebanon, the timings of b t and b e shifted earlier by about 40 days compared to the oceanic areas i.e. b t December-January and b e in February-March (Fig. 6b-c). These coastal regions also present shorter growing periods compared to the open ocean areas (b d~1 00 days, Fig. 6d) and higher peak amplitude (b a~5 -7 mg m −3 , Fig. 6e). Coastal waters are the areas where the probability of occurrence of an autumn growing period is higher (> 0.7 probability of occurrence). Nevertheless, a second peak is not limited to neritic regions. Some oceanic areas such as the NWMS, the Strait of Gibraltar and the northern Aegean Sea also present secondary growing periods (Fig. 6f). Fig. 7 shows the mean Chl value during the period 1998-2014 over the Mediterranean Sea (Fig. 7a) and the 8-day seasonal climatological cycle of five oceanic (Fig. 7 z1-z5) and five coastal (Fig. 7 z6-z10) areas during that period. Coastal areas show higher b a and background Chl values compared with oceanic areas. Coastal climatological cycles also present higher standard deviation (see dashed lines). Marked seasonal climatological differences are observed between the oceanic regions. For example, the NWMS (Fig. 7 z1) presents a well-defined seasonal peak in April (1.3 ± 0.2 mg m −3 ) whereas a weaker (~0.1 mg m −3 ) and earlier growing period occurring in February in the Levantine, Ionian Sea and SWMS (Fig. 7 z3-z5). The South Adriatic Sea (Fig. 7 z2) present two intense phytoplankton growing periods, one in April (0.30 ± 0.03 mg m −3 ) and the second one in November (0.16 ± 0.01 mg m −3 ).
With the exception of the Gulf of Lion and Alboran Sea, coastal regions present a large interannual variability (standard deviation of 0.7, 1 and 0.7 mg m −3 for North Adriatic, Nile estuary and Gulf of Gabes respectively). The North Adriatic Sea (Fig. 7 z7) is markedly different to other coastal areas, exhibiting two phytoplankton growing periods (May and November) with similar characteristics in duration and intensity (~2 mg m −3 ). The Nile estuary (Fig. 7 z8) has its Chl peak around February (1.11 ± 1.25 mg m −3 ) and is the coastal area with the highest standard deviation. Finally, the Gulf of Gabes (Fig. 7 z9) presents one single peak (1.6 ± 1.4 mg m −3 ) with the particularity to occur earlier compared to the rest of the coastal areas, starting in June, peaking in September and terminating in February. Finally, even though the Alboran Sea (Fig. 7 z10) present relatively large variations in the seasonal cycle, the interannual variations are low, as reported in the low Standard Deviation.

Trends in phytoplankton phenology
Trends in phytoplankton phenology metrics are displayed in Fig. 8. A positive trend in b i is observed in most of the Mediterranean Sea (~70% of the pixels) indicating a mean delay of +9 ± 5 days decade −1 . In shelf regions and some areas such as the Alboran Sea or the Rhodes gyre the main growing period starts earlier (mean − 6 ± 5 days decade −1 ; Fig. 8a). b e is delayed mainly in the south-western Mediterranean Sea by +9 ± 7 days decade −1 on average, whereas an early b e is observed in the eastern basin (mean − 6 ± 5 days decade −1 ; Fig. 8b). As shown in Fig. 8c, b t is generally earlier with a mean value of +12 ± 8 days decade −1 in the 71% of the Mediterranean basin. In the eastern and north-western Mediterranean Sea b d has decreased by −12 ± 7 days decade −1 on average, and increased mainly in the Algerian basin (+11 ± 7 days decade −1 ; Fig. 8d). The trend in b a , displayed in Fig. 8e, shows an overall increase in 86% of the western basin of +0.27 ± 0.29 mg m −3 decade −1 , reaching +1 mg m −3 decade −1 in the NWMS and Alboran Sea. The exception to this positive pattern is observed in coastal areas such as the Gulf of Lion and less intensely south of the Balearic Islands (mean decrease in b a of −0.15 ± 0.24 mg m −3 decade −1 ). In contrast, the eastern basin presents a decreasing b a of −0.03 ± 0.08 mg m −3 decade −1 in 66% of the basin. However, some positive trends are observed in coastal areas, particularly in the Nile Estuary and the north Aegean Sea (mean increase of +0.08 ± 0.23 m −3 decade −1 ).

Patterns of seasonal and non-seasonal phytoplankton chlorophyll variation
In the present analysis, seasonality is shown as the main contributing component to Chl variability in the Mediterranean Sea, representing from 50 to 80% of total variance. Lower values are observed in oceanic compared with coastal areas where high frequency episodes may be locally important (Fig. 3a). Our estimations are consistent with the values reported by Vantrepotte and Mélin (2009) who estimated that the contribution of seasonality in the global ocean was 64% of the total variance, and obtained values > 80% in the Mediterranean Sea. The latter value is consistent with our estimations for low Chl oceanic regions such as Levantine and Ionian Sea (Table 1). However, it is less representative of the northern regions of the Mediterranean Sea (i.e. NWMS and North Adriatic). The coastal areas present a higher proportion of high frequency episodes (irregular variability) as shown in Fig. 3b. Seasonal Chl peaks in these areas are highly irregular and characterized by a series short-lived (1 to 3 weeks) enhanced Chl episodes during the productive season. Moreover, these events are not restricted to the main growing season (i.e. winter) but, instead, are also observed during other periods. This suggests that Chl is driven by other forcing mechanisms such as intense winds, river outflows and exchanges with adjacent Seas (i.e. in Alboran and Northern Aegean Seas). Indeed, high frequency variability is enhanced in the Northern Mediterranean where winds are more intense during winter (Zecchetto and De Biasio, 2007).
Low frequency variations only represent a small fraction of variance (i.e. < 15%) and interestingly, they are intensified in key bottleneck areas. The channels of Ibiza, Mallorca, Otranto, Sardinia-Sicily and northern Aegean with its connection with the Marmara Sea, present higher variability. According to Astraldi et al. (1999), the channels and straits of the Mediterranean Sea are crucial for intercepting the major cells of the Mediterranean Sea circulation and often present complex circulation patterns (Fig. 3c). In the case of the Balearic channels, it is a transitional area where strong and highly variable meridional exchanges take place between the north and the south of the western Mediterranean (Lopez-Jurado et al., 1996). Jordi et al. (2009) observed deviations in the non-seasonal Chl patterns in the region influenced by the Balearic current and southeast of Mallorca and Ibiza. Also, Pinot et al. (1995) reported that winter intermediate waters formed further north, accumulate to the north of the channels in late spring obstructing the circulation through the Ibiza Channel. This process is determinant for the regulation of the interannual meridional water exchange in the sub-basin and has consequences in the distribution on the functioning of the marine ecosystem (Balbín et al., 2014;Pinot and Ganachaud, 1999). In the case of the northern Aegean Sea, previous studies have reported that variations in water input from the Black Sea can have the same influence, or even higher, compared with variations that may be observed in relation to water discharge from any Mediterranean rivers (Zervakis et al., 2004). Exchanges through the Marmara Sea have also been reported to present a clear interannual component having the similar influence compared with water discharge from large rivers discharges (Macias et al., 2018;Maderich et al., 2015). For instance, an increase in nutrient availability is associated with an increase in Chl amplitude or peak.
Low frequency shifts in phytoplankton abundance and composition in the Mediterranean Sea have been attributed to long-term alterations in the nutrient stoichiometry (Mercado et al., 2005). It has been observed that changes in the availability of reduced N and the N:P-ratios in inorganic and organic forms could be important factors in the control of bloom succession and ecosystem organization (Pujo-Pay et al., 2006). These nutrient variations can be motivated by climate driven variations. For example, decadal changes in winter cooling or evaporation can vary the intensity of vertical motions controlling the exchange of properties between the surface and deeper layers can influence the uplift of mineralized nutrients in the photic layers enhancing biological production (Polovina et al., 1995). García-Comas et al. (2011) observed decadal periodicity in secondary production forced by winter hydrographic conditions in the NWMS. Gregg and Conkright (2002) observed, for the global ocean, decadal changes in ocean Chl likely to be associated with climate oscillations. Furthermore, temperature driven variations in top-down controls of phytoplankton have been regarded as an important climate-related factor in the NWMS (Molinero et al., 2005). Indeed, large scale teleconnection patterns such as the North Atlantic Oscillation (NAO) are found to be the most important modes of atmospheric variability related to Chl concentration and distribution in the Mediterranean Sea (Basterretxea et al., 2018;Katara et al., 2008). Fig. 7. a) Mean Chl in the Mediterranean Sea over the period of study. Black squares delimit the ten areas from which seasonal climatology was derived. White contour lines separate coastal waters from offshore areas (200 m isobath). Panels z1 to z10: Mean climatology of 8-day Chl time series for a composite year, relative to the period 1998-2014, of oceanic (z1-z5) and coastal areas (z6-z10). Note that different scales are used in some subplots.
However, our analysis shows more spatially localized differences in the long-term trends (Fig. 3c), which appear to be closely related with hydrographic variations, such as changes in eddy kinetic energy or intensified exchange through channels and straits. Large interannual variations in eddy kinetic energy at regional scale have been observed in the Mediterranean Sea (Pujol and Larnicol, 2005). Sannino et al. (2009) have discussed the importance of decadal variations of the Atlantic water inflow through the Gibraltar Strait in determining decadal variations in the western Mediterranean circulation patterns, water column stratification and convection events; and Civitarese et al. (2010) have emphasized the relevance of decadal scale variations through the Otranto channel on de biogeochemistry of the Adriatic Sea.
In the present analysis, the basin scale differences in phytoplankton Chl concentration are marked and appear to be enhanced associated with the positive trend in the western basin and in the Adriatic Sea (Fig. 5). Increasing Chl concentration in the western basin was also observed by Vantrepotte and Mélin (2009) and further confirmed in more recent studies Coppini et al., 2013). Nevertheless, sub-basin and regional differences appear to stand out from this general tendency. For example, changes are intensified along the NWMS, in the Alboran Sea and along the path of the Algerian Current in the SWMS region. Increased phytoplankton Chl concentration in the NWMS (+0.107 ± 0.004 mg m −3 decade −1 , Table 1) has been investigated by Marty et al. (2002) who report an increase in small-sized phytoplankton in response to the lengthening of the summer stratification. A marked positive trend is also observed in the northern Adriatic (+0.241 ± 0.022 mg m −3 decade −1 , Table 1). This trend differs from previous studies, which reported Chl decrease in the Adriatic associated with a decline in productivity due to a reduction in precipitations and a decrease in nutrient supply from river run-off Mélin et al., 2011;Mozetič et al., 2009). Discrepancies with these studies may be explained by the inclusion in the present analysis of an additional six years (i.e. 2009-2014) displaying an important enhancement in phytoplankton Chl concentration. The main driver of this recent increase is uncertain, but it could be related to an increase in the Po River flow including higher phosphates and dissolved nitrogen concentrations (Giani et al., 2012). Long-term variations are typical from large river systems which are affected by climatic factors (e.g. Donner and Scavia, 2007;Goolsby and Battaglin, 2001) and the linear trends used for open waters may not capture well the overall trend.
In summary, Fig. 5 reveals that trends in Chl concentration of opposite sign are observed between the western and eastern basin. It is plausible that, these differences could be attributed to the different forcing mechanisms that regulate phytoplankton productivity in both basins. The western basin is more responsive to short-term processes such as wind pulses, deep winter convection and river discharges. Conversely, as suggested by Siokou-Frangou et al. (2010), the functioning of the eastern basin resembles that of subtropical regions. This region is expected to be particularly affected by long-term changes in evaporation and precipitation that can influence vertical mixing and, subsequently, modify nutrient availability (Ozer et al., 2017;Schroeder et al., 2017).

Patterns of phytoplankton phenological variation
The underlying mechanisms driving the variability of phenological indicators (Fig. 6) at each region can be complex (Thackeray et al., In the open ocean, the development of an upper mixed layer following the summer solar heating, and the convective processes caused by winter cooling and strong winds are the main driving forces of the seasonal signal in phytoplankton Chl concentration (Lavigne et al., 2013;Siokou-Frangou et al., 2010;Tanhua et al., 2013). However, variations in open waters, the phenology appear to be more closely related to regional scale than to basin scale processes. In particular, three regions have been shown to present distinct phytoplankton phenological characteristics: the NWMS, the Alboran Sea and the Adriatic Sea. The NWMS is a well-studied area displaying an intense seasonal peak (> 1 mg Chl m −3 , Fig. 7 z1) that is generated by the nutrients provided by deep water convection processes driven by cold and dry winds causing loss of surface water buoyancy and mixing with large depths (Macias et al., 2018;Mayot et al., 2016;Severin et al., 2014). In the case of the Alboran Sea, productivity is regulated by the input of Atlantic waters through the Strait of Gibraltar. Mixing in the Gibraltar sill, the development of strong geostrophic cross-frontal circulation and upwelling of the Atlantic jet near the Spanish coast increase the fertility of these waters (Estrada, 1996;Oguz et al., 2014;Solé et al., 2016). The influence of Atlantic waters may extend beyond the Alboran Sea by indirectly enrichening the waters of the Algerian current through the generation of instabilities, as suggested by Mayot et al. (2016). Finally, the phenology of the Adriatic Sea is driven by interactions between the outflows of major river discharges and wind induced convection processes (e.g. Gačić et al., 2002;Mauri and Poulain, 2007). These three areas present high Chl amplitude but short duration in the spring growing period. Also, the probability of occurrence of a secondary growing period during water column de-stratification in autumn is higher than in other open ocean areas (Fig. 6 d-f).
There are some differences in phenological metrics between the south-west and the eastern basin. In the specific case of the SWMS in the Algerian basin (Fig. 7 z5) seasonal patterns are similar to those in the eastern Mediterranean, as previously observed by D'Ortenzio and Ribera d'Alcalà (2009). Indeed, even the contribution of seasonality to total Chl variability is relatively similar (73% contribution, Table 1). However, the main difference is that the SWMS is a more productive area fueled by a shallower nitracline, and weaker mixing efficiency in the eastern basin, as reported by Lavigne et al. (2015). Furthermore, the phytoplankton growing period last about~30 days more in the SWMS compared with the oligotrophic areas of the eastern basin ( Fig. 6 c & d). In the ocean eastern oligotrophic areas such as the Ionian (Fig. 7 z4), and Levantine Sea (Fig. 7 z3), the increase in winter Chl concentration is very low (~0.1 mg m −3 ). Barbieux et al. (2018) reported that in these areas, phytoplankton Chl is essentially constant throughout the year. However, field evidence reveals that the eastern Mediterranean is not a stationary system but, instead, it presents significant seasonality with primary production and phytoplankton Chl generally increasing in response to winter phosphate and/or nitrate increase, particularly in coastal areas (Azov, 1986;Balkis, 2009;Ignatiades et al., 2002).
Coastal areas are shown to display different phenological pattern compared to the open ocean regions. Fueled by land sources and/or by coastal upwelling (Estrada, 1996), the amplitude of the main growing period b a in coastal areas exceeds that in oceanic waters (Figs. 6 & 7a). Secondary growing periods are also frequent in these areas such as the Northern Adriatic, the Gulf of Lion, the Nile delta, the northern Aegean Sea and the Gulf of Gabes (Fig. 6f). The first four areas are Regions Of Fresh water Influence (ROFI) regulated by the discharge cycles of the Po, Rhone, Nile and other smaller flows such as the Thracian rivers (Nestos and Aliakmon rivers). However, while the Po and Rhone rivers often display discharge peaks in fall and spring related to precipitation and snow melt periods (Mariotti et al., 2002;Montanari, 2012;Sempéré et al., 2000), this is not the case of Thracian and Nile rivers (Milliman and Farnsworth, 2011). Regarding the upwelling influence, it will be particularly important in coastal areas where wind fields pulses are predominant. These areas are characterized by local decrease in SST accompanied by a deep-water mixing, which brings nutrient-rich deep waters into the coastal euphotic layer associated with an increase in phytoplankton Chl. These mainly occurs in the Gulf of Lion, west coast of Greece in the Aegean Sea and in the Algerian Coast (Bakun and Agostini, 2001). In the shallow Gulf of Gabes, frequent tidal movements can cause vertical mixing and increase productivity levels (Sammari et al., 2006). In coastal regions displaying autumnal growing periods (see Fig. 6f) other factors providing new nutrients, such as surface runoff from terrestrial inputs or sediment resuspension, can be regarded as main drivers of phytoplankton variability as observed by Arin et al. (2013).

Trends in phytoplankton chlorophyll concentration and phenology
Phytoplankton growing periods are influenced by global, regional and local driving-factors and they are likely to be altered by climate change. The Chl and phenological trends presented here for the Mediterranean Sea may lead to significant variations in the functioning of the marine ecosystem as result of the alteration of food web structures. The observed trends in phytoplankton Chl concentration and phenology present notable variations at basin and regional scale. Perhaps the most striking observation regarding the trends in phenology is the delay in the timings of initiation b i and peak b t (Fig. 8a and c respectively). These trends do not follow global general trends described for a warming ocean, i.e. earlier b i and b t and b a decline (Hays et al., 2005;Poloczanska et al., 2013). In fact, our analysis reveals a delay in the b i and b t , mainly in oceanic waters. Consistently with our results, a delay in the timing of the phytoplankton growing period in the Mediterranean Sea is coherent with predicted changes in phytoplankton phenology in mid-latitudes reported by Henson et al. (2018) using an ecosystem model-based analysis. A second major pattern observed in the phenological trends is the earlier growing period termination b e (Fig. 8b) and a reduction of duration b d in the eastern basin (Fig. 8d). These phenological changes are coupled with a weakening in the winter growing period amplitude b a (Fig. 8e), which can further reduce the limited productivity of these ultraoligotrophic waters and increase the period in which picoplankton dominates the algal community. These structural changes will have major consequences on the transfer efficiency to higher trophic levels. In the Mediterranean Sea, most fish species are late springearly summer spawners, indicating that their spawning seasonality is coupled to plankton dynamics: fisheries stocks spread their larvae when food availability arises just after the phytoplankton bloom (Álvarez et al., 2012;Tsikliras et al., 2010). The hypothesis is that if the timing of the onset of the phytoplankton growing period advances, fisheries spawning seasonality may also initiate earlier or may become out of phase (mismatch, such as shown in Koeller et al., 2009). Furthermore if the phytoplankton growing period becomes shorter in the eastern basin, the fish assemblage could be less abundant or diverse, causing possible changes in the ecosystem structure and functioning. While an earlier termination b e of the growing period in the eastern basin may be explained by the limited phytoplankton growth (low Chl concentration) caused by scarce nutrient availability in this region, delays in the initiation b i of the main growing period are more difficult to explain on the basis of mixing and subsequent stratification and, more complex processes may be need to examined. Trophic interactions may play an important role modulating the seasonal patterns of Chl concentration in some oceanic regions (Gomez et al., 2017). The delays observed in the timing of initiation b i (Fig. 7a) might be explained by the grazing pressure exerted by winter zooplankton, as reported in other systems (Smetacek and Cloern, 2008). The relative importance of bottom-up and top-down controls at each biogeographical area may explain the different responses observed in the phytoplankton phenological indices. In fact, changes in timing and intensity are not uniform over the globe owing to oceanographically distinct regions and latitudinal drivers (Friedland et al., 2018).
The trends in Chl peak b a observed in the eastern basin are consistent with studies reporting ocean regions, which are becoming more oligotrophic (e.g. Gregg and Rousseaux, 2014). However, globally, this trend has been recently questioned in the light of new and revised data showing no significant trend in global annual median Chl (Gregg et al., 2017). The trend in b a in the western basin is positive. This particularly apparent in the NWMS an in the northern Adriatic (Fig. 8e). The dynamics of phytoplankton in these NWMS is mainly driven by the deep convective formation of Western Mediterranean Water and influenced by the North Atlantic Oscillation NAO (Macias et al., 2017;Mayot et al., 2017)). Based on ecosystem model analysis, Macias et al., 2018 predicted that in the near future (~2030) the spring bloom may occur later in the year and last longer. In the Adriatic Sea, convection is an important mechanism controlling the spring phytoplankton growing period (Gačić et al., 2002); however, the spatial distribution of the positive trend seems to follow the area of influence of Po River. Both winter mixing and river flow in the Mediterranean are closely related to the NAO (Brandimarte et al., 2011;Rixen et al., 2005;Struglia et al., 2004). Patterns of decadal variability driven by NAO could be influencing the trends in the amplitude of the Chl peak b a in these regions (e.g. Basterretxea et al., 2018).
Phytoplankton functional group or species-specific responses to changing environmental conditions are likely to differ from total phytoplankton Chl responses, offering insights into mechanisms of trophic interactions, which are not possible with measurements and analyses of bulk Chl alone (Agirbas et al., 2015;Cabré et al., 2016;Ji et al., 2010a;Racault et al., 2014a). The changes in Chl concentration and species composition have complex implications in the dynamics and diversity of higher trophic levels that the scientific community is only in the beginning to anticipate. Likewise, mismatch between production cycles and spawning season influences the recruitment success of some fish and crustaceans species (Asch, 2015;IOCCG, 2009;Koeller et al., 2009;Platt et al., 2003). Phenological changes can propagate to the rest of the food web and alter the functioning of the marine ecosystem with important consequences for commercially important species (Druon et al., 2015;Kassi et al., 2018).

Conclusions
We have analysed the regional patterns of seasonal and non-seasonal variations in the Mediterranean Sea using ocean-colour remotesensing data. Our results show that seasonal variations are dominant in oceanic areas whereas that non-seasonal variations play an important role in some specific areas of the Mediterranean Sea. In particular, irregular variations may be found in areas of intense wind forcing and/or well strong stablished current flows, and large interannual variations are mainly found in straits and water-exchange zones.
The timings of initiation, peak and termination of the main phytoplankton growing period show differences at regional scale, and between open ocean and shelf waters. The duration of the spring growing period in the western basin (~170 days) is longer than in the more oligotrophic eastern basin (~150 days). The duration generally shortens towards the coast where Chl is enhanced and where the probability of occurrence of a second growing period is higher (i.e. in fall).
Using the longest, continuous, bias-corrected ocean-colour record, reaching nearly two decades of observations (17 years), we observed significant trends in the phenological indices. The duration of phytoplankton growing period has generally increased (+10 ± 7 days decade −1 on average during the period of study 1998-2014) in the western Mediterranean basin whereas an average decrease of −12 ± 7 days decade −1 is observed in the eastern basin. Exceptions to these trends are found in some coastal areas where more complex processes may control phytoplankton phenology. The observed variations in the phytoplankton phenological metrics will be relevant to investigate further and may help us advance understanding of the regulation of the energy flow in the Mediterranean Sea trophic chain.