Asynchronized spawning responses of small pelagic fishes to a short-term environmental change

: We provide substantial evidence on how short-term changes in environmental conditions activate and deactivate spawning activities in small pelagic fishes. An ichthyoplankton survey was conducted along the southernmost part of the Canary Current upwelling ecosystem in May 2013, covering the area twice within 20 d. This period coincided with a strong environmental change from a cold productive upwelling regime to a warmer and less productive upwelling relaxation event. This change triggered a shift in spawning activity from European anchovy Engraulis encrasicolus to round/flat sardinella Sardinella spp. We used zero-altered negative bi - nomial regression models with a generalized additive structure based on integrated nested Laplace approximations to link early larval distribution patterns to the 2 different regimes. The models confirmed 2 species-specific temperature spawning windows, suggesting a spawning pause of anchovy during upwelling relaxation while simultaneously activating spawning in sardinella. Observing immediate spawning responses to the 2 environmental regimes underlines the as - sumption that windows of spawning opportunity are the main drivers of small pelagic fish fluctuations in upwelling regions. The duration of a specific environmental condition can, therefore, increase or decrease the chances for reproductive success. The observations of this study may explain why certain small pelagic fish species can dominate over others during a particular period and might also apply to other upwelling regions of the world oceans where upwelling and relaxation events alternate.


INTRODUCTION
In upwelling ecosystems, fish communities are dominated by small pelagic fishes such as anchovies, sardinellas, sardines, mackerels and horse mackerels (Cury et al. 2000). To coexist and to maintain persistence, these species have evolved specific life-traits. One of these is to safeguard reproductive success, which is dependent on when and where the adults release their offspring. In small pelagic fishes such as anchovy and sardinellas, the timing of spawning is critical, as it is linked to the environmental conditions that determine the fate of the offspring and, consequently, successful recruitment (Hjort 1914, Cushing 1969, Houde 2008. Fish species are adapted to upwelling processes in eastern boundary ecosystems (Bakun 1996) and may have evolved spawning strategies around the onset and relaxation of upwelling conditions (Cury & Roy 1989. Here, upwelling processes result from a deep nutrient-rich, cold-water uplift, boosting plankton productivity that propagates through the food web (Hutchings et al. 1995). Small pelagic fishes take advantage of these rich planktonic food resources during and after an upwelling event, which is often linked to their spawning initialization (Roy et al. 1992, Mbaye et al. 2015. Populations of different small pelagic fishes coexist in large ecosystems in which a given physical state can be dominant within seasonal to decadal timeframes (Schwartzlose et al. 1999, Sánchez-Garrido et al. 2019), favouring one species over the others. In the California Current upwelling ecosystem, for instance, Pacific sardine Sardinops sagax and northern anchovy Engraulis mordax peak spawn at 15°C (Lluch-Belda et al. 1991). Data on later stages of larvae, however, show that only Pacific sardine larvae seem to thrive when the ecosystem changes into a warmer state, while northern anchovy larvae only maintain high abundances under colder conditions. This can result in changing dominance regimes between Pacific sardine and northern anchovy in the long-term driven by multidecadal climate fluctuations (Lindegren et al. 2013). In contrast, modelling studies of European pilchard Sardina pilchardus and European anchovy E. encrasicolus from the Canary Current upwelling ecosystem suggest that the stock size asynchrony differs from that observed in the California Current upwelling ecosystem (Sánchez-Garrido et al. 2021). Instead of temperature favouring one stock over the other, food availability and its utilization seem to be the main factors modulating their recruitment success. All these multidecadal observations, however, do not explain the more detailed processes occurring on shorter timescales, and there is a lack of knowledge on how the spawning responses of small pelagic species are adapted to fastchanging environments and how resources are utilized for spawning. Hence, field-based observations from dynamic upwelling ecosystems could prove helpful in identifying asynchronized traits in different small pelagic species over shorter timescales.
The Canary Current upwelling ecosystem is partitioned into a weak permanent (26−35°N), a permanent (21−26°N) and a seasonal upwelling area (Cropper et al. 2014). The permanent upwelling areas are dominated by the more temperate-associated species European anchovy and European pilchard, while during the non-upwelling phases, especially in the seasonal upwelling area in the south, sardinellas (round sardinella Sardinella aurita and flat sardinella S. maderensis combined as sardinella) migrate northward and become the most dominant small pelagic fish species (Braham & Corten 2015). The seasonal upwelling region is particularly dynamic along the coastal area off Senegal, The Gambia and Guinea-Bissau (Roy 1998, Capet et al. 2017. It is characterized by upwelling during winter and early spring, while stratification starts to increase with a ceasing of upwelling periods in late spring (Marcello et al. 2011). At this time, small pelagic fishes such as anchovy and sardinellas occur frequently in the region (Zeeberg et al. 2008, Diankha et al. 2015, Ba et al. 2016. Migration patterns of both species are suggested to be driven by changes in the upwelling and stratification regimes. While anchovies normally prefer periods of intense upwelling, sardinellas prevail in areas dominated by seasonally altering upwelling and stratified periods (Binet et al. 1998. Migration patterns of these species are assumed to be driven by upwelling and relaxation conditions (Zeeberg et al. 2008, Braham & Corten 2015, which seasonally shift from southern Senegal to northern Mauritania (Mittelstaedt 1983). The upwelling expands southward from October to May and relaxes northwards from April to September (Arístegui et al. 2009).
In this paper, we investigate the influence of shortterm changes in environmental conditions on the spawning activities of small pelagic fishes at the southernmost boundary of the Canary Current upwelling ecosystem, covering the upwelling area off Senegal, The Gambia and Guinea-Bissau. The study area is an important spawning and nursery ground for small pelagic fishes , Ndour et al. 2018) and serves as an ideal environment to investigate rapid ecological responses to environmental changes. The upwelling process in the region can be observed as a cold-water tongue between October and May (Marcello et al. 2011) stretching from the Cap-Vert peninsula 300−400 km southward (Roy 1998, Capet et al. 2017. As the southernmost stretch of the upwelling plume of the Canary Current upwelling system, this region is also known for variations in upwelling intensity (Ndoye et al. 2014), with alternations between upwelling and relaxation events occurring on a 3−10 d basis (Ndoye et al. 2014).
We carried out an oceanographic survey during the expected transition period from an upwelling to a stratified regime, covering the area twice within 20 d. We conducted pelagic fish and ichthyoplankton sam-pling to explore the alternation of spawning responses of the most abundant small pelagic fishes -European anchovy E. encrasicolus and the 2 most common sardinella species (S. aurita and S. maderensis). Off Senegal, May is the month of transition from the upwelling to the stratified regime and was therefore considered a suitable time to in vestigate the rapid and short-term spawning re sponses of fishes to changes in the environment. Our main question in this paper was how an up welling (cold) and subsequent relaxation event (warm) are utilized for spawning, which may also provide a deeper understanding of small pelagic fish stock alternations occurring worldwide.

Sampling design
A 20 d survey was conducted in May 2013 onboard the R/V 'Dr. Fridtjof Nansen'. Sampling took place within the upwelling period (Marcello et al. 2011) and consisted of 2 consequent coverages of the study area (Leg 1: 3−12 May; Leg 2: 13−22 May). Both legs started from the south and progressed northwards. Leg 1 covered the area between central Guinea-Bissau and north Senegal at about 15.7° N (Fig. 1), while leg 2 extended between southern Guinea-Bissau and just north of Cap-Vert. Concurrent receptions of daily remote sensing imagery (MODIS-Aqua; https://oceancolor.gsfc.nasa.gov) for monitoring environmental conditions revealed that Leg 1 coincided with the period of upwelling and Leg 2 with a subsequent relaxation event. Thus, the ichthyoplankton sampling coverage of each leg was adjusted to the 2 environmental conditions recorded.

Ichthyoplankton sampling and identification
The ichthyoplankton sampling grid comprised 128 stations (Leg 1, n = 62; Leg 2, n = 66), located along transects 15 nautical miles apart, with an interval of approximately 10 nautical miles between stations. Due to the vessel's draft, diel sampling started at about 20 m bottom depth, continuing offshore perpendicularly to the coast and ending at about 1000 m bottom depth. At each station, a Multinet Midi with a square 0.25 m 2 mouth-opening area and mesh size of 405 μm (type: 438 130; Hydro-Bios, www.hydrobios. de) was deployed. The Multinet was lowered down to either a maximum depth of 200 m or to 5−10 m above the seafloor at stations with bottom depths < 200 m.
The Multinet was then towed obliquely to the surface, using 1 net per sampling occasion, thereby collecting ichthyoplankton in the cod-end that was fitted with a 180 μm mesh. The Multinet was deployed at 3 selected transects (T1 in Leg 1; T2 and T3 in Leg 2), using up to 5 nets per sampling event to enable depth-stratified diel sampling (at 200−100, 100−75, 75−50, 50−25 and 25−0 m). This data provided additional insights regarding the depth niche of the fish larvae and the environmental conditions associated with their distribution.
At all stations, when the Multinet had reached its maximum lower sampling depth, the first net was opened and the Multinet was retrieved with an average (± SD) towing speed of ~1.6 ± 0.4 m s −1 . After each haul, the net and cod-end contents were gently washed down, and fish larvae were immediately sorted from the fresh sample. Anchovy and sardinella eggs and larvae were identified based on Russell (1976) and Olivar & Fortuño (1991), respectively. The difficulty in discriminating between Sardinella species in the area (i.e. S. aurita, S. maderensis) led us to group all eggs and larvae as Sardinella spp. Total length of larvae was measured before fixation. After sorting for fish larvae, eggs were counted from a subset of the whole sample using a Motoda plankton splitter (Motoda 1959). Larvae were preserved in 96% ethanol and/or 4% borax buffered formaldehyde solution.

Hydrological sampling
Vertical temperature and salinity profiles were obtained using a Seabird 911+ CTD with rosettemounted water collectors at each of the Multinet stations. Salinity measurements were validated using a Portasal salinometer (model 8410A), analysing the salinity of the deepest Niskin bottle sample from each ichthyoplankton station. Dissolved oxygen was measured with a CTD-mounted oxygen sensor. Oxygen was validated using the Winkler redox titration method (Winkler 1888) from 10 Niskin bottle samples comprising depths between 5 and 2000 m. A CTDattached Chelsea Mk III Aquatracka fluorometer measured in situ fluorescence on a relative scale.

Biological fish sampling -gonad staging
To verify that mature fish were present during the survey and to explore the distribution of adult stocks of the species of interest (E. encrasicolus, S. aurita, S. maderensis), a total of 15 fishing trawls were conducted (1−8 in Leg 1; 9−15 in Leg 2). The location of each trawl was based on prior acoustic recordings of pelagic fish in the area. Acoustic recordings were conducted using a calibrated ES38B (38 kHz) transducer connected to a Simrad ER60 Echosounder. The large scale survey system (LSSS) (Korneliussen et al. 2006) was used to post-process the raw acoustic data and allocate to species categories (sardinella, anchovy, etc.) based on the proportional contribution of the species in the trawl catches and their acoustic signature. Data were exported from LSSS as nautical area scattering coefficient (NASC; m 2 nautical mile −2 ). The method is described in more detail in Text S1 in the Supplement at www.int-res.com/articles/suppl/m696 p085_supp.pdf.
After each trawl, total catch weight for each species was recorded. Information on the sex, maturity stages, body and gonad weight were collected from the first randomly selected 20−30 individuals of each target species. A 6-phase maturity scale was used to describe their reproductive state (1: immature; 2: maturing virgin and recovering spent; 3: mature before spawning; 4: mature pre-spawning; 5: mature in spawning; 6: post-spawning). These were further considered as immature (i.e. 1, 2, 6) vs. mature (i.e. 3, 4, 5). The gonadosomatic index (GSI) for each individual was calculated according to: (1)

Remote sensing data analysis and wind
To describe the environmental regimes during the 2 legs, mean sea surface temperature (SST) and mean sea surface chlorophyll a (SSCa) maps were produced using data provided by Copernicus Marine Services (https://marine.copernicus.eu/). The OSTIA global foundation SST product, with daily gap-free maps at 0.05 × 0.05° resolution (Donlon et al. 2012), was used to display SST. SSCa maps were produced using the ocean colour product OC-ACRI-NICE-FR with daily cloud-free 4 × 4 km grids (Colella et al. 2020). The daily data products were merged to calculate their mean values for each leg and plotted using the 'akima' package (Akima 2016) in R version 4.0.2 (R Core Team 2020) for interpolating purposes.
Additional meteorological data for wind was recorded using a weather station equipped with a wind observer sensor on top of the vessel's mast. Hourly average data on wind direction and speed were used to confirm wind conditions that can be related to the upwelling and relaxation phases during the 2 legs.

Ichthyoplankton data analysis and modelling
Egg and larval distributions were mapped along the study area as depth-integrated abundances standardized to the number of individuals in a 10 square meter surface area. Maps were compiled using the R packages 'maps', 'mapdata' (Brownrigg 2016) and 'marmap' (Pante & Simon-Bouhet 2013). To evaluate the spawning responses of anchovy and sardinella, we used a generalized additive modelling approach with integrated nested Laplace approximations (INLA) for approximate Bayesian inference  in R version 4.0.2 using the R package 'R-INLA' . Assuming that the different environmental regimes activate or deactivate spawning, we tested the response of larval fish count data to different environmental co-variables.
As a response variable, we used only larval fish stages, excluding eggs that persist for too short a time in the environment and any post-larval stages that would persist for too long in the environment. For instance, in areas of SST > 20°C, the duration of the egg stage is generally very short (<2 d) (Somarakis et al. 2002). Standard length measurements for both species were made on larvae collected along T1 (Leg 1) and throughout the sampling grid during Leg 2. Thus, these measurements should only be considered as a subset of standard length frequencies.
During the 10 d of each leg, successive spawning events may have taken place at different sites within the surveyed area. As environmental co-variables, we tested average values of temperature, fluorescence, oxygen concentration and salinity from CTD casts within the surface 1−25 m depth-layer, where larvae were most abundant (as revealed from their vertical distribution patterns from the 3 additional vertical stratified samplings; Fig. S1). Additional covariables were bottom depth, latitude, longitude, and filtered water volume during the Multinet sampling.
We formulated separate models for each species to answer the 2 following questions: (1) what were the main environmental variables related to the presence or absence of the larval fish (i.e. positive vs. negative sites), and (2) what were the main environmental variables related to the number of larval counts at the positive sites? We also aimed to identify ranges of temperature that may shape the distribution of anchovy and sardinella larvae and determined whether the ranges differed between the 2 taxa. We also assessed if there is an important nonlinear relationship between temperature and larval count data. This information allowed for the detection of spawning windows for both species and to relate these to upwelling and relaxation conditions. Data exploration, according to Zuur et al. (2010), preceded model formulation. Collinear co-variables were excluded based on Pearson correlation (r) > 0.7 and a variance inflation factor (VIF) > 3 (Table S1). For the sardinella model, we dropped salinity, oxygen, filtered water volume and longitude to achieve a VIF < 3. For anchovy, we used longitude instead of bottom depth (Pearson correlation of both co-variables: r = 0.74; VIF = 4.1).
Our data set was also characterized by spatial and temporal dependency (Carroll & Pearson 2000), which is typical for quasi-synoptic sampled data from ocean expeditions (McManus & Woodson 2012). To address these dependencies in the models, we allowed a separate spatial random field inside R-INLA for each leg (Zuur & Ieno 2018). We constructed a Matérn stochastic partial differential equation model (Lindgren et al. 2011) based on penalized complexity priors (Simpson et al. 2017) with a range (r) and a standard deviation parameter (σ) for each species model. This was implemented as a mesh based on a constrained refined Delaunay triangulation (Fig. S3) and defines the spatial domain of the sampling area. Based on the distances between sampling sites, we assumed a probability of stations that are correlated as: P(r < 30 km) = 0.05 and P(σ > 2) = 0.05 for the sardinella model and P(r < 100 km) = 0.05 and P(σ > 2) = 0.05 for the anchovy model.
To answer our 2 main questions (1 and 2 above), we used a zero-inflated altered negative binomial (ZANB) approach with a generalized additive structure and a spatial random field for each leg inside the presence model (Zuur & Ieno 2018). This model, also known as a hurdle model, combines a Bernoulli model that tests the presence−absence structure of the response variable together with a zero-truncated negative binomial model testing the presence-only part. The zero-inflated approach was used due to the high rate of zeros in the count data (sardinella: 41% zeros; anchovy: 65% zeros).
The ZANB is defined as follows (Zuur & Ieno 2018): ( 3) where (4) log (π i ) = Intercept + covariables i where the expected value Count i is μ i with a negative binomial distribution, π i , a Bernoulli distribution and a parameter k that allows for additional variation. The expected value E (Count i ) is μ i and π i with μ i as a function of a log-link and π i as a function of a logitlink; p 0 is the probability of larval fish presence. The response variable count data is either the number of sardinella or anchovy larvae sampled per station (i).
A negative binomial distribution was used due to the large variation in the presence count data. All co- variables were used as standardized continuous linear terms except temperature, which was used as a smoothing function in the models for both species. Latitude was used as an additional smoother in the anchovy model to address large-scale spatial correlation. v i , j is the random spatial term in the negative binomial model part that is al lowed to change be tween the 2 legs ( j) to address the smaller spatial correlations as well as temporal dependency (Zuur et al. 2017). A cubic regression spline smoother was used in both models for temperature with either 4 knots (knots − 1 = degrees of freedom) for the sardinella model or 5 knots for the an chovy model. The latitude smoother was only applied in the anchovy model with 3 knots. Limiting the degrees of freedoms defines the maximum number of regressions in a smoother curve (see Table 1: Temperature 1-4, Latitude 1-2). It im proves the interpretability of non-linear relationships. We applied model selection using the Deviance Information Criterion (DIC) and the Watanabe-Akaike Information Criterion (WAIC). The final models for each taxon were: and log(μ anchovy i,j ) = Intercept + ƒ( where ƒ() represents the smoothing functions. After applying the models, Pearson residuals were tested for homogeneity using residual plots (Fig. S4a,b), normality was assessed using histograms (Fig. S4c,d) and spatial independence was determined using variograms (Fig. S5) as well as Pearson residuals plotted on a map (Fig. S6). Residuals were plotted against all co-variables used and not used in the model to detect missing relationships or non-linear relationships to eliminate possible smoother terms. These plots highlighted the necessity of using smoothing terms for the above-mentioned co-variables temperature and latitude. Under-or overdispersion (p ≈ 1) of the final models were assessed using a simulation study simulating 1000 data sets that produced zero-inflated negative binomial distributed data from the final models to compare the dispersion parameter from the original model with those produced by the simulated data sets (Fig. S7a,b). These simulated data sets were also used to examine if they produce a similar number of zeros, as the original data are zero-inflated (Fig. S7c,d). The aim of using the spatial random field, the smoothing functions and the hurdle model approach was to minimise violations of the underlying assumptions of the statistical techniques (Zuur et al. 2010).
To evaluate and interpret model performance, we plotted observed versus fitted values from the final models. Additionally, we simulated n = 128 (based on number of sampling stations) new count data from the final models based on regular spaced temperature values from CTD measurements (between minimum and maximum mean temperature values from 1−25 m depth) and additional mean values of the other covariables that were used in the models. The relationships of temperature and simulated count data from the models were plotted with 95% credible intervals (CIs) to address the high variation in the presence data.

Hydrology
Environmental conditions differed considerably between the 2 survey legs. A rapid change was observed from Leg 1 with upwelling conditions to Leg 2 with an upwelling relaxation event (Fig. 2). Dominant north-westerly winds with mean speeds of 16 m s −1 were present from 1−8 May (Fig. S2). The winds diminished and subsequently changed to westerly winds with mean speeds of 10 m s −1 from 8−15 May. The upwelling plume in Leg 1 north of the Cap-Vert peninsula reached mean minimum SSTs of 20°C (Fig. 2a). The upwelling plume southward from the southern side of Cap-Vert resembled a typical tongue, stretching towards 12.5°N at the border between Senegal and Guinea-Bissau. The plume was characterised by increasing mean SSTs from about 20.5− 22.5°C and occurred roughly between the 20 and 100 m isobaths. Farther in-and offshore, SSTs increased gradually with values reaching 24°C in the study area. This southward extension of the up welling plume occurred until 10 May, when the up welling process ceased. During Leg 2, the upwelling plume subsided, and warmer conditions with a relaxation of the upwelling were observed (Fig. 2b). SSTs < 23°C were only found around and north of Cap-Vert. Gradual increases in SSTs occurred on the zonal axis from north to south with a maximum SST of 26°C.
In Leg 1, the upwelling plume along the whole stretch from 12.5−16°N was accompanied by high values of mean SSCa concentrations ranging be - distribution of egg and larval abundances (no. per 10 m 2 ) of anchovy (Engraulis encrasicolus) and sardinella (Sardinella aurita and S. maderensis) during (e,g,i,k) Leg 1 and (f,h,j,l) Leg 2. Isobaths as in Fig. 1 tween 15 and 25 mg m −3 north of Cap-Vert and 5− 15 mg m −3 south of Cap-Vert (Fig. 2c). This dense stretch of phytoplankton was absent in Leg 2, with higher concentrations only on the coast north of Cap-Vert (Fig. 2d). During Leg 2, remaining patches of phytoplankton occurred from the coast off The Gambia southwards as well as a patch in the southernmost part of the study area off southern Guinea-Bissau, with mean concentrations between 5 and 10 mg m −3 .

Spawning responses: upwelling vs. relaxation conditions
In total, we collected 10 023 anchovy eggs and 1050 larvae as well as 2144 sardinella eggs and 2760 larvae. During Leg 1, anchovy eggs were collected at 4 stations south of Cap-Vert (Fig. 2e) inside the upwelling plume with the lowest SSTs (< 21°C) and abundances ranging from 8−19 785 eggs per 10 m 2 . During Leg 2, anchovy eggs were absent from the collections (Fig. 2f), suggesting a spawning pause in this area or a possible migration under the upwelling relaxation event.
Abundances of anchovy larvae in Leg 1 ranged between 3 and 492 larvae per 10 m 2 , with higher values mainly in the core part of the up welling north of approximately 13°N (Fig. 2g). In Leg 2, the abundance of larval anchovies ranged from 3−337 larvae per 10 m 2 with distributions limited to approximately 14° N; only 2 stations with low abundances extended farther southward along the coast (Fig. 2h). Confirmed stations with anchovy larvae present decreased from 34% (n = 20) to 27% (n = 15) from Leg 1 to Leg 2.
In Leg 1, 5 stations confirmed sardinella eggs (19−100 eggs per 10 m 2 ), at SSTs between 21 and 22°C at the edges of the upwelling zone; only one station confirmed eggs at < 21°C where an chovy eggs also occurred (Fig. 2i). Relaxation of the upwelling in Leg 2 induced a more frequent occurrence of sardinella eggs (8−1812 eggs per 10 m 2 ) collected at 19 stations. The eggs were present at SSTs between 22.5 and 25°C, mostly in the area where the upwelling plume was observed in Leg 1 (Fig. 2j). While sardinella larvae in Leg 1 were found in lower numbers (2−527 larvae per 10 m 2 ), mainly at the margins and apart from the upwelling core (Fig. 2k), the relaxation of the upwelling in duced an intense spawning event represented by in creased larval abundances (2−952 per 10 m 2 ) with a higher presence at SSTs between 22.5 and 26°C (Fig. 2l). The proportion of stations where sardinella larvae were present increased from 33% (n = 19) to 84% (n = 47) from Leg 1 to Leg 2.
Mean (± SD) standard length measurements of anchovy and sardinella larvae were 8.8 ± 2.8 mm (n = 442) and 7.1 ± 2.1 mm (n = 1183), respectively (Fig. 3). Anchovy larvae sampled across T1 were mostly distributed in the upper 25 m of the water column, with lower concentrations in deeper strata (25−100 m) and larvae being absent at >100 m. Abundances were limited in deeper layers at the offshore stations. Sardinella larvae were concentrated in the 0−25 m layer in both T2 and T3 (Fig. S1).

Model results
The hurdle models for both species included either an important temperature smoother in the Bernoulli 92 Fig. 3. Histograms of the total length (mm) of anchovy Engraulis encrasicolus and sardinella (S. aurita and S. maderensis) larvae; number of larvae measured (n) for each taxon and the mean of the length distribution (red dashed line) or the negative binomial part (Table 1). This confirms that early anchovy and sardinella larvae oc cupy specific temperature regimes that represent upwelling and relaxation conditions (concurrent ad vection of tropical water mass), respectively. Our sampling covered areas beyond the minimum and maximum range of larval sardinella occurrence and resulted in a temperature smoother that emulated a typical dome-shaped relationship (Fig. 4a). This relationship resembled a temperature window ranging from 20−23.2°C, peaking at about 22.2°C for sardi nella. In contrast, anchovy larvae were mainly present at temperatures below 19°C (Fig. 4b). Another small peak was observed for anchovy at a range between 21.8 and 23.2°C, representing a few stations surrounding the Cap-Vert peninsula.
The fitted values indicated a decent fit of both models with some variation around the mean (Fig.  4c,d). The fit of the models was mainly generated by the spatial random fields that correspond to one or more unknown co-variables without generating a perfect fit that would assume a model overfit. Models without a spatial random field resulted in a violation of homogeneity, normality and spatial dependency in the Pearson residuals (data not shown). The spatial random field differed between species and legs (Fig. 5). In the sardinella model, the maximum distance of spatial correlation was 32.0 ± 9.5 km, while in the anchovy model the maximum distance was 103.0 ± 34.5 km. That means stations inside these dis-tances were correlated, which could not be addressed by the co-variables in use. In the sardinella model, the smaller scaled spatial correlation resulted in negatively correlated stations in Leg 1 north of 13°N with some positive stations farther south (Fig. 5a). This pattern was reversed in Leg 2, where mainly positively correlated stations were observed north of 13°N (Fig. 5b). The larger maximum distance of spatial correlation in the anchovy model resulted in larger scaled correlations, with negative correlations at coastal stations north of Cap-Vert (Fig. 5c). Farther offshore stations were positively correlated. A band of negative correlations was found between 14 and 13.5°N that extended southward through stations close to the coast. The rest of the stations were positively correlated. In Leg 2, there was a split at 14.5°N, with positively correlated stations northwards and negatively correlated stations southwards (Fig. 5d).
The 2 final models for each taxa were slightly underdispersed (0.79 for the sardinella model and 0.74 for the anchovy model; Fig. S7a,b). The result of this underdispersion was an underestimation of the total count shown by plotting observed versus fitted larval counts (Fig. 4c,d). However, the underdispersion was inside the range of possible dispersal values given by the simulation studies. The 1000 simulated data sets produced enough zero counts for the sardinella model and slightly over estimated zeros in the anchovy model (Fig. S7c,d). While the mean larval counts were very low for anchovy and higher for the sardinella model, both models resulted in large 95% CIs due to the large variation of the count data and the use of the negative binomial approach in the presence part of the model. These mean values and their CIs indicated the asynchronous occurrence of sardinella and anchovy larvae at different temperature windows. Only a very small overlap at temperatures between 19 and 20°C assumed a shared temperature regime (Fig. 6). The remaining unshared temperature windows resembled the temperature smoothers (Fig. 4a,b) and highlighted the asynchronous spawning behaviour of both species inside differing temperature ranges.

Adult stocks distribution and sexual maturity
Sardinella was the most widely distributed fish group observed acoustically over the shelf (depths > 20 m) during both legs (Fig. 7a,f); however, the distribution between the 2 legs differed. In Leg 1, sardinella had a wider but lower density distribution south of Cap-Vert (mean NASC = 227, excluding zero values) with the highest concentrations at the inshore border of the survey area (~20 m depth). During Leg 2, the distribution of sardinella south of Cap-Vert was closer to the coast. The distribution area was smaller compared to Leg 1 and the density in creased (mean NASC = 1257, excluding zero values), with the highest concentrations found at 30−50 m bottom depth. Targeted trawling based on these acoustic observations amounted to a maximum catch rate of 6062 kg h −1 at Stn 7 for S. aurita. This species was found at all stations south of Cap-Vert in Leg 1, while S. maderensis was caught only at the most northern (Stn 6) and southern trawl stations (Stn 2) (Fig. 7b,c). In Leg 2, S. maderensis was mainly collected south of 14°N, while S. aurita was found in higher numbers south of and close to Cap-Vert and off Guinea Bissau (Fig.  7g,h). Anchovies were ob served acoustically only in Leg 1 (Fig. 7d), mostly in a small area south of Cap-Vert. Trawling in that area (Stn 7) also yielded high numbers (11 747 kg h −1 ). The acoustic backscatter signal was weak in other regions and trawling only 94 Fig. 4. Results of the temperature smoothers from the (a) sardinella and (b) anchovy models as well as observed versus fitted value plots of (c) sardinella and (d) anchovy with 95% credible intervals resulted in a low catch rate of anchovies at Stn 2. Trawl catches re vealed the presence of mature adults for all caught species (Fig. 8). GSI values were found relatively higher for the mature fished specimens of S. made ren sis in Leg 2 compared to Leg 1 (Fig. 8).

Asynchronous spawning responses: anchovy vs. sardinella
Short-term environmental fluctuations are difficult to predict and trace in the field. Thus, fish spawning responses under such circumstances cannot easily be evidenced. The present study in the Canary Current upwelling ecosystem is likely the first of this kind. The spatiotemporal design of our survey allowed us to explore the spawning responses of European anchovy Engraulis encrasicolus and sardinella (Sardinella aurita and S. made rensis) during a rapid shift in environmental conditions. The choice of the survey period in May increased the chances of sampling during the changing regime from upwelling to a warmer relaxation regime at the southernmost part of the Canary Current up welling ecosystem. Keeping track of the en vironmental conditions by means of remote sensing allowed a distinct sampling allocation during 2 different consecutive environmental regimes. These regimes involved a change from cold upwelling conditions to a warmer upwelling relaxation event likely coupled with northward advection of tropical waters (Marcello et al. 2011). This shift seemed to trigger a spawning pause in anchovy and a possible migration away from the study area. Simultaneously, low spawning activity of sardinella during the upwelling was re placed by intense spawning activity during upwelling relaxation, replacing anchovy as the main reproducing constituent.
In subtropical and tropical oceans, small pelagic fish species are generally confined to regions dominated by the oceanographic regimes most favourable to their feeding and reproduction strategies (Bakun 1996). In the Canary Current up welling ecosystem, an cho vies and sardines prefer regions of intense up welling, while sardinellas prevail in regions dominated by seasonally altering upwelling and stratified periods (Binet et al. 1998. In syn chrony with the annual trade-wind cycle, the boundary between these 2 regimes shifts seasonally from southern Senegal to northern Mauritania (Mittelstaedt 1983). The southward expansion of the up welling re gime occurs from October to May (Arístegui et al. 2009), while the stratified regime expands northwards during the rest of the year. The timing to conduct such a survey was therefore appropriate to investigate rapid and shortterm spawning responses of fishes to changes in the environment. To enable quick reproduction, parts of an adult school must be in spawning condition. That means adult schools require a fraction of mature individuals that progressively carry ripe gonads. This would allow protracted spawning with multiple spawning batches per year (Alheit 1989). Anchovies, sardinellas and sardines have fast maturation that allows them to produce multiple generations within a year (Alheit 1989, McBride et al. 2015. In some periods, there is a peak in the proportion of mature spawners (Shin et al. 1998). In Senegalese waters, spawning peaks of S. aurita are observed from February to May and from October to December (Baldé et al. 2019), while for S. maderensis the reproductive period seems to peak between April and October and from January to February (Ba et al. 2016). Mature adults of both sardinella species were caught in the fish trawls of our survey. During Leg 1, S. aurita presented higher GSI, while the more tropical S. maderensis showed elevated GSI during Leg 2 (Fig. 8), suggesting an aligned maturation to the warmer conditions. As soon as spawning regimes are present, some adults seem to spawn at their present location, which leads to shared habitats for larvae and adults. This coincidence was also indicated by recordings of both life stages at similar locations during the survey. For E. encrasicolus, no local information is available, and due to its Atlantic-wide occurrence, there is variation in peak spawning seasons (Table S2). For example, in the Mediterranean, the reproductive period of an chovy lasts from spring to autumn, with peaks in egg production in June and July (Palomera 1992, Somarakis et al. 2004, while for the population in the Bay of Biscay, spawning occurs earlier, from March to August, peaking in May and June (Motos et al. 1996). In our study, most adult anchovies were caught in a single trawl during Leg 1 in the cold and productive waters south of Cap-Vert. The adult spawning activity in this regime was verified by both elevated GSI values as well as the presence of eggs and larvae in the area. This suggests that, for anchovies, different life history stages share the same habitats by occupying upwelling conditions along the Senegalese coast. Regardless of the peak spawning seasons, the ability to spawn at any time of year increases the chances for reproduction at various locations. Spawning events are estimated to occur more than 50 times per year for anchovy in the Black Sea (Lisovenko & Andrianov 1996). Thus, any time a multiple spawner encounters adequate spawning conditions, a spawning event may occur. Meeting spawning criteria is often linked to changes in ambient temperature (Quaatey & Maravelias 1999, Ettahiri et al. 2003, Takasuka et al. 2008, Abdelouahab et al. 2021. For instance, in the Mediterranean, anchovies prefer temperate conditions between 15 and 22°C for spawning (García & Palomera 1996), while eggs of S. aurita can be found at 20−27°C (Sabatés et al. 2009). Spawning responses for both species in our study area were found to be strongly related to temperature. During the upwelling conditions (Leg 1), anchovy spawning Legend in (a) refers to the hydroacoustic maps and legend in (b) to the catch rates meet the typical spawning conditions for this species, the rapid increase of temperature in Leg 2 seemed to pause its spawning activity. In contrast to anchovy, adult sardinella had an extended distribution but showed only limited spawning in Leg 1. The limited spawning was re stricted to the boundaries of the upwelling core, indicating avoidance of the lower temperature regime found inside. Only with the onset of upwelling relaxation and the increase of temperatures was elevated spawning initiated for sardinella species. These observations suggest that short-term changes in the environment are immediately utilized for spawning. The temperature ranges confirm thermal preferences for sardinella adults in Senegalese waters (i.e. 21−25 °C; Diankha et al. 2015), which is now also indicated for their offspring. The temperature window of sardinella larvae comprises a range of 21−24°C. Recruitment peaks of S. aurita, the more subtropical species, are mainly found in such environments, while those of S. maderensis typically occur at water tempera-tures > 24°C (Diankha et al. 2018). This would suggest that most larvae collected during the 2 legs were S. aurita larvae. However, young larvae (preflexion and flexion stages) represented the bulk of total larvae collected, and therefore we combined all sardinella larvae together. Thus, spawning and larval residence may indicate speciesspecific temperature windows shared by the adults and their offspring.
The question arises of what will happen when the upwelled water masses in that region exceed temperature levels beyond spawning conditions of anchovy, for instance under the significant warming trend being observed in the southern part of the Canary Current upwelling ecosystem (Gómez-Letona et al. 2017). In May, the upwelled water masses appear to constitute the upper limit of anchovy spawning conditions, and with rising ocean temperatures, this may lead to a shortened period of spawning opportunities. On the other hand, if upwelled water masses are already within limits of the spawning conditions for sardinella, this will inevitably increase spawning opportunities and would essentially lead to an enhanced disparity between the anchovy and sardinella stocks regionally. Such fluctuations are known to occur in many other upwelling areas, where natural climate cycles control stock sizes and hence the dominance structure in small pelagic fish communities (Lluch-Belda et al. 1991, Schwartzlose et al. 1999, Shannon et al. 2004, Lindegren et al. 2013, Mhlongo et al. 2015, Checkley et al. 2017. Species-specific preferences may be better met at particular environmental conditions and may inevitably favour one species over another, and therefore, dominance structures between coexisting small pelagic fishes in large ecosystems can often be asynchronized. Such asynchrony, for instance, occurs in the 3 major upwelling ecosystems, California, Humboldt and Benguela, where the small pelagic fish community is dominated by Engraulis spp. and Sardinops sagax (Checkley et al. 2009).
Asynchronous fluctuations in abundance between small pelagic fishes, however, do not always occur (Siple et al. 2020 (Parrish et al. 1989), S. pilchardus becomes a seasonal resident south of Mauritania during the upwelling and is replaced by the 2 main sardinella species (S. aurita and S. maderensis) in summer (Braham & Corten 2015). This specificity separates the Canary Current upwelling ecosystem from the other 3 upwelling ecosystems and may therefore deviate from the latter in terms of the asynchrony of population sizes among the various small pelagic fish species due to the different small pelagic fish community.
Reproduction of small pelagic fish is strongly associated with the availability of food resources (mainly zooplankton) (Quaatey & Maravelias 1999, Somarakis et al. 2004, Abdelouahab et al. 2021. Income breeders, such as anchovy and sardinella, allocate energy for spawning primarily from feeding (Somarakis et al. 2004), while foraging by the planktivorous larvae ensures their survival. Feeding studies on adults and larvae of S. aurita in the Mediterranean Sea have revealed that zooplankton prey items (mostly crustaceans) comprise the main diet of the species (Tsikliras et al. 2005, Morote et al. 2008, while phytoplankton is rarely found in the stomachs of adults (Tsikliras et al. 2005) and very early larval stages (Morote et al. 2008). Mesozooplankton has also been found to be a main source of dietary carbon for the European anchovy (e.g. Costalago et al. 2012, Nikolioudakis et al. 2014). Using average fluorescence values within the 0−25 m surface layer as a proxy for trophic conditions, we detected only a slight negative correlation for sardinella. During Leg 2, when the fluorescence values were much lower over the surveyed area, the larvae of sardinella were found at higher numbers and showed a wider distribution. Although no adult anchovies or their eggs were recorded in Leg 2, early larval fish were still sampled around Cap-Vert (average length: 8.4 ± 2.5 mm). These findings reveal that the proxy used did not reflect the actual food resource fuelling the reproduction of both species. It is therefore reasonable to assume that the duration of the survey has allowed for active zooplankton feeding and reproduction, hence channelling energy from primary producers to both adult spawners and their offspring.
In the southern part of the Senegalese−Mauritanian upwelling system, spawning of S. aurita occurs during the upwelling season, suggesting that reproductive strategy in the area has evolved as a trade-off between retention processes and food availability (Mbaye et al. 2015). Since S. aurita larvae are mainly found in the upper layers of the water column (uppermost 20 m), they are exposed to the risk of advective loss at the offshore side of the upwelling plume . However, since the presence of larvae coincides mainly with upwelling relaxation phases, the risk of an offshore loss may be reduced. In comparison, anchovy larvae occur mainly in the uppermost 10 m during the upwelling , which suggests a different dispersal strategy. With upwelling conditions, anchovy larvae drift with the productive water masses, whereas sardinella larvae are exposed to increased retention conditions.

Modelling evaluation and interpretation
In our models, we used larval fish data as a proxy for the spawning activities of anchovy and sardinella. Anchovy eggs are known to develop to hatching within 36 h at 22°C (Bernal et al. 2012), whereas for S. aurita eggs there is an estimated time before hatching of 24 h at 24−29°C (Ditty et al. 1994). Even under the lower temperature regime of the upwelling (Leg 1), the eggs would be expected to have a fast development and short period of occurrence. This is not the case for fish larvae, as they have a longer persistence in the water column (Ditty et al. 1994). Since each leg of the cruise lasted 10 d, spawning would likely be underestimated if just considering eggs. Therefore, we compromised by using the larval fish stages, even though dispersal processes may have occurred and resulted in a wider distribution pattern than actual egg spawning.
The differing temperature preferences for sardinella versus anchovy were identified by means of observation and statistically confirmed by the hurdle models which included a spatial random field. The spatial random field accounts for variation in the count data and adjusts the model to be inside a valid statistical framework for regression-type analyses (Zuur et al. 2010). Spatial autocorrelation accounts for any potential co-variables that were not measured during the survey. It represents a 'black box' designated as 'w', which is the spatial correlation of Multinet MIDI stations that corresponds to spatially correlated larval abundances (Fig. 5) to improve the credibility of the models (Hoeting 2009). The complexity of population dynamics involves many traits that are difficult to detect in the field (Houde 1987). Let us assume different schools of fish allocated heterogeneously in an ecosystem. In those schools, different numbers of mature fish may be present, and the feeding rate between schools as well as the swimming speed and direction may also differ. When eggs are released at one location on a specific day by a fraction of the school, and several other fractions of mature specimens release eggs in the following days, this will result in spatially and temporally correlated patches of fish eggs and larvae (McGurk 1986). In other words, stations closer to each other in space and time are likely more similar than stations farther apart. These traits are difficult to measure in fieldbased studies. Therefore, using spatial random fields, which is now possible with software packages like R-INLA, is a promising development for ecological studies as this method can cope with spatially and temporally autocorrelated data sets. Without considering spatial and temporal autocorrelation, it is likely to infer relationships (between response variables and co-factors) that are deemed important. However, by neglecting dependencies, we may underestimate credible intervals (CIs) (Hoeting 2009). Here, CIs are large, which makes sense when considering the high variability in the positive count data and the usage of the dependency structure in our models. R-INLA is a great advancement that allows calculating models with low computational cost (seconds to minutes in many cases) and high accuracy  while accounting for dependency structures in ecological data sets.
In summary, our study combined survey data to unravel how fast-changing environments are utilized for spawning in dominant small pelagic fishes (anchovy and sardinella). Modern statistical techniques suggest species-specific spawning windows that indicate an important asynchronized life history trait. This life trait could potentially occur in all world oceans where cold upwelling conditions alternate with warm relaxation events. These findings provide a deeper insight into the population dynamics of small pelagic species in upwelling ecosystems.
Data availability. The data collected within the framework of the EAF-Nansen Programme can be requested through the FAO webpage https://www.fao.org/in-action/eaf-nansen/ data-access-requests/en/.