On the robustness of an eastern boundary upwelling ecosystem exposed to multiple stressors

The resistance of an east border upwelling system was investigated using relative index of marine pelagic biomass estimates under a changing environment spanning 20-years in the strongly exploited southern Canary Current Large marine Ecosystem (sCCLME). We divided the sCCLME in two parts (north and south of Cap Blanc), based on oceanographic regimes. We delineated two size-based groups (“plankton” and “pelagic fish”) corresponding to lower and higher trophic levels, respectively. Over the 20-year period, all spatial remote sensing environmental variables increased significantly, except in the area south of Cap Blanc where sea surface Chlorophyll-a concentrations declined and the upwelling favorable wind was stable. Relative index of marine pelagic abundance was higher in the south area compared to the north area of Cap Blanc. No significant latitudinal shift to the mass center was detected, regardless of trophic level. Relative pelagic abundance did not change, suggesting sCCLME pelagic organisms were able to adapt to changing environmental conditions. Despite strong annual variability and the presence of major stressors (overfishing, climate change), the marine pelagic ressources, mainly fish and plankton remained relatively stable over the two decades, advancing our understanding on the resistance of this east border upwelling system.

The Canary Current Large Marine Ecosystem (CCLME) is the largest eastern boundary upwelling ecosystem (EBUE) globally. It is ranked second to third in terms of primary productivity globally 1 , and is of major economic and social importance in providing sustainable livelihoods, fish-protein supplies, and revenue for the coastal populations and states of West African countries. Because upwelling-favorable winds (UW) have different seasonal dynamics along the coast, important environmental differences exist in different parts of the CCLME 2 . The upwelling is permanent in the central part of the system (∼ 21-26°N), and is less intense and more variable off Morocco (26-33°N), peaking in summer. Upwelling is seasonal from Guinea to Senegal (10-16°N). Finally, it gradually changes from winter to spring in the productive season between Mauritania and Senegal (16-21°N) 3 . The seasonality of the upwelling in this region is associated with seasonal movements of the intertropical convergence zone (ITCZ), which drives latitudinal changes in the trade winds 4,5 . The wind-forced upwelling of deep, nutrient-rich water is responsible for high phytoplankton primary productivity in coastal waters 6 , which provides food for higher trophic levels.
Marine organisms, including phytoplankton to small pelagic fishes, are sensitive to environmental changes 7 . For many fish species, the interaction between fishing and climate impacts the life history parameters of populations (growth, maturation, recruitment), migration, spatial distribution, and food web complexity and stability 8 . The fisheries of certain pelagic fish species, such as the Atlantic herring (Clupea harengus) and sardine (Sardina pilchardus), have fluctuated for many hundreds years, frequently exhibiting alternating population dominance 9,10 . Previous studies conducted off the southern part of the CCLME 'sCCLME' (Sahara Bank to South Senegal) highlighted changes to the abundance and distribution limits of a number of fish species in relation to environmental changes, such as cooling or warming waters 11,12 . Long-term changes to zooplankton and micronekton abundance and distribution have, however, received less focus in the CCLME. Yet, these components of the ecosystem are highly sensitive to environmental changes (e.g., nutrient concentration, salinity, oxygen and temperature) 13,14 .
Because of the short life cycles of these species, intermediate communities (zooplankton and micronekton) respond rapidly to mesoscale physical processes associated with currents and frontal structures linked to coastal Figure 1. Study area was limited to the continental shelf of the south part of the Canary Current Large Marine Ecosystem. The dashed line represent the Cap Blanc limit dividing the study area into two part. We conductedpelagic acoustic sea survey on the same vessel and used the same protocol to process the total acoustic density of the Nautical Area Scattering Coefficient (NASC, m 2 nmi −2 ) from 1995 to 2015, as a proxy of pelagic abundance. The map was performed with Matlab R 2018a 20 (https ://fr.mathw orks.com/produ cts/new_produ cts/ relea se201 8a.html).    www.nature.com/scientificreports/ the two areas. Positive anomalies were higher in the north, reaching 5 m s −1 . In comparison, negative anomalies were around 2 m s −1 in the north and south of Cap Blanc (Fig. 3). Sea surface chlorophyll-a concentration (SSC) showed the opposite pattern, with a significant increase in the north sCCLME and a very significant decrease in the south (Fig. 2, Table 1). Low SSC positive (+ 2 mg m 3 ) and negative (< − 1 mg m 3 ) anomalies were recorded in the north of Cap Blanc. In comparison, in the south, high positive (+ 9 mg m 3 ) and negative (− 4 mg m 3 ) anomalies were recorded. Barycenter displacement analysis. No significant trends were detected for the barycenter of the acoustic density (expressed in Nautical Area Scattering Coefficient (NASC or s A , m 2 nmi −2 )), which was used as a proxy of pelagic abundance 'PA' , for both groups in the acoustic surveys (i.e., plankton group (PG) and pelagic fish group (PFG)) from 1995 to 2015 (Fig. 4). The barycenter of PA showed strong inter-annual fluctuations, particularly in the area south of Cap Blanc where the position of the PG and PFG barycenters was highly variables, sometimes overlapping. There was no significant displacement for both groups in both surveyed areas (north and south of Cap Blanc).
Analysis of pelagic abundance and anomaly trends. PA were higher for PFG compared to PG, with a quantitative ratio of 1/10. PA was higher in the south compared to the north of Cap Blanc for both groups (Fig. 5). No significant trend was observed in PA during the study period ( Table 2, p-value > 0.05) in the sCCLME. PA was particularly high in 2015 in the area south of Cap Blanc for PFG, whereas abundance dropped by 50% in the north (Fig. 5, Fig. S2). PG showed more stability over the long-term. High positive PA anomaly (> 2 × 10 7 m 2 nmi −2 ) was observed for PFG in 2002 in the area north of Cap Blanc. High negative anomalies . For PFG, PA was also higher at night in the surface waters of both areas. The day-night difference in PA was significant (p-value < 0.000) for all groups (Table 3). Δ DVM (i.e., night PA-day PA) mainly produced  Table 2. Spearman trend tests for the barycenter displacement of pelagic abundance and Δ DVM trend (i.e., differences between day and night pelagic abundance for the plankton group (PG) and pelagic fish group (PFG)) in the areas north and south of Cap Blanc. r: Spearman correlation coefficient is shown in the first row of each cell and the p-value is shown in parentheses. www.nature.com/scientificreports/ values, with higher NASC densities at night compared to by day (Fig. 7). Δ DVM was higher in the south compared to the north of Cap Blanc for both acoustic groups, and mostly for PG. The DVM trend was also evaluated over the 20 years using Δ DVM , and no significant trends were found ( Table 2). Significant DVM day and night differences were reported in both areas for all groups (Table 3).

Discussion
In recent years, many changes or trends have been documented in physical forcing at regional and global scales 21 . Most large marine ecosystems (LME) were subjected to continuous warming between 1982 and 2006 22 . A mean warming trend was detected over the last three decades in the whole CCLME region, with substantial warming since 1995 12,23 . Previous study 24 detected significant SST (p < 0.001) warming trends in the CCLME from 1993 to 2014, with markedly higher SST in the southern part (around 0.4 °C decade −1 ) compared to the northern part (about 0.2 °C decade −1 ). Our results support these previous findings, whereby SST warming in the sCCLME was stronger in the area south of Cap Blanc. Bakun 25 hypothesized that, in the context of the climate change, global warming will enhance land-sea temperature gradients, enhancing UW. We found that UW significantly increased  www.nature.com/scientificreports/ in the area north of Cap Blanc (5.5 m s −1 decade −1 ), whereas no UW trend was recorded in the south. Benazzouz et al. 26 documented the same pattern, with UW strongly increasing between 26°N and 33°N from 1982 to 2011, with this trend being much weaker between 21°N and 26°N. Previous studies demonstrated that warming oceans, including in the CCLME, are associated with a decline in average chlorophyll-a concentration 'Chl-a' , which is used as proxy of phytoplankton biomass and primary productivity [27][28][29] . This declining trend might be explained by a physically mediated effect of upper-ocean warming on vertical stratification, which indirectly affects phytoplankton by limiting nutrient supply to the sunlit layer 28 . It might also be explained by the positive effect of warming on the metabolic rates of plankton 30 . Previous observations in the CCLME indicated that Chl-a declined during 1998-2007 31,32 . In other areas, there is evidence that increasing wind speeds overcome any potential increase in stratification, due to warming, increasing Chl-a production 33 . In our study, SSC, which was used as a proxy for Chl-a, increased significantly in the area north of Cap Blanc, and decreased in the south. Spatial variation in UW and SSC was consistent in the north of Cap Blanc, where increasing SSC coincided with increasing UW. However, in the area south of Cap Blanc, SSC and UW trends were not spatially correlated. Indeed, coastal upwellings depend on a combination of wind intensity, coastal topographic effects, and alongshore geostrophic flow (related to large-scale circulation patterns) 2 . Moreover, there is a non-monotonic relationship between UW and biological responses, where strong UW limits productivity and phytoplankton levels 34,35 . Moderate wind speed provides optimal conditions for the development of coastal phytoplankton populations 36,37 . The PA is 10 times greater for PFG compared to PG, supporting existing studies suggesting that there is a factor of 10 between fish and zooplankton biomass in the marine food web 38 . Marine ecosystems have generally inverted biomass pyramids 39,40 with the highest biomasses in the highest trophic levels i.e. the producers are small organisms with least biomass and the biomass gradually increase towards the top of the pyramid. It is due to higher reproducibility and shorter lifespan of phytoplanktons, that although their biomass is less at any time, they frequently replenish themselves to meet the increased needs of zooplankton and larger fish. The increase of prey growth rate, the conversion efficiency, the prey carrying capacity, or the predator life span robustly facilitates the development of inverted biomass pyramids 41 . In the sCCLME, no significant trend was observed for the PA barycenter for all trophic groups. This absence of spatial shift in PA might be related to the high productivity of the sCCLME 2 . Future advances in remote technologies are expected to facilitate the greater discrimination of pelagic organisms 42 .
The geographic range and abundance of pelagic species are strongly associated with the temperature tolerance of organisms (thermal window) 43,44 . Temperature is a key driver of changes to communities at different spatial scales. Increases in temperature lead to the dominance of species adapted to warmer waters in the pelagic community 45 . Out of the most important species in the sCCLME, two categories emerged based on their temperature preferences: (1) cold water species, such as Sardina pilchardus, Scomber colias, and Trachurus trachurus, and (2) warm water fishes, such as the Sardinella aurita, Sardinella maderensis, and Trachurus trecae 46 . The abundance of S. aurita, a dominant pelagic fish species in the sCCLME, responded to warming. Low catches of S. aurita were associated to temperatures below 21 °C in the waters off Mauritania, whereas high catches were associated with temperatures above 21 °C in the waters off Senegal 12,47,48 . Range shifts in marine species linked to SST warming have been documented across all ocean regions 7 . Sarré et al. 49 reported northward shifts in the distribution of sardinella and other small pelagic species in the sCCLME, which contrasted with our results, despite falling within the same study period. The authors attributed the northward shift of these low thermal tolerance species to high warming trends in the southern part of sCCLME. The population dynamics and physiological rates of zooplankton are strongly linked to temperature 14,50 . In the north-east Atlantic, the distribution of two species of copepods (Centropages chierchiae and Temora stylifera) shifted northwards in response to rising sea temperature 51,52 . This northward shift of copepod assemblages in response to global warming was estimated to be around 260 km per decade 50 .
The PA of PG and PFG was higher in the area to the south versus north of Cap Blanc, which was attributed to oceanographic conditions, particularly Chl-a 53 . The area south of Cap Blanc was characterized by significant SST warming associated with a stable UW, which provided a favorable habitat for pelagic species 7 . Moreover, previous studies 54 showed that marine organisms aggregate as dense layers of PA under stable physical conditions. In our case, the area north of the Cap Blanc was characterized by high levels of UW, which did not favor the presence of a dense PA layer.
Although, previous studies 49 in North West Africa demonstrated a northward shift of S. aurita, which is a key PFG species, other pelagic fish species exhibited phenotypical adaptations to changes in temperature 55 . Many fish species have geographic ranges spanning a wide temperature gradient, suggesting a capacity for acclimation or adaptation to changing temperature 56 . Indeed, temperature is a limiting factor to fish species with tropical affinity, but not to species with wider temperature tolerance. For example, S. maderensis, which tolerates high changes in temperature 47 , did not exhibit a northward shift 49 . However, oceanographic surveys are typically implemented in November-December during the southward migration of S. aurita from Morocco to Senegal 44 . This phenomenon might partly explain the higher PFG PA detected to the south of Cap Blanc in the current study. The few studies investigating the response of marine fishes to climate change in the sCCLME have focused on fish, and usually single species, rather than pelagic organisms as a whole. Yet, the population dynamics of zooplankton are tightly linked to temperature, which exhibit rapid responses to changing environments. Zooplankton species also respond to SST warming by shifts in their distribution 51,52 and phenology 57 , or by changing their morphological traits, such as reducing body size 58,59 .
Relative spatial coherence was documented between environmental and PA anomalies (Fig. S5)  www.nature.com/scientificreports/ in the area south of Cap Blanc in 2015 (paralleled by a 50% drop in the north in the same year) could not be attributed to any environmental change, as no environmental anomaly was observed. The adaptive DVM behavior 60,61 has been reported, and was higher for PG, which had high Δ DVM . Small pelagic fishes disperse widely at night and aggregate by day with small bathymetric changes 62 . In comparison, planktonic organism performs true DVM; specifically, they occupy deep waters during the day and migrate toward the surface at night to feed 63 . Such behaviors explain the higher extent of DVM (i.e., high Δ DVM ) reported for PG. In the area north of Cap Blanc, we documented higher PA during day compared to at night for PG in the upper water column (20-40 m). This phenomenon indicates an inverse DVM behavior, whereby plankton ascend in at dawn and descend at dusk. This inverse DVM pattern, named DVM type II, was previously documented for planktonic fish larvae in the sCCLME 63 . The more marked DVM (higher positive Δ DVM ) in the south of Cap Blanc suggests that this area is mostly populated by communities performing DVM, particularly PG. Yet, the absence of a clear DVM trend allow to confirm that key functional pelagic groups remain stable in system over the study period.
Resistance is the ability for an ecosystem to remain unchanged when being subjected to a disturbance or disturbances. The absence of significant long-term trends in PA for both trophic level groups (PG and PFG), and absence of longitudinal shifts, demonstrates the strong resistanceof the sCCLME ecosystem to environmental stressors (overfishing, climate change). Indeed, despite the well-known overfishing statute of the main pelagic fish stocks and environmental change documented in the sCCLME, particularly for SST and UW shifts, no trend emerged for long-term variation in PG or PFG PA and DVM behavior. The sCCLME, particularly in the area south of Cap Blanc, presented strong inter-annual environmental variability, which did not impact long-term variation in PG and PFG. Moreover, as a tropical system, the sCCLME is characterized by high species diversity 64 , which increases system resistance. Diversity, connectivity, and adaptive capacity represent ecological properties that underlie ecosystem resilience 65 . Biological diversity in sCCLME might contribute towards resistance by increasing the likelihood that some species and/or functional groups are resistant to stressors. This phenomenon allows species to compensate for one another within a community, and facilitates ecological processes vital for recovery and adaptation. The phenotypic plasticity of key species in the PFG, such as S. aurita 55,66 and S. maderensis 47 , to the changing climate facilitates the rapid evolution of traits better suited to new conditions. This work presents a first step in understanding the resistance of the sCCLME at the pelagic level. Future studies are encouraged to investigate the ecological responses of the sCCLME to the effects of climate change at individual, population, and community levels.

Methods
Survey area. The research vessel R/V Dr. Fridtjof Nansen (DFN) conducted regular surveys in the marine area off West Africa. We focused on data collected from the southern part of the CCLME (hereafter named sCCLME; Fig. 1). Cap Blanc (20.77°N) is considered to represent a "faunistic limit" for the planktonic population 63 , and the northernmost distributional area for several small pelagic fishes in North-West Africa 67 . Thus, we divided our surveys in the sCCLME into two parts: area north of Cap Blanc vs. area south of Cap Blanc, extending from Cap Cantin (31° 65′ N) in Morocco to Casamance (12° 31′ N) in southern Senegal, respectively.
The sCCLME is rich in zooplankton, with copepods dominating. Copepod species account for 60 to 95% of total zooplankton abundance in the sCCLME, and constitute the bulk of PA of macro and mesozooplankton 68,69 . At a higher trophic level, the bulk of pelagic biomass is represented by small pelagic fishes (e.g., S. aurita, S. maderensis), which are also key for safeguarding food security in West Africa.
Acoustics and environmental data processing. Acoustic data. Fisheries acoustic surveys were performed annually in November or December (upwelling season) from 1995 to 2004, and again in 2006, 2011, and 2015 ( Fig. 1), on-board the research vessel R/V Dr. Fridtjof Nansen (DFN). The way in which total acoustic density is distributed divides the study area into two parts: north and south, separated by Cap Blanc. These surveys were conducted using a 38 kHz echosounder at depths ranging from 10 to 500 m isobaths. We deployed an ES38-B transceiver, towed at a depth of 5.5 m, with the following settings: absorption coefficient of 8.7 dB km −1 , pulse length of 1.024 ms −1 , and maximum used transmission power of 2000 W 70 . The echo sounder was calibrated following classic calibration procedures 71 . Acoustic data were processed with an offset of 10 m below the sea surface to avoid integrating air bubbles during post processing.
Environmental data. To obtain a homogenized overview of the area, environmental data were collected from satellite remote sensing products: upwelling wind (UW), sea surface chlorophyll-a concentration (SSC), and sea surface temperature (SST). UW is the upwelling-favorable meridian wind speed at the sea surface, and is used as an upwelling index 3 . UW was extracted from the daily CCMP (Cross-Calibrated MultiPlatform) wind product V2.0 at a 0.25° spatial resolution (available at http://www.remss .com/measu remen ts/ccmp/) from 1988 to 2017. SSC is an adequate proxy of primary productivity models on which biomass is directly based 32 . SSC was collected daily from the SeaWIFS (1997)(1998)(1999)(2000)(2001)(2002) and AQUA-MODIS sensor (2003-2018) (available at https ://ocean color .gsfc.nasa.gov/). SST data were extracted from daily day-time series of the pathfinder AVHRR dataset version 5.2, from 1982 to 2018, at 4 km resolution (available at https ://www.nodc.noaa.gov/Satel liteD ata/pathfi nder 4km53 /). The spatial averages of satellite data were computed using a fixed distance from the coast of 100 km. The monthly average of all data (UW, SSC, and SST) were calculated, and the spatial averages were computed for the sCCLME area. Annual climatological anomalies were estimated to elucidate fluctuations of the different environmental variables and support the trend analysis. www.nature.com/scientificreports/ Data analysis. Acoustic data were echo-integrated per 0.1 nautical miles (nmi) using 'Matecho' 72 . To obtain homogenous data, only continental shelf data were considered (i.e., 10-150 m). The relative acoustic density, expressed as Nautical Area Scattering Coefficient (NASC or s A , m 2 nmi −2 ) 73 , was used as a proxy of marine organism abundance to assess spatial variability. In the sCCLME area, the abundance of small pelagic fish is closely correlated to the upwelling of nutrient-rich waters, which is the basis of food production for these species 2,48 . Given that small pelagic fish mainly feed on phyto-and zooplankton, which are very sensitive to environmental changes, changes to primary production might affect the abundance of stocks and, thus, both the structure and functioning of pelagic ecosystems. Consequently, a "thresholding" approach was applied, using volume backscattering strength (S v , expressed in dB), to separate marine organisms into three echo classes: "low acoustic trophic level" [− 80 ≤ S v < − 65 dB] 74,75 , "high acoustic trophic level" [− 65 ≤ S v < − 20 dB] 76-78 and the "all trophic level" [− 80 ≤ S v < − 20 dB]. Low and high acoustic trophic level classes were assumed to represent plankton (i.e., meso and macrozooplankton) vs. pelagic fish, respectively. Usually, small pelagic fishes occurring in the CCLME have a S v , higher than ~ − 68 dB 79 , while macrozooplankton backscatter at < − 70 dB 16,80 . Low, high, and total acoustic trophic levels were, thus, designated as plankton (PG), pelagic fish group (PFG), and total group (TG), respectively.
To detect spatial changes in the spatial distribution of pelagic marine organisms in the areas north and south of Cap Blanc, we analyzed the barycenter displacements of PA over the time-series. The barycenter is the average position of acoustic density for 1 year, weighted by the biomass of each zone. For instance, if "B ai " represents the biomass in zone "Z i " for year "a" and "A i " the coordinates of the center of "Z i ", the coordinates of the barycenter "C a " of the biomass would fit the equation below as: And the coordinates in latitude "X a " and longitude "Y a " of C a are: All statistical analyses were completed in a free software environment for statistical computing and graphics 81 . Figures and maps were also performed with R and Matlab R2018. Spatial and temporal variability of the acoustic abundance time-series were investigated for: "zooplankton, " "pelagic fish, " and "total. " The statistical significance of the annual position (latitude and longitude) of changes to the barycenter and PA trend were tested using Spearman tests, which are considered valid for testing the significance of trends 82 . The Spearman rank correlation can be used to determine if a variable is increasing or decreasing with time. Spearman's rho is the product-moment correlation between the ranks of paired data. To test for trend, one member of the pair is the time of observation, the other is the observed variable 83 Trends in the annual means of environmental variables (i.e., SST, UW, and SSC) were assessed in the north and south of Cap Blanc using Spearman tests, which are valid for supplementing linear models with a non-parametric models, such as the Spearman correlation. Diel changes to PA in the sCCLME were first analyzed by comparing PA (log 10 transformed) for day vs night in the vertical dimension (i.e., the mean value for each 1-m depth step over the whole period study). Then, Δ DVM (night PA-day PA) was computed for each 1-m depth step and each year. Diel change in the annual PA for each acoustic group was tested using the Wilcoxon test. The Spearman test was used to assess DVM trends for PG and PFG in the north and south of Cap Blanc. Diel transition periods were removed from the analyses to prevent diel vertical migrations causing bias to changes in acoustic density. Transition periods were defined using sun altitude, with sunset and sunrise corresponding to a sun altitude of 18° and + 18°, respectively. Time series anomalies for environmental variables and acoustic groups were calculated by computing differences between actual (annual or monthly) and long-term average values.

Data availability
The public cannot access our acoustic data because they belong to the partners who funded the oceanographic cruise. Environmental data are already available as indicated in the methodology part.