© D. Ahmed-Babou et al., Published by EDP Sciences 2023

1 Introduction

Ecological monitoring and information systems are increasing worldwide. They are key to providing long-term observations, which are particularly important for reference point determination and investigating tipping points and regime shifts. Aside from their maintenance, one of the main challenges associated with these monitoring systems is how to extract scientific knowledge. For example, how to detect spatial patterns that are persistent in time from large sets of georeferenced observations and how to evaluate their temporal stability in conjunction with abundance variability.

This study aimed to develop a generic approach for summarizing the overall spatiotemporal information of a long time series of spatial distributions into a portfolio of distribution maps with two elements, one time-dependent and the other, possibly low dimensional, space-dependent. The approach is based on min-max autocorrelation factors (MAF) (Switzer and Green, 1984) that can be seen as an extension of empirical orthogonal functions (EOFs, Lorenz, 1956; Wikle et al., 2019) to deal with the presence of spatial structure in the data. MAF have been widely applied in ecology (e.g. Petitgas et al., 2020; Solow, 1994; Woillez et al., 2009). However, MAF were initially designed to filter out small-scale noise from a set of images. The percentage of variance explained by the selected factors is thus unknown. A recent paper by Bez et al. (2022) suggested that MAF can be appropriately reformulated into Empirical Orthogonal Maps (EOMs) that order the factors according to the percentage of variance explained. The aim of the study by Bez et al. (2022) was to identify the main spatial patterns that shape the dynamics of a stock. The present work defines how the EOMs are formulated and selected to construct a portfolio of maps consisting of a set of maps describing the main spatial patterns governing the biomass distribution of a species. It then investigates the mutual temporal fluctuations of the elements of the portfolio.

The ideas developed in this paper are illustrated for common octopus (Octopus vulgaris) off Mauritania over the past thirty years. Octopus represents a key species of the Mauritanian marine ecosystem both in ecological terms (Boyle & Boletzky, 1996; Caddy, 1983) and in economic terms. It accounts for more than 70% of cephalopod landings and one third of demersal catches (Khallahi et al., 2020). In 2019, octopus generated 71% of total export value ($360 million) of fishery products (Société Mauritanienne de Commercialisation Des Poissons “SMCP”, 2020).

The current management and conservation of octopus is based on biological and ecological knowledge, in particular its short lifespan of 12–14 months (Guerra, 1979; Hernández-López et al., 2001; Mangold and von Boletzky, 1973), its rapid growth (Domain et al., 2000; Mangold, 1983; Sanchez et al., 1998), its high fecundity of 100,000 and 500,000 eggs laid by female (O'dor et al., 1978) and its sedentary lifestyle (Hatanaka, 1979; Caverivière et al., 1999). Octopus reproduce only once; males die after mating and females after hatching of offsprings. Therefore, it is essential to ensure a sufficient proportion of the stock has the chance to reproduce to ensure stock renewal. The sedentary nature of the species implies the need for spatialized management to limit the effects local depletion by fishing and thus creates the need for a good knowledge of its spatiotemporal distribution. In Mauritanian waters, two biological rest periods of 1–2 months each have been implemented each year for all demersal fisheries since 2008. These two rest periods reflect the existence of two annual octopus cohorts with distinct spawning seasons (Hatanaka, 1979; Bez et al., 2022). The variability in abundance of this species is mainly a consequence of recruitment fluctuations that are partly dependent on environmental conditions, especially the intensity of upwelling (Caverivière et al., 1999; Otero et al., 2008) and on the two spawning periods taking place under distinct environmental conditions (Faure et al., 2008).

Since 1982, a dedicated monitoring system has been developed to support scientific advice to octopus management bodies (Gascuel et al., 2007). It is based on regular scientific surveys with standardized sampling protocols providing information on octopus spatial distributions through time. Using these data, persistent spatial patterns or tipping points of change in spatial distribution can then be tracked. In this study, the portfolio of octopus distribution maps in Mauritanian waters was characterized using data from 63 scientific surveys covering the period 1987– 2019. The size of the portfolio is discussed as a potential proxy for ecosystem resilience.

2 Materials and methods

2.1 Scientific monitoring surveys

Between 1987 and 2019, the Mauritanian Institute of Oceanographic Research and Fisheries (IMROP) performed 97 scientific demersal bottom trawl surveys on the continental shelf. However, only 63 surveys covered the entire continental slope and were therefore kept for analysis in this study (Fig. 1a). The sampling protocol was standardized over the entire period and followed a stratified random protocol based on three latitudinal strata (Fig. 1). The mean number of hauls was 102 hauls per survey (min = 57, max = 205, but only four exceed 121 hauls), equally distributed in each strata. The survey area was divided into five bathymetric strata (<30 m, 30–80 m, 80–200 m, 200–400 m and 400–600 m). The sampling covered very coastal areas due to the probable presence of coastal demersal species. The depth covered ranged from 10m to around 600m. The area of the Banc d'Arguin National Park (PNBA) was not covered by the sampling because it is a marine protected area closed to motorized fishing. The observations consisted of the octopus densities expressed in terms of number of individuals per square meter obtained by dividing the catches by the surface of the area swept by the trawl. These were considered as indicator of abundance assuming that gear catchability remained constant. More details on these data are available in Gascuel et al. (2007), Meissa et al. (2013) and Meissa and Gascuel (2015).

For each survey, haul positions were randomly selected and mutually independent. The hauls were thus not performed at the same locations each year. However, the method used in this study required that the observations are located at the same positions (regular or not in space) during the entire time series. In order to fulfill this requirement, the data were kriged on a 0.1° × 0.1° grid covering the sampling area before analysis (341 grid cells). Kriging allows estimating the mean octopus density within each grid cell (Chilés and Delfiner, 2009, provided that there were enough hauls in the vicinity. Variogram models were fitted, survey by survey, to the empirical variograms (for the sake of parsimony they are not shown). The input data for the analyses were thus 63 kriging maps and the corresponding abundance indices.

A time series of upwelling strength over the same period was built from a monthly dataset of the meridian component (North-South) of composite wind speeds over the Mauritanian zone, 50–100 km from the coast. This proxy is considered the best proxy linearly linked to the upwelling index for a North-South oriented coast (Demarcq and Faure, 2000), which is the case for the study area.

A spatialized time series of seasonal average sea surface temperature (two maps per year, one for the hot season and one for the cold season) on the same spatiotemporal grid as the survey data was built from the National Oceanic and Atmospheric Administration (NOAA) database (NOAA, 2020).

Four climatic seasons are generally distinguished in the Mauritanian zone (Dobrovine et al., 1991). A cold season from January to May, a cold-warm transition season from June to July, a hot season from August to October, and a hot-cold transition season from November to December. For this study, these four seasons were aggregated into two main seasons of five and seven months respectively: a hot season from June to October, including 28 of the 63 available surveys and a cold season for the rest of the year (November–May), including the remaining 35 surveys (Fig. 1). The precise timing of the surveys within each season was opportunist rather than following a predefined sampling protocol. The temporal spreading of the surveys was however of the same level in each season (similar coefficients of variation; Fig. 1). A portfolio of maps was considered for each season. All the calculations were performed in R using the package RGeostats (MINES ParisTech / ARMINES (2022)).

Fig. 1 (a) Spatial distribution of sampling stations during the study period. PNBA = Parc National du Banc d'Arguin. The three latitudinal strata (North, 20°36N−19°15N; Central, 19°15N−17°45N; and South, 17°45N−16°03N) are represented by dashed lines. (b) Timing of the surveys (month) for each of the hot and cold seasons. The mean and coefficient of variation of the timing of the surveys by season are given in blue. (c) Number of surveys per year.

2.2 Setting up the portfolio

A portfolio is defined as a selection of relevant empirical orthogonal maps (EOMs). Below the definition of EOMs is recalled and the method to select those included in the portfolio is explained.

2.2.1 Building EOMs

The workflow for creating empirical orthogonal maps (Bez et al., 2022) is depicted in Figure 2. EOMs are a variation of min-max autocorrelation factors (Switzer and Green, 1984), where the factors are ordered according to their percentage of contribution to the total input variance. The first EOM explains most of the input variance, and so on down to the last EOM that explains little of the input variance.

We denote x_i, i = 1, … , S the geographical positions of the S grid nodes that were systematically informed (isotopy) at time t_j, j = 1, … , T, withT being the number of input maps (S = 341, T = 63). Input data can be formatted as an S × T matrix denoted by Z = Z [i, j] = z (x_i, t_j) where each line in Z corresponds to the time series of the value of a given grid cell. Computing EOMs consists of two sequential PCAs. The first PCA, which is nothing but an EOF based on Z, produces an S × T matrix, denoted Y, with T uncorrelated columns y (⋅ , t_j):

(1)

In the spatio-temporal context, each column of Y is a set of georeferenced points, i.e. a map, uncorrelated to each other. The absence of correlation refers to pointwise correlations, therefore, to the absence of correlation between map values for the same geographical point.

The second PCA aims at building T new factors F_j (x_i) = F (x_i, t_j) , j = 1, … , T, without mutual spatial correlation for a given spatial distance r (the unit spatial lag corresponds to the grid cell size of the input maps), that is, with the following null scalar products:

(2)

This second PCA is based on the eigen-decomposition of the matrix of variogram and cross-variogram values between standardized EOFs from the first PCA. Finally, we obtain a set of T new maps called empirical orthogonal maps (EOMs) ordered by decreasing percentage of explained variance. This realizes the decomposition of the initial spatio-temporal variables z (x_i, t_j) into a product of two elements, one element purely spatial made up of the EOMs F_k (x_i) and the other element made up of their temporal weights ψ_k (t_j).

(3)

where m (t_j) is the average abundance of the input map for time t_j. The weights measure the importance of a given EOM for explaining the survey time series. In other words, the initial set of maps is decomposed into a weighted linear combination of maps without correlation at short scale (blue box in Fig. 2). At this stage, the decomposition is performed with no loss of information. EOMs can be computed for raw or standardized maps. The pseudo-algorithm provided in supplementary material highlights the main steps for computing EOMs when the first PCA is applied to standardized maps.

The signs of the values of the EOMs and their weights are conventional. Therefore, the interpretation of EOMs and their weights must be established jointly. For example, the interpretation of an EOM whose weights are all negative in the decomposition is strictly equivalent to its symmetry with respect to 0 when taking the symmetry of both the EOM and its weights.

By construction, EOMs are mutually uncorrelated for distance 0 and distance r. Assuming that they are also orthogonal for all distances implies that they are fully spatially orthogonal, and that they can be considered as basis for map decomposition. This assumption could be of concern if we were interpolating (i.e. kriging) several EOMs obtained from irregular points in space. In such a case, rigorous interpolation should be done by cokriging the different EOMs. In this study, interpolation was done before the decomposition and we were simply taking linear combinations of the interpolated input maps. Thus, two contrasting situations can be considered: “interpolate first and decompose then” (as in the present case study) where the assumption has no consequence, and “decompose first and interpolate then”, in which case, co-kriging should be recommended to provide a spatially interpolated version of the EOM.

Fig. 2 Workflow of the method. Upper blue block: building T different EOMs from a set of T maps. Lower orange block: selecting the P EOM entering the portfolio. In this case, the Q-first EOMs explain a given percent of the overall variance. Amongst them, the second EOM is not retained as its weights are not sufficiently correlated to the temporal fluctuations of abundance.

2.2.2 Building the portfolio

The selection of EOMs for the portfolio is based on two nested steps (dashed orange box in Fig. 2). In step one, as it is standard practice in PCA analyses, the number of selected EOM is set to ensure a sufficient part of the variance of the input maps is retained. In the present work, we selected the Q EOMs explaining at least 60% of the input variance. The 60% threshold is arbitrary. It is a compromise between complexity (i.e. the number of EOMs) and explained variance. In our case study, beyond 60% the number of factors increased considerably without a significant gain in explained variance.

In step two, among the Q EOMs from step 1 only those are retained whose spatial patterns are connected to the population dynamics of the stock, that is their temporal weights ψ_k (t_j) are correlated with abundance. This is achieved using three complementary criteria. The first criterion is a statistically significant linear correlation between the EOMs' weights and abundance (test level 5%). The second criterion is based on a multivariate time series analysis using auto- and cross-correlograms to identify EOMs whose weights vary smoothly over time (auto-correlogram of weights) and in a consistent manner with regards to the temporal variation of abundance (cross-correlogram between weights and abundances). For this, a linear model of co-regionalization (Chilés and Delfiner, 2009; Goulard and Voltz, 1992) consisting of a nugget effect and a linear part is fitted to the empirical correlograms. This allows modeling and quantifying in a consistent manner the random part and the temporally structured part of the correlations. EOMs to be retained have a small random part (small percentage of nugget effect in the respective auto-correlograms). The third criterion is based on maximizing the covariation in time between abundance and EOM weights. This criterion is evaluated using the slopes of the cross-correlograms between weights and abundance. In linear models of co-regionalization, the (absolute) value of the slope of the cross-correlogram must be smaller than the square root of the product of the slopes of the two auto-correlograms (Chilés and Delfiner, 2009). Thus, the third criterion is that the ratio between the slope of the cross-correlograms and this upper limit is near one. The final selected set of P EOMs (P ≤ Q) constitute the portfolio of maps summarizing the spatial dynamics of the species for a given season (Fig. 2). Note that when a portfolio is made up of several maps, their temporal joint dynamic is accessible through the auto- and the cross-correlograms that were used for the selection of the EOMs entering the portfolio.

For each of the two seasons, this methodology is used to build portfolio for octopus and for SST separately.

2.2.3 Climatology

A climatological spatial distribution pattern can be obtained computing weighted averages of the P maps in the portfolio:

(4)

where is the mean temporal weight of the P selected EOMs such that:

(5)

and where is the mean abundance:

(6)

3 Results

3.1 Portfolio of distribution maps for the hot season

The hot season is the main octopus breeding period (June to October). The first EOM alone restored 47% of the variance of the observed abundance data (Fig. 3), and the first three EOMs recover 60% of the initial variance. Amongst them, only the first one was finally selected, the two others ones being not sufficiently correlated with abundance, both statistically and temporally (small coordinate values in Fig. 3), and having poorly organized temporal evolutions (small feature size in Fig. 3).

The hot season's portfolio thus consists only of the first EOM, which is characterized by a north-south gradient, with a high-density area in the north (especially in its wide part), a moderately dense area in the center, and a low density area in the south (Fig. 4, left). The temporal evolution of the weights of the first EOM shows three different phases. There is a sharp decrease over the first five years followed by stability at low level for about 15 years and then a slightly higher level. Most of the signal describing the spatial distribution of octopus abundance is represented by this pattern. This distribution pattern provides information on space occupancy of octopus during this season. It summarizes the historically known pattern of octopus distribution with decreasing abundance from north to south. The years where this pattern had a large weight were years of high octopus abundance during the hot season such as the beginning of the study period. The spatial distribution represented by this EOM is linked to the main breeding season of the species, which takes place at the end of the hot season (September-October). During reproduction, octopus seems to display an affinity to milder water temperatures, which occur to the North of the study area (Fig. 5, left).

Fig. 3 Criterion used to define the distribution map portfolio for octopus in Mauritanian waters for the hot (left) and the cold (right) season. The inset panels represent the percentage of variance explained when including progressively more and more EOMs. The threshold of 60% was used for selecting the candidate EOMs for building the portfolio. The main graph shows the position of EOMS with respect to their relationship with octopus abundance. The x-axis corresponds to the correlation between EOM weights and abundance (the closer to +/− 1, the better). The y-axis represents the relative values of the slopes of the cross-correlograms between each EOM weights and abundance, with respect to the maximum possible value (the closer to +/− 1, the better). The size of EOF numbers is inversely proportional to the percentages of nugget effect in the auto-correlograms of the EOMs' weights (the larger the better). The EOMs that allow recovering 60% of the input variance are indicated in red and connected for better visibility.

Fig. 4 Portfolio of octopus during the hot (left − red) and the cold (right − blue) seasons over the entire study period. Each plot is made of an EOM with the temporal evolution of its weights, smoothed by a loess method with a span parameter value of 0.5 and associated 95%-confidence interval. Pie chart indicate the percentage of variance explained by the EOM (arrows represent the cumulative percentage of variance explained by the current EOM and the previous ones and the percentage of variance explained by the portfolio). The EOMs that were not selected for the portfolio are represented in grey.

Fig. 5 Portfolio of SST during the hot (left − red) and the cold (right − blue) seasons over the entire study period. Each plot is made of an EOM with its weights (smoothed by nonparametric loess with a span parameter value of 0.5 and associated 95%-confidence interval) and a pie indicating the percentage of variance explained by the EOM (arrows represent the cumulative percentage of variance explained by the current EOM and the previous ones).

3.2 Portfolio of distribution maps for the cold season

During the cold season the spatio-temporal complexity and thus the size of the portfolio was larger than during the hot season (Fig. 3, right). No single EOF spatial pattern was found to be dominant. Indeed, the first EOM explained only 14.7% of the variance of the input distribution data, while the first nine EOMs restored 60% of the variance. Four of these top nine EOMs, respectively the 1st, 3th, 4th and 5th, were linked to the spatiotemporal dynamic of the observed mean gridded survey abundance (density). Their weights depicted significant temporal cross correlation with mean abundance. Together, the four finally selected EOMs to build up the portfolio represented 35.4% of the initial spatiotemporal variance.

The spatial pattern of the first EOM for the cold season was similar to that of the hot season, characterized by an overall north/south gradient with a slightly different temporal evolution (Fig. 4, left). However, this EOM showed a denser abundance off the northern area with a more pronounced southward expansion of moderately dense areas. The second EOM of the portfolio (EOM3) indicated a higher abundance in the center with a concentration of abundance south of Nouakchott and two lower density areas in the far south and off the Banc d'Arguin. The EOM4 displayed across-shelf gradient patterns, opposing low density coastal areas to denser offshore areas, particularly in the north and center but also, to a lesser extent, in the south. The last EOM in the portfolio (EOM5) showed two hot spot areas in the north and center. This spatial pattern reinforced, in the north, the spatial pattern of the first EOM but complemented it in the central and southern zones.

Given the criteria used for their selection, i.e. correlation with abundance, the EOMs selected for the portfolio had similar long-term temporal variations in their weights (Fig. 4, right). The weights decreased at the beginning of the period, then flattened and slightly increased during 2010–2015. This pattern was reflected in the slope of the respective correlograms (Fig. 6). However, as indicated by the (quasi systematic) absence of nugget effects in the cross-correlograms, their fluctuations at short temporal scales were not correlated (Fig. 6). A noticeable exception concerned the short term negative cross-correlation between EOM1 and EOM3. This reflected the opposite fluctuations of the empirical time series in the middle of the period.

The model used in this study quantified the joint temporal dynamics of the EOMs of the portfolio and mean seasonal abundance. It decomposed their mutual temporal correlation into two parts: a nugget effect and a slope. The nugget effect quantified short-term joint fluctuations (e.g. the average variations between two consecutive cold seasons). The slope characterized long-term (multi-annual) joint evolutions. In this context, a change in the overall correlation can be due to either short or long-term processes, or both. However, by construction, the EOMs of the portfolio have been selected because they shared long-term dynamics with that of abundance. By consequence, it remains that within the portfolio the various levels of correlation between EOM weights and abundance rely strongly on the short term. This was consistent with the strong correlation found between EOM' weights and abundance and the value of the nugget effect (Fig. 6). In particular during the cold season, the weights of EOM3 were the most strongly correlated to abundance of octopus specifically because, in addition to the fact that they shared similar long-term variations, their short-term fluctuations were also correlated. The cross-correlogram between EOM3 and abundance was indeed proportional to the auto correlogram of abundance (Fig. 6). This means that the residual from the regression of the weights of EOM3 against abundance was pure noise, i.e., residuals with no temporal auto-correlation (self-krigeability; Chilés and Delfiner, 2009). In practice, this means that the weights of EOM3 were proportional to abundance over time with some uncorrelated noise (Fig. 7). The third EOM represents thus a central spatial pattern of the spatio-temporal dynamic of octopus. It explained around 7% of the spatial distribution of octopus that is density dependent.

The climatological spatial pattern associated with the cold season was characterized by two dense areas rather offshore in the northern zone and coastal in the central zone (Fig. 8). A northern zone off Cape Blanc had its epicenter located between the bathymetric lines of 80 and 200m. A second zone was found south of the Cap Timiris (in the middle of the central zone) with an extension to a bathymetric level close to that of the northern zone.

The weights of the EOMs that were part of the portfolio were not statistically different (Student p-value = 0.008) for surveys carried out in December and in March (Fig. 9). The ranges of fluctuation and the mean densities appeared smaller for December surveys. However, for the two instances where a December survey was followed by a March surveys the following year, the evolutions of the weights were mixed up (increasing for EOM5, decreasing for EOM3, one increasing and the other one decreasing for EOM1 and EOM4).

During the cold season, the two months with the largest number of surveys were December and March (Fig. 1b). Although not significant, the temporal weights for EOM 1 to 5 were larger in March compared to December (Fig. 9). March was the month used for analyzing relationships with the upwelling intensity index. Correlations with the upwelling index for the same month (i.e. considering EOM weights during the months of the surveys and the upwelling indices for the same month during the period 1987–2019) and with a time lag of one to seven months before showed similar behavior for each of the EOMs of the portfolio. The correlation is significantly negative only for EOM3 for a 6 months delay (5% level test; Fig. 10 left). The upwelling indices in September belonged to the lower part of the interval of fluctuation and corresponded to low upwelling intensity (Balguerías et al., 2002). The eight instances available for this study (octopus survey in March and upwelling index available six months before) covered the full range of possible values observed in September from 1988 to 2019. In particular, they were not distributed in the tails and hence did not represent a particular situation.

Fig. 6 Analysis of the temporal dynamics of the principal spatial patterns (EOM 1, 3, 4, 5) and of the abundance for the cold season. Diagonal panels: auto-correlograms; Lower triangle panels: cross-correlograms. For all panels, the x-axis represents time lags in years (from 0 to 30 years) and the y-axis represents correlations. The horizontal black dashed line is at 1. For off diagonal panels, the horizontal dotted lines represent the correlations. The panels of the last line concerns the cross-correlograms of EOM weights with abundance and the auto correlogram for abundance (last panel). The horizontal red dashed lines represents the correlation between EOM weights and abundance. The linear model of coregionalisation is represented as a continuous line (nugget effect + slope). The isolated panel represents the correlation between abundance and EOM weights as a function of the nugget effect in the model.

Fig. 7 Cold season. Time series of the abundance (black) and temporal weights of EOM3 (red).

Fig. 8 Climatology for the cold season. This corresponds to the mean of the EOMs selected in the portfolio weighted by their mean weights. The scale is consistent with the input octopus densities (i.e. nb of individual per m2).

Fig. 9 Portfolio for the cold season. Distribution of the EOMs' weights as a function of the EOM numbers for March and December surveys. The weights for the surveys carried out in March (N = 10) are represented in black. The weights for surveys performed in December (N = 10) are given in red. The average values are superimposed. For the two instances where a December survey is followed by a March surveys three months later, the evolutions of the weights are represented by doted segments.

Fig. 10 Portfolio of the cold season. Correlation between the weights of EOM3 and the upwelling index. Left: lagged correlations between EOM3 weights in March and the upwelling index 6 months earlier. Right: boxplot of upwelling indices for all months except September and September (i.e. 6 months before March). In red the values for the nine years when a survey was performed in March with concomitant upwelling observations.

4 Discussion

The approach of portfolios of distribution maps proposed in this paper relies on a set of criteria to select the empirical orthogonal maps to be included in the portfolio. These criteria may be case specific and it is hard to avoid subjective choices. First, similarly to any PCA analysis, one has to choose a threshold for the variance to be explained by the set of selected maps. For a target value of 60% of explained variance, it can happen, as in our case that three and nine maps can be enough for summarizing spatial distributions of octopus in two seasons. The aim of the study being to identify the spatial patterns which were correlated with abundance, we then selected those EOMs whose weights had significant temporal cross-correlations with abundance. This resulted respectively in one and four spatial patterns in each seasonal portfolio explaining 47% and 35 % of the overall variance observed over the past thirty years. In other words, this analysis demonstrated for the first time that half of the information brought by the initial twenty-eight maps recorded between 1987 and 2019 during the hot season could be summarized by a single distribution map whose temporal evolution increased and decreased with octopus abundance. During the cold season, the situation was more complex and the variability could be less easily reduced. Four maps madee up the portfolio for that season and their recombination based on their respective weights allowed recovering only a third (35%) of the initial variance. There was thus a more clear and persistent signal for connecting spatial pattern and stock abundance during the hot season, and thus more possibilities for efficient spatial planning during this season.

Despites its importance for explaining the overall variance of the distribution maps observed during the cold season, EOM2 was not selected for the portfolio. Its weights were too variable over time and not sufficiently correlated with abundance, both in the short and in the long term (Fig. 3). In other words, EOM2 represented a spatial pattern that was important for explaining the variance of the original set of distributions (hence its rank in the decomposition) but without temporal coherence and without connection to abundance. EOM2 was symptomatic of a coastal hotspot of abundance south of Cap Timiris (Fig. 11) that pulses at random from time to time without contributing much to the abundance. The portfolio of the cold season explained not more than a third of the overall variance, which indicates that the spatial distribution of octopus during the cold season were highly variable and poorly linked to abundance variations.

The derived climatology aimed at representing the average spatial distribution of octopus over the study period accounting for the effects of abundance variations on the spatial distribution. The two seasonal climatologies based on one map for the hot season and on the average of four maps weighted by their average weights for the cold season, showed similar spatial patterns (not shown). This confirms the known yearlong persistence of a gradient of octopus density from north to south (Pease, 1973) whose intensity is related to abundance (stronger in years with higher densities). However, during the cold season, abundance tended to be more offshore with a secondary area of distribution in the middle of the Mauritanian EEZ. The fact that the portfolio was larger during the cold season could be linked to an affinity to optimal water temperatures.

The spatio-temporal decomposition of the time series of thirty-four SST maps covering the period 1986–2020 for each season provided interesting insights on the spatial distribution of octopus and on the two portfolios generated in the present study. First, during the hot season, 88% of the variance of the SST input maps was represented by a single map. During the cold season, the signal was more blurred. Four maps were needed to recover 67% of the variance (the first one incorporating 38%). This difference in number of maps was similar to the difference between the two portfolios obtained for octopus even though the levels of variance explained by the portfolio (48% and 35%) were notably smaller than the percentage of variance explained by the first x EOMs of SST (? %). Second, the main spatial patterns in both seasons was a north-south gradient. However, the gradient did not mean the same thing in both seasons. According to the literature, an optimal window for octopus would likely be centered on 21 °C (Villanueva, 1995). The observed mean SST was around 24 °C and 19.5 °C during the hot and the cold seasons respectively so that the optimal window is respectively below and above the mean during the two seasons. Values around 21 °C corresponded to negative areas in SST EOM1 for the hot season (Fig. 5). More precisely, the SST in the northern area of EOM1 during the hot season represented an area that was −2 °C below the seasonal average. Weighted by 1.5 in the decomposition of the true SST distributions, this area was thus precisely where one finds systematically SST values around 21 °C during the hot season. There was thus a very strong match between octopus densities and SST main spatial patterns during the hot season. In the cold season, this gradient was complemented by other patterns accounting for more than 10% each, that were consistent with the patterns observed in the octopus portfolio (more southern, more offshore). Indeed, given the mean SST, the optimal window, if relevant, corresponded to positive areas in the main maps. In particular, the fourth most important pattern presented in the cold season time series corresponded to a clear and strong westwards cross-shelf gradient that can be connected to a similar medium signal found in the octopus portfolio.

The fact that the arrangement of maps in terms of the percentage of variance explained (EOM) strongly matched the arrangement based on the strength of their spatial structures (MAF) is worth mentioning (Fig. 12). The situation could be that the EOMs explaining most of the variance could be maps without spatial structure, i.e. with hot spots appearing here and there from EOM to EOM. This would be the case for systems with strong spatiotemporal fluctuations and poor persistence of their spatial distribution over time. In this study, the maps explaining the most of the variance of octopus spatial patterns are indeed also those with the strongest spatial structures.

In an ideal situation, one would work with a regular spread of scientific surveys over time. This was not the case in the present study as the timing of the surveys was somewhat opportunistic. For instance, during the cold season, surveys took place mainly in December (N = 10) or in March (N = 10) and April (N = 8). However, the weights of the EOMs were not associated with specific months and their differences were not statistically significant. Therefore, the four EOMs selected for the portfolio of the cold season (November–May) were considered together as the set of principal spatial patterns that best summarized octopus spatiotemporal distributions during that season. However, the analysis of the joint temporal dynamics of the weights of the different EOMs highlighted the particular characteristics of EOM3. While all EOMs of the portfolio shared their long-term trends with that of abundance, the short-term (i.e. year-to-year) variations of the weights of EOM3 were also correlated with abundance. EOM3 represented thus a major spatial pattern of the spatiotemporal dynamics of octopus. It captured a part of the spatial distribution of octopus that is density dependent. The particular role played by EOM3 was reinforced by the analysis of the correlation with the upwelling intensity index. Amongst the four EOMs making up the portfolio of the cold season, this was the only one whose temporal weights showed significant correlation with the upwelling index when considered six months before. Even though to be interpreted with care given it was based on only eight observations, this result is interesting. While the phenology between upwelling and primary and secondary productions is well documented, the demonstration of a delayed impact on octopus distribution is new. The March surveys caught individuals that weight on average one kilogram (data not shown). These individuals mostly correspond to pre-spawning adults of 9-10 months. They come from the spawning that took place around May the year before. Given the duration of larvae and para-larvae developments (Otero et al., 2007, Domain et al., 2000), the end of the pelagic phase for these individuals falls around August-September. Thus September during which upwelling is usually not intense (Fig. 10), coincides with the period of settlement for the individuals spawned in May. Our results suggest that the less intense the upwelling is in September, the more spatially concentrated the octopus recruitment is in the following March, with highest occurrence of octopus off the coast in the central Mauritanian ZEE as indicated by the spatial pattern of EOM3 (Fig. 4).

The present study suggested that a portfolio of distribution maps and in particular its size can provide novel insights. During the hot season, only one EOM was needed to summarize the set of observations for octopus compared to four EOMs for the cold season. The two seasons are of similar durations (five and seven months respectively) and surveys data were spread across both. Thus, the difference in portfolio size between the two seasons is not attributable to variations in the timing of the monitoring. Instead, the difference in map portfolio size could be considered a consequence of a real difference in the spatial dynamics of octopus during the two seasons. The portfolio size quantified the effective number of spatial patterns of octopus seasonal distributions. It was reduced to a single spatial pattern during the hot season in the present study. In contrast, the spatial distribution of octopus during the cold season was a mixture of four principal maps, with the importance (weight) of one being strongly related to the fluctuation of abundance.

Ecological considerations arising from portfolio size can be of two kinds. Portfolio size measures the complexity of spatial occupancy patterns and the associated complexity of possible external parameters governing them over time. On the one hand, a large portfolio can characterize systems where one process governing biomass distribution can decrease without impacting the overall system if other processes offset it. On the other hand, systems associated with a small portfolio size have fewer buffer effects and could be more sensitive to external changes (putting all the eggs in one basket is risky). Following these lines, a general perspective of this work could be to quantify the size of many species distribution portfolios and relate them to knowledge about the resilience of the ecosystem the species belong to.

Fig. 11 Left: Empirical Orthogonal Map EOM2 of octopus for the cold season. This EOM was not selected for the portfolio. Right: temporal evolution of its weights smoothed by a loess method with a span parameter value of 0.5 and associated 95%-confidence interval.

Fig. 12 Comparison between the MAF and the EOM arrangements for the hot (left) and the cold (right) season for octopus in Mauritanian waters. The green connections concern the Q first EOMs explaining 60% of total variance.

Funding

No funding.

Conflict of interest statement

No conflict of interest.

Author contribution statement

DAB and NB conceived the ideas, analysed the data and led the writing of the paper. HD and BM contributed critically to the draft. All authors gave final approval for publication.

References

All Figures

Fig. 1 (a) Spatial distribution of sampling stations during the study period. PNBA = Parc National du Banc d'Arguin. The three latitudinal strata (North, 20°36N−19°15N; Central, 19°15N−17°45N; and South, 17°45N−16°03N) are represented by dashed lines. (b) Timing of the surveys (month) for each of the hot and cold seasons. The mean and coefficient of variation of the timing of the surveys by season are given in blue. (c) Number of surveys per year.

In the text

	Fig. 2 Workflow of the method. Upper blue block: building T different EOMs from a set of T maps. Lower orange block: selecting the P EOM entering the portfolio. In this case, the Q-first EOMs explain a given percent of the overall variance. Amongst them, the second EOM is not retained as its weights are not sufficiently correlated to the temporal fluctuations of abundance.
In the text

Fig. 3 Criterion used to define the distribution map portfolio for octopus in Mauritanian waters for the hot (left) and the cold (right) season. The inset panels represent the percentage of variance explained when including progressively more and more EOMs. The threshold of 60% was used for selecting the candidate EOMs for building the portfolio. The main graph shows the position of EOMS with respect to their relationship with octopus abundance. The x-axis corresponds to the correlation between EOM weights and abundance (the closer to +/− 1, the better). The y-axis represents the relative values of the slopes of the cross-correlograms between each EOM weights and abundance, with respect to the maximum possible value (the closer to +/− 1, the better). The size of EOF numbers is inversely proportional to the percentages of nugget effect in the auto-correlograms of the EOMs' weights (the larger the better). The EOMs that allow recovering 60% of the input variance are indicated in red and connected for better visibility.

In the text

Fig. 4 Portfolio of octopus during the hot (left − red) and the cold (right − blue) seasons over the entire study period. Each plot is made of an EOM with the temporal evolution of its weights, smoothed by a loess method with a span parameter value of 0.5 and associated 95%-confidence interval. Pie chart indicate the percentage of variance explained by the EOM (arrows represent the cumulative percentage of variance explained by the current EOM and the previous ones and the percentage of variance explained by the portfolio). The EOMs that were not selected for the portfolio are represented in grey.

In the text

Fig. 5 Portfolio of SST during the hot (left − red) and the cold (right − blue) seasons over the entire study period. Each plot is made of an EOM with its weights (smoothed by nonparametric loess with a span parameter value of 0.5 and associated 95%-confidence interval) and a pie indicating the percentage of variance explained by the EOM (arrows represent the cumulative percentage of variance explained by the current EOM and the previous ones).

In the text

Fig. 6 Analysis of the temporal dynamics of the principal spatial patterns (EOM 1, 3, 4, 5) and of the abundance for the cold season. Diagonal panels: auto-correlograms; Lower triangle panels: cross-correlograms. For all panels, the x-axis represents time lags in years (from 0 to 30 years) and the y-axis represents correlations. The horizontal black dashed line is at 1. For off diagonal panels, the horizontal dotted lines represent the correlations. The panels of the last line concerns the cross-correlograms of EOM weights with abundance and the auto correlogram for abundance (last panel). The horizontal red dashed lines represents the correlation between EOM weights and abundance. The linear model of coregionalisation is represented as a continuous line (nugget effect + slope). The isolated panel represents the correlation between abundance and EOM weights as a function of the nugget effect in the model.

In the text

	Fig. 7 Cold season. Time series of the abundance (black) and temporal weights of EOM3 (red).
In the text

	Fig. 8 Climatology for the cold season. This corresponds to the mean of the EOMs selected in the portfolio weighted by their mean weights. The scale is consistent with the input octopus densities (i.e. nb of individual per m2).
In the text

Fig. 9 Portfolio for the cold season. Distribution of the EOMs' weights as a function of the EOM numbers for March and December surveys. The weights for the surveys carried out in March (N = 10) are represented in black. The weights for surveys performed in December (N = 10) are given in red. The average values are superimposed. For the two instances where a December survey is followed by a March surveys three months later, the evolutions of the weights are represented by doted segments.

In the text

Fig. 10 Portfolio of the cold season. Correlation between the weights of EOM3 and the upwelling index. Left: lagged correlations between EOM3 weights in March and the upwelling index 6 months earlier. Right: boxplot of upwelling indices for all months except September and September (i.e. 6 months before March). In red the values for the nine years when a survey was performed in March with concomitant upwelling observations.

In the text

	Fig. 11 Left: Empirical Orthogonal Map EOM2 of octopus for the cold season. This EOM was not selected for the portfolio. Right: temporal evolution of its weights smoothed by a loess method with a span parameter value of 0.5 and associated 95%-confidence interval.
In the text

	Fig. 12 Comparison between the MAF and the EOM arrangements for the hot (left) and the cold (right) season for octopus in Mauritanian waters. The green connections concern the Q first EOMs explaining 60% of total variance.
In the text