Issue 
Aquat. Living Resour.
Volume 36, 2023



Article Number  15  
Number of page(s)  14  
DOI  https://doi.org/10.1051/alr/2023008  
Published online  16 May 2023 
Research Article
Portfolio of distribution maps for Octopus vulgaris off Mauritania
^{1}
IMROP, Nouadhibou, Mauritania
^{2}
Ministry of Fisheries and Maritime Economy, Nouakchott, Mauritania
^{3}
MARBEC, IRD, Univ Montpellier, Ifremer, CNRS, INRAE, Sète, France
^{*} Corresponding author: dedah.ahmedbabou@ird.fr
Handling Editor: Verena Trenkel
Received:
19
October
2022
Accepted:
30
March
2023
This study introduces the concept of portfolios of distribution maps, which consist of the reduced set of empirical orthogonal maps that best explain spatial biomass distributions of a given species over time. The approach is demonstrated for the distributions of common octopus (Octupus vulgaris) off Mauritania over the last thirty years. The maps in the portfolio are the subset of empirical orthogonal maps that allowed to recover 60% of the spatiotemporal biomass distribution variance and whose temporal weights were significantly correlated with abundance. For octopus during the hot season, one single map explained half of the overall variance of the distribution data, while during the cold season, the portfolio of octopus distribution maps consisted of four maps, with the temporal weights of the second map being negatively correlated with upwelling intensity six months before. The size of each portfolio represents the number of distinct spatial patterns describing octopus spatial distributions. Assuming that specific but hidden processes explain each biomass spatial distribution of the portfolio, the size of a map portfolio might be interpreted as a proxy for system resilience. A small portfolio could reflect systems that are more fragile.
Key words: Spatiotemporal distribution data / Maps portfolio / temporal dynamics / Octopus vulgaris
© D. AhmedBabou et al., Published by EDP Sciences 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Ecological monitoring and information systems are increasing worldwide. They are key to providing longterm 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 timedependent and the other, possibly low dimensional, spacedependent. The approach is based on minmax 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 smallscale 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ándezLó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 (NorthSouth) 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 NorthSouth 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 coldwarm transition season from June to July, a hot season from August to October, and a hotcold 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 minmax 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}):
In the spatiotemporal 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:
This second PCA is based on the eigendecomposition of the matrix of variogram and crossvariogram 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 spatiotemporal 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}).
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 pseudoalgorithm 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, cokriging 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 Qfirst 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 crosscorrelograms to identify EOMs whose weights vary smoothly over time (autocorrelogram of weights) and in a consistent manner with regards to the temporal variation of abundance (crosscorrelogram between weights and abundances). For this, a linear model of coregionalization (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 autocorrelograms). 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 crosscorrelograms between weights and abundance. In linear models of coregionalization, the (absolute) value of the slope of the crosscorrelogram must be smaller than the square root of the product of the slopes of the two autocorrelograms (Chilés and Delfiner, 2009). Thus, the third criterion is that the ratio between the slope of the crosscorrelograms 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 crosscorrelograms 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:
where is the mean temporal weight of the P selected EOMs such that:
and where is the mean abundance:
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 northsouth gradient, with a highdensity 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 (SeptemberOctober). 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 xaxis corresponds to the correlation between EOM weights and abundance (the closer to +/− 1, the better). The yaxis represents the relative values of the slopes of the crosscorrelograms 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 autocorrelograms 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 spatiotemporal 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 acrossshelf 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 longterm 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 crosscorrelograms, their fluctuations at short temporal scales were not correlated (Fig. 6). A noticeable exception concerned the short term negative crosscorrelation 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 shortterm joint fluctuations (e.g. the average variations between two consecutive cold seasons). The slope characterized longterm (multiannual) joint evolutions. In this context, a change in the overall correlation can be due to either short or longterm processes, or both. However, by construction, the EOMs of the portfolio have been selected because they shared longterm 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 longterm variations, their shortterm fluctuations were also correlated. The crosscorrelogram 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 autocorrelation (selfkrigeability; 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 spatiotemporal 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 pvalue = 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: autocorrelograms; Lower triangle panels: crosscorrelograms. For all panels, the xaxis represents time lags in years (from 0 to 30 years) and the yaxis 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 crosscorrelograms 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 m^{2}). 
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 crosscorrelations 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 twentyeight 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 spatiotemporal decomposition of the time series of thirtyfour 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 northsouth 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 crossshelf 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 longterm trends with that of abundance, the shortterm (i.e. yeartoyear) 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 prespawning adults of 910 months. They come from the spawning that took place around May the year before. Given the duration of larvae and paralarvae developments (Otero et al., 2007, Domain et al., 2000), the end of the pelagic phase for these individuals falls around AugustSeptember. 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
 Balguerías E, HernándezGonzález C, PeralesRaya C. 2002. On the identity of Octopus vulgaris Cuvier, 1797 stocks in the Saharan Bank (Northwest Africa) and their spatiotemporal variations in abundance in relation to some environmental factors. Bull Mar Sci 71: 147–163. [Google Scholar]
 Bez N, Renard D, AhmedBabou D. 2022. Empirical Orthogonal Maps (EOM) and distance between empirical spatial distributions. Application to Mauritanian octopus distribution over the period 1987–2017. Math Geosci 55: 113–128. [Google Scholar]
 Boyle PR, Boletzky SV. 1996. Cephalopod populations: Definition and dynamics. Philos Trans Royal Soc Lond Ser B 351: 985–1002. [Google Scholar]
 Caddy JF. 1983. The cephalopods: Factors relevant to their population dynamics and to the assessment and management of stocks. Adv Assess World Cephalopod Resour 231: 416–449. [Google Scholar]
 Caverivière A, Domain F, Diallo A. 1999. Observations on the influence of temperature on the length of embryonic development in Octopus vulgaris (Senegal). Aquat Liv Resour 12: 151–154. [CrossRef] [Google Scholar]
 Chilés JP, Delfiner P. 2009. Geostatistics: Modeling Spatial Uncertainty. John Wiley & Sons. [Google Scholar]
 Domain F, Jouffre D, Caverivière A. 2000. Growth of Octopus vulgaris from tagging in Senegalese waters. J Mar Biol Assoc UK 80: 699–705. [Google Scholar]
 Demarcq H, Faure V. 2000. Coastal upwelling and associated retention indices derived from satellite SST. Application to Octopus vulgaris recruitment. Oceanolog Acta 23: 391–408. [Google Scholar]
 Dobrovine B, Ould Mohamed Mahfoud M, Ould Sidina D. 1991. La ZEE mauritanienne et son environnement géographique géomorphologique et hydroclimatique. Bull Sci CNROP 24. https://aquadocs.org/handle/1834/518. [Google Scholar]
 Faure V, Inejih C, Demarcq H, Cury P. 2008. The importance of retention processes in upwelling areas for recruitment of Octopus vulgaris: The example of the Arguin Bank (Mauritania). Fish Oceanogr 9: 343–355. [Google Scholar]
 Gascuel D, Labrosse P, Meissa B, Taleb Sidi MO, Guénette S. 2007. Decline of demersal resources in NorthWest Africa: An analysis of Mauritanian trawlsurvey data over the past 25 years. Afr J Mar Sci 29: 331–345. [CrossRef] [Google Scholar]
 Goulard M, Voltz M. 1992. Linear coregionalization model: Tools for estimation and choice of crossvariogram matrix. Math Geol 24: 269–286. [Google Scholar]
 Guerra A. 1979. Fitting a von Bertalanffy expression to Octopus vulgaris growth. Investigación Pesquera 43: 319–326. [Google Scholar]
 Hatanaka H. 1979. Studies on the fisheries biology of common octopus off the northwest coast of Africa. Bull Far Seas Fish Res Lab 17: 13–124. [Google Scholar]
 HernándezLópez JL, CastroHernández JJ, HernándezGarcía V. 2001. Age determined from the daily deposition of concentric rings on common octopus (Octopus vulgaris) beaks. Fish Bull 99: 679–684. [Google Scholar]
 Khallahi B, Taleb H, Barham CB, Habibe BM, Kane EA, Bouzouma ME (Eds). 2020. Aménagement des ressources halieutiques et gestion de la biodiversité au service du développement durable. Rapport du Neuvième Groupe de Travail de l'IMROP, Nouadhibou. Mauritanie 11–14 février 2019. 246p. [Google Scholar]
 Lorenz EN. 1956. Empirical orthogonal functions and statistical weather prediction. Massachusetts Institute of Technology, Department of Meteorology Cambridge. Statistical forecasting project, Scientific report 1. [Google Scholar]
 Mangold K. 1983. Octopus vulgaris. In Boyle, PR (ed.), Cephalopod Life Cycles (I): Species Accounts. Orlando: Academic Press, 335–363. [Google Scholar]
 Mangold K, von Boletzky S. 1973. New data on reproductive biology and growth of Octopus vulgaris. Mar Biol 19: 7–12. [Google Scholar]
 Meissa B, Gascuel D, Rivot E. 2013. Assessing stocks in datapoor African fisheries: A case study on the white grouper Epinephelus aeneus of Mauritania. Afr J Mar Sci 35: 253–267. [CrossRef] [Google Scholar]
 Meissa B, Gascuel D. 2015. Overfishing of marine resources: Some lessons from the assessment of demersal stocks off Mauritania. ICES J Mar Sci 72: 414–427. [CrossRef] [Google Scholar]
 MINES ParisTech / ARMINES. 2022. RGeostats: The Geostatistical R Package. Version: [12.0.1] Free download from: http://cg.ensmp.fr/rgeostats [Google Scholar]
 National Oceanic and Atmospheric Administration. 2020. Global Coral Bleaching Monitoring, 5km, V.3.1, Monthly, 1985Present. https://coastwatch.pfeg.noaa.gov/erddap/griddap/NOAA_DHW_monthly.html?sea_surface_temperature%5B(20210615 T23: 00:00Z)%5D%5B(89.975):(89.975)%5D%5B(179.975):(179.975)%5D&.draw=surface&.vars=longitude%7Clatitude%7Csea_surface_temperature&.colorBar=%7C%7C%7C%7C%7C&.bgColor=0xffccccff [Google Scholar]
 Otero J, González ÁF, Sieiro MP, Guerra Á. 2007. Reproductive cycle and energy allocation of Octopus vulgaris in Galician waters, NE Atlantic. Fish Res 85: 122–129. [CrossRef] [Google Scholar]
 Otero J, ÁlvarezSalgado XA, González ÁF, Miranda A, Groom SB, Cabanas JM, Casas G, Wheatley B, Guerra Á. 2008. Bottomup control of common octopus Octopus vulgaris in the Galician upwelling system, northeast Atlantic Ocean. Mar Ecol Progr Ser 362: 181–192. [Google Scholar]
 Pease N. 1973. Fisheries of Mauritania, 1971. https://repository.library.noaa.gov/view/noaa/30110 [Google Scholar]
 Petitgas P, Renard D, Desassis N, Huret M, Romagnan JB, Doray M, Woillez M, Rivoirard J. 2020. Analysing temporal variability in spatial distributions using minmax autocorrelation factors: sardine eggs in the Bay of Biscay. Mathe Geosci 52: 337–354. [Google Scholar]
 Sanchez FJ, Iglesias J, Moxica C, Otero JJ. 1998. Growth of octopus (Octopus vulgaris) males and females under culture conditions. International Council for the Exploration of the Sea (ICES), CM 1998/M: 47. [Google Scholar]
 Solow AR. 1994. Detecting change in the composition of a multispecies community. Biometrics 50: 556–565. [CrossRef] [PubMed] [Google Scholar]
 Société Mauritanienne de Commercialisation Des Poissons. SMCPSEM. 2020. Bulletin de statistiques 2019 Version 42. http://dtf.gov.mr/images/2019/SMCP2019.pdf. [Google Scholar]
 Switzer P, Green AA. 1984. Min/max autocorrelation factors for multivariate spatial imagery. Technical report 6, Department of Statistics, Stanford University, CA. [Google Scholar]
 Villanueva R. 1995. Experimental rearing and growth of planktonic Octopus vulgaris from hatching to settlement. Can J Fish Aquat Sci 52: 2639–2650. [CrossRef] [Google Scholar]
 Wikle CK, ZammitMangion A, Cressie NAC. 2019. Spatiotemporal statistics with R. Chapman and Hall/CRC. [Google Scholar]
 Woillez M, Rivoirard J, Petitgas P. 2009. Using min/max autocorrelation factors of surveybased indicators to follow the evolution of fish stocks in time. Aquat Liv Resour 22: 193–200. [CrossRef] [EDP Sciences] [Google Scholar]
Cite this article as: AhmedBabou D, Demarcq H, MEISSA B, Bez N. 2023. Portfolio of distribution maps for Octopus vulgaris off Mauritania. Aquat. Living Resour. 36: 15
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 Qfirst 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 xaxis corresponds to the correlation between EOM weights and abundance (the closer to +/− 1, the better). The yaxis represents the relative values of the slopes of the crosscorrelograms 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 autocorrelograms 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: autocorrelograms; Lower triangle panels: crosscorrelograms. For all panels, the xaxis represents time lags in years (from 0 to 30 years) and the yaxis 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 crosscorrelograms 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 m^{2}). 

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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.