Automatic fish detection in underwater videos by a deep neural network-based hybrid motion learning systemSalman, Ahmad; Siddiqui, Shoaib Ahmad; Shafait, Faisal; Mian, Ajmal; Shortis, Mark R; Khurshid, Khawar; Ulges, Adrian; Schwanecke, Ulrich
doi: 10.1093/icesjms/fsz025pmid: N/A
Abstract It is interesting to develop effective fish sampling techniques using underwater videos and image processing to automatically estimate and consequently monitor the fish biomass and assemblage in water bodies. Such approaches should be robust against substantial variations in scenes due to poor luminosity, orientation of fish, seabed structures, movement of aquatic plants in the background and image diversity in the shape and texture among fish of different species. Keeping this challenge in mind, we propose a unified approach to detect freely moving fish in unconstrained underwater environments using a Region-Based Convolutional Neural Network, a state-of-the-art machine learning technique used to solve generic object detection and localization problems. To train the neural network, we employ a novel approach to utilize motion information of fish in videos via background subtraction and optical flow, and subsequently combine the outcomes with the raw image to generate fish-dependent candidate regions. We use two benchmark datasets extracted from a large Fish4Knowledge underwater video repository, Complex Scenes dataset and the LifeCLEF 2015 fish dataset to validate the effectiveness of our hybrid approach. We achieve a detection accuracy (F-Score) of 87.44% and 80.02% respectively on these datasets, which advocate the utilization of our approach for fish detection task. Introduction Monitoring the effect of preventive and recovery measures requires the estimation of fish biomass, and abundances by sampling their populations in water bodies like lakes, rivers and oceans on a regular basis (Jennings and Kaiser, 1998). This requires observation of the interaction of different fish species with changing environmental conditions. This is an essential process, especially in those regions of the world where certain species of fish are either threatened or at the risk of extinction due to habitat loss and modification, industrial pollution, deforestation, climate change, and commercial overfishing (Tanzer et al., 2015). There is a well-established and increasing interest in using non-destructive fish sampling techniques by marine biologists and conservationists (McLaren et al., 2015). Underwater video-based fish detection approaches have been used to achieve non-destructive and repeated sampling for many years (Harvey and Shortis, 1995; Shortis et al., 2009). Manual processing of underwater videos is labour intensive, time consuming, expensive and prone to fatigue errors. In contrast, automatic processing of the underwater videos for fish species classification and biomass measurement is an attractive alternative. However, high variability in underwater environments due to changes in lighting conditions, clarity of water, and background confusion due to vibrant seabed structure pose great challenges towards automatic detection of fish. These factors result in a compromise on accuracy, which supports the continuing practice of less cost effective and cumbersome manual sampling and tagging of fish. In general, automatic fish sampling involves the following three major tasks: (i) Fish detection, which discriminates fish from non-fish objects in underwater videos. Non-fish objects include coral reefs, aquatic plants, sea grass beds, sessile invertebrates such as sponges, gorgonians, ascidians, and general background. (ii) Fish species classification, which identifies the species of each detected fish from the predetermined pool of different species (Siddiqui et al., 2017). (iii) Fish biomass measurement, using length to biomass regression methods (Froese, 2006). This article addresses the first task and the interested reader is referred to the literature for details of the following two steps in the overall process. Various approaches have been followed for fish detection and consequently their assemblage estimation using image and video processing algorithms. Broadly speaking, these approaches can be divided into two categories based on the medium available for sampling, namely constrained and unconstrained sampling. In the former case, early attempts were made that involved detection of fish using information of their shape and colour (Strachan and Kell, 1995) or 3D modelling of fish to acquire features like height, width, or thickness (Storbeck and Daan, 2001). Harvey and Shortis (1995) presented an approach to acquire underwater images of fish under controlled conditions. This was achieved by making fish swim through a chamber with controlled illumination. Unconstrained underwater fish detection and classification does not assume any specific environmental conditions and, therefore, faces difficulty in achieving the required accuracy due to high variations in the aforementioned conditions. To address this problem, Spampinato et al. (2008) presented an image processing based method for fish detection and counting by capturing the texture pattern of fish in the natural underwater environment. They were able to achieve an average accuracy of about 84% on five underwater videos. In the past, several attempts have been made to solve the same problem in underwater videos using machine learning. Principal component analysis (Turk and Pentland, 1991), linear discriminant analysis (Mika et al., 1999), and sparse representation-based detection (Hsiao et al., 2014) presented some ways to capture fish-dependent features through mathematical modelling, which assumed independence of modelled fish with surrounding environments in videos. In other words, information like fish colour, texture, and shape was extracted with the prior assumption that foreground fish instance was easily distinguishable from the background. In reality, it is challenging to differentiate fish within underwater video/images due to camouflage with the background, poor visibility, and loss of contrast as a result of light attenuation through the water medium, low light, and water turbidity. In pursuit of suppressing the effects of environmental variability, Kernel Descriptors in Kernel density estimation (KDE) approach with colour information for background pixel modelling in images were used by Sheikh and Shah (2005). In contrast, texture-dependent features computed via local binary patterns for background modelling was proposed in Yao and Odobez (2007). Background modelling is a popular technique to segment moving foreground objects from the background in video sequences. An approach using motion-based fish detection in videos was presented by Hsiao et al. (2014). This method implements background subtraction by modelling background pixels in the video frames using Gaussian mixture models (GMMs). Although training the GMM, it is assumed that subsequent frames of video lack fish instances. Motion is detected in the video frames (apparently from fish) when a certain region of the frame does not fit into the trained background model. This approach produces fish detection results with an average success rate of 83.99% on several underwater videos collected near southern Taiwan. A similar scheme was proposed on covariance modelling of background and foreground (fish instances) in the video frames using colour and texture features of fish (Palazzo and Murabito, 2014). Using a dataset of four underwater videos with a high variation in luminosity, strong background movements, dynamic textures, and rich background, they were able to achieve an average detection accuracy of 78.01%. Presently, GMM- and KDE-based fish detection approaches are considered state-of-the-art (Spampinato et al., 2014). We will compare the performance of various state-of-the-art techniques with our proposed approach in a later section. All of the above-mentioned machine learning and feature extraction approaches fall into the category of shallow learning architectures (Bengio, 2009). These techniques are unable to accurately model the complexity of fish-dependent features in the presence of highly variable and diverse environments, and therefore these video or image-based fish detection techniques exhibit low performance in real-world scenarios (Siddiqui et al., 2017). In the last decade, deep learning has been at the centre of attention for many researchers developing detection and classification algorithms in computer vision. Marked by their ability to extract and model highly nonlinear data, deep architectures have been utilized in numerous tasks related to computer vision, including facial recognition, speech processing, generic object detection, and classification in video and still images producing state-of-the-art results (Lin et al., 2015; Ren et al., 2017). In realizing deep architectures, multilayer deep neural networks are among the most successful schemes capable of extracting task-dependent features in the presence of variability in the images. Most commonly used variants of deep neural networks include deep convolutional neural networks (CNNs) which are parametric neural network models capable of extracting task-specific features and are widely used in computer vision problems like object recognition in images and facial recognition (LeCun et al., 2015). Deep learning is being used lately to solve fish-related tasks (Moniruzzaman et al., 2017). An important work using CNN was proposed by Sung et al., (2017) to detect fish in underwater imagery with 65.2% average accuracy on a dataset containing 93 images having fish instances. The system was trained on raw fish images to capture colour and texture information for localizing and detecting fish instances in the images. In a similar work, deep region-based CNN (R-CNN) were used for the abundance estimation of fish from 4909 underwater images recorded in the coast of Southeast Queensland, Australia. In this work, an accuracy of 82.4% was reported using the R-CNN system tuned for locating and detecting fish instances in an image with a unified network framework. Despite the high accuracy achieved by the deep learning based fish species classification, the task of vision-based automatic fish detection in unconstrained underwater videos is still under extensive investigation as most of the previous attempts reported results on relatively small datasets with a limited variety in the surrounding environment. Therefore, it is important to judge the robustness and effectiveness of any system in a large dataset with a high degree of environmental variation. In this article, we address fish detection in highly diverse and unconstrained underwater environments and propose a hybrid system based on explicit motion-based feature extraction followed by a final detection phase using deep CNNs. In the first step, we use background subtraction by modelling still pixels of the video frames using GMMs. These models represent pixels related to a range of coral reefs, seabed features, and aquatic plants. Foreground objects are segmented from the background based on the motion in the scene that does not match the background model. To enhance the quality of the extracted features in each video frame, we concatenate the GMM candidate output blobs with the moving candidates generated by optical flow, a well-established approach used for motion detection in videos (Brox et al., 2004). However, due to poor image quality, noise and background confusion, the detection remains far from perfect. To address this problem, we tune the parameters of GMM and optical flow systems to generate high recall by trying various values of the number of Gaussian distributions, initial variance, blob size and sensitivity in case of GMM, as well as pyramid size, number of pyramid layers, and window size in case of optical flow. The details of these parameters are given in Zivkovic and Heijden (2006) for GMM and in Beauchemin and Barron (1995) for optical flow. Specifically, in this step, all entities that exhibit even a slight movement are detected as fish. In the second step, we discriminate all the candidate regions in the video frames as fish and non-fish entities using a CNN architecture arranged in a hierarchical fashion to fine tune the detection system. Our CNN is trained using a supervised training style in which the GMM and optical flow blobs acts as the input while ground truth blobs (given in the training data) acts as the desired output. We worked on two different datasets; the Fish4Knowledge Complex Scenes Dataset, where the aim is fish detection with videos arranged into seven different categories based on the variation in the underwater environment; and the LifeCLEF 2015 (LCF-15) dataset, which is also designed for the detection of freely swimming fish in video sequences. These datasets contain marine scenes and species; unfortunately, there is no public domain benchmark datasets available containing underwater recordings in fresh water bodies. The contribution of this work is to overcome the main challenge faced by the conventional motion detection and image classification approaches using deep learning. These deep learning modules are trained to select the relevant information from the data and minimize confusion which contributes to false alarms or missed detections. This approach improves the detection and classification accuracies especially in the data marked by high environmental variability like unconstrained underwater videos of fish. Our novelty lies in the proposed hybrid setup to mine the relevant motion information content by pooling the information generated by GMM and optical flow and refining the outcome by deep CNNs. Our approach is capable of detecting fish in the video in its stationary or moving state with region-based feature localization. This equips our detection system with motion-influenced temporal information that is not available otherwise, in order to enhance detection performance in cases where fish is occluded or camouflaged in the background. Material and methods Dataset We use two benchmark datasets in our study, both of which are specially designed to provide a resource for testing algorithms for detection and classification of fish in images and video sequences and have been used for benchmarking a number of approaches. The first dataset is used for the fish detection task and is a collection of 17 videos under different environmental conditions (http://f4k.dieei.unict.it/datasets/bkg_modeling/). The second dataset is taken from the LCF-15 fish task (http://www.imageclef.org/lifeclef/2015/fish). This dataset contains 93 underwater videos comprising 15 different fish species. Both datasets are derived from a very large fish database called Fish4Knowledge (Fisher et al., 2016). With over 700 000 underwater videos in unconstrained conditions, the Fish4Knowledge dataset has been collected over a period of 5 years to monitor the marine ecosystem of coral reefs around Taiwan. This region is home to one of the largest fish biodiversity environments in the world with more than 3000 fish species. The first dataset, dubbed FCS (Fish4Knowledge with Complex Scenes) hereinafter, comprises seven sets of selected videos recorded in typical underwater conditions addressing complex variability in the scenes. Thereby, the environmental variations provide a major challenge for fish identification and are categorized as follows: Blurred, comprising three low contrast, blurred videos. Complex background, composed of three videos with rich seabed structures that provide a high degree of background confusion. Crowded, in which three videos with a high density of moving fish in each video frame imposes specific challenges for fish detection techniques, especially when it comes to high recall and precision in the presence of occluding objects. Dynamic background, in which two videos are provided with rich textures of coral reef background and moving plants. Luminosity variation composed of two videos involving sudden luminosity changes due to surface wave action. This phenomenon can induce false positives in detection due to moving light beams. Camouflage foreground, two videos are chosen, addressing the challenge of detecting fish camouflaged in the presence of textured and colourful background. Hybrid, in which two videos are selected to show a combination of all the above-mentioned conditions of variability. Table 1 summarizes the technical information regarding both datasets used in this article. For the FCS dataset, complexity is specifically depicted for all seven environmental conditions. The LCF-15 dataset is used to detect fish instances in the video i.e. to count all the fish in the video regardless of their species. Of the 93 videos given in LCF-15, 20 are used for training the computer vision or machine learning modules, while the remaining 73 videos are set aside for testing/validating the developed algorithms. In total, there are 9000 annotated fish instances available in the LCF-15 training set, and 13 493 annotated instances for the test videos. All these videos are manually annotated by experts. Apart from videos, there are 20 000 labelled still images in LCF-15, where each image comprises of a single fish. These images can also be used to supplement the training set if required. Thus, in total there are 42 493 labelled fish instances in videos and still images in the LCF-15 dataset. The FCS dataset is also designed and used for the fish detection task. Therefore, ground truth is available for all moving fish, frame by frame in each video. There are a total of 1328 fish annotations available for the FCS dataset. Figure 1 shows some video frames extracted from FCS and LCF-15 datasets exhibiting the variation in the surrounding environment, fish patterns, shape, size, and image quality. Figure 1. Open in new tabDownload slide Sample images to illustrate the high variation in underwater environment. The first two rows depict seven categories of the FCS dataset from left to right top to bottom being Blurred, Complex background, Dynamic background, Crowded, Luminosity variation, Camouflage foreground, and Hybrid. The last row shows an excerpt from different videos of the LCF-15 dataset. Table 1. Information about LCF-15 and FCS fish datasets. Dataset . No. of videos . Format . Resolution . Frames/sec . No. of labelled fish instances . LCF-15 93 FLV 640 × 480, 320 × 240 24 42 493 FCS 17 FLV 640 × 480, 320 × 240 24.5 1 328 Dataset . No. of videos . Format . Resolution . Frames/sec . No. of labelled fish instances . LCF-15 93 FLV 640 × 480, 320 × 240 24 42 493 FCS 17 FLV 640 × 480, 320 × 240 24.5 1 328 Open in new tab Table 1. Information about LCF-15 and FCS fish datasets. Dataset . No. of videos . Format . Resolution . Frames/sec . No. of labelled fish instances . LCF-15 93 FLV 640 × 480, 320 × 240 24 42 493 FCS 17 FLV 640 × 480, 320 × 240 24.5 1 328 Dataset . No. of videos . Format . Resolution . Frames/sec . No. of labelled fish instances . LCF-15 93 FLV 640 × 480, 320 × 240 24 42 493 FCS 17 FLV 640 × 480, 320 × 240 24.5 1 328 Open in new tab Proposed algorithm To perform fish detection, we propose a hybrid system based on the initial motion-based feature extraction from videos using GMM and optical flow candidates. These feature images are combined with raw greyscale images and fed to the CNN system to mark final detected fish. Therefore, our proposed hybrid fish detection system is made up of three components i.e. GMM, optical flow and a CNN. Gaussian mixture modelling In machine learning, GMM is an unsupervised generative modelling technique to learn first and second order statistical estimates of input data features (Stauffer and Grimson, 1999; Zivkovic and Heijden, 2006). This approach and its variants are frequently used in computer vision and speech processing tasks. GMM represents a probability density function Pxt at time t of data x as a weighted sum of multiple individual normal distributions ηxi for pixel i . Thereby, each density is characterized by the mean and covariance of the represented data. Using a combination of individual Gaussian densities, any probability distribution can be estimated with arbitrary precision (Reynolds and Rose, 1995). In our case, each pixel value with a fixed location in the video frame acts as a feature. Multiple such values from successive frames are combined to form a feature vector. As elaborated in Figure 2, we end up with a total number of feature vectors that equals the total number of pixels in a video frame. The GMM requires a certain amount of training data to effectively estimate the mean and covariance of an object class. For fish detection in videos, there are two classes i.e. background and foreground. Ideally, the background in underwater videos should cover everything in the frame but moving fish. For example, seabed structure, coral reef, aquatic plants, and wave action causing variation in light intensity are categorized as background. Freely moving fish, on the other hand, constitute as foreground. The GMM is used to learn the background features in a statistical model using mean and covariance values of pixel data and separate them from the foreground pixel distribution. In other words, any random and sudden change in the pixel value of a frame causes a mismatch with the background model of that pixel and hence, a motion is assumed to be detected. The statistical pattern of foreground (fish in our case) movement is usually different from the pattern of fixed objects like seabed structures, coral reefs and also objects with confined movement like to and fro motion of plants and refracted light rays from the surface. The outputs of the GMM are the candidate blobs marked by bounding boxes localizing the moving objects in a frame (see Figure 2). Figure 2. Open in new tabDownload slide (A) Illustration of background subtraction and foreground segmentation using GMM which detects any change in the foreground by matching it with the background model. (B) Motion detection in an optical flow setup to estimate the direction of moving objects in two dimensions (x, y) for consecutive frames in time t of a video sequence. The video frames that are used to train the GMM should not include any fish instance but only the background. However, it is very difficult to capture such videos in a natural environment as fish can appear in any number of frames. When a GMM is trained on videos that do not have pure background but also some fish, the fish will also be modelled as background resulting in misdetections in the upcoming test frames. Optical flow To compensate for this shortcoming of GMM, we additionally extracted optical flow features which are purely generated by motion occurring in the underwater videos (see Figure 2). Optical flow is a 2D motion vector in the video footage caused by the 3D motion of the displayed objects (Warren and Strelow, 1985). There are various methods to estimate optical flow. We opted for a simple yet effective method where motion is detected between two successive video images taken at times t and t+Δt at every position using Taylor series approximation with partial derivatives with respect to spatial and temporal coordinates (Beauchemin and Barron, 1995). A region of interest (ROI) in a video frame at time t and coordinates x, y can be represented in terms of intensity as Ix,y,t . After any motion in the next frame, the intensity becomes Ix+Δx,y+Δy,t+Δt where the notation Δ represents the change in coordinates and time. Based on the motion constraint, optical flow can be determined as described in, for example, Beauchemin and Barron (1995). Optical flow depends on the analysis of the consecutive frames to estimate the difference in the intensity vector of a pixel at a particular location. However, such an analysis is also prone to false motion detection apart from fish when applied to a dynamic background with moving aquatic plants and abrupt luminosity variation due to disturbance at the water surface. The parameters of the GMM and optical flow algorithm are chosen such that even the smallest movements are detected. In other words, the sensitivity of the algorithms is maximized, producing a high rate of false alarm in addition to detecting fish instances leading to high recall rates. In the next step, the precision of the system is further increased by fine-tuning and refining regions in the frames to localize moving fish. This requires a robust detector to categorize fish motion in complex and variable environments. We propose the use of a R-CNNs (hereinafter referred to as R-CNN) trained on images, created by combining candidate regions generated by the GMM and optical flow together with the original greyscale images in a supervised learning setup. Region-based convolutional neural network A deep CNN is a nonlinear parametric neural network model capable of extracting and learning complex yet abstract features of the input data. Variations in the lighting condition, size, shape, and orientation of the fish, poor image quality and significant noise are the factors that introduce nonlinearity into the data (Bengio, 2009). Since all of these challenges are encountered in the videos recorded in an unconstrained underwater scenario, it is difficult for conventional machine learning algorithms to model data features in the presence of such nonlinearity. However deep neural architectures, especially CNNs, learn to extract invariant and unique features of the objects of interest in data when properly trained with a sufficient amount of labelled data (LeCun et al., 2004; Simonyan and Zisserman, 2014). The deep architecture exemplified by the R-CNN employed in our study is a hierarchical parametric model composed of two modules. The first module is a generic deep CNN trained for generic object recognition on a very large dataset called ImageNet (Deng et al., 2009). Smaller than the first module CNN, the second module is another CNN, which acts as the object detector and called region proposal network (RPN) (Ren et al., 2017). It selects candidate regions in the feature space of the input image in which a motion is likely to have occurred. The entire system is used for detecting moving objects as depicted in Figure 3. The first module utilizes the concept of transfer learning (Siddiqui et al., 2017). It learns characteristic feature representation of the object of interest in the input image in order to recognize and classify the objects in the imagery. In transfer learning, a CNN pretrained on totally different, yet relevant dataset, is utilized as a generic feature extractor for the dataset of interest. In our case, the CNN was trained on the vast ImageNet dataset that contains 1.2 million images of a very large and diverse number of objects. This dataset is not related or designed for fish species recognition or fish detection in underwater videos. However, it provides a high degree of variability to detect generic objects with different backgrounds in input images based on their texture, size and shape features. Once the network is trained, it can be applied to a different dataset, in our case on underwater video imagery of fish, as a feature extractor. Transfer learning is suitable for the applications where a large amount of training data is not available to train the deep CNNs (Siddiqui et al., 2017). This is exactly the problem in the current underwater datasets. Training on such relatively small datasets (see Table 1) overfits a deep CNN to generate better performance on training dataset and fails on previously unseen test datasets. In other words, the training dataset is so small that the CNN is able to memorize it and produce good results only on the training dataset. We utilize a deep CNN known as ResNet-152 as the pretrained model (He et al., 2016). The parameters of this network are further refined by including examples of our fish dataset video imagery in training. This network is composed of an input layer, various hidden layers and an output layer to process an input image to obtain its output feature representation (LeCun et al., 2004). Starting from an input layer that represents the pixels of an image, the hidden layers are interconnected by a set of weights that are tuneable as a result of the training procedure. Thereby, each hidden layer represents a higher-level form of feature representation. There are several types of hidden layers used in our network, e.g. a convolution layer that performs the mathematical operation of convolution between image pixels (values of the input layer) or feature values (values of the hidden layers) with the weight vectors. Convolution is generally used in image processing for noise reduction or detecting features such as edges or other patterns (LeCun et al., 1989). In a CNN, convolution is followed by a nonlinear activation layer to induce nonlinearity in the feature vectors. There are several types of nonlinear functions, e.g. ReLUs (rectified linear units), Sigmoid and Hyperbolic Tangent (LeCun et al., 2004; Simonyan and Zisserman, 2014; He et al., 2016). The choice of the nonlinear function depends on the data distribution and nonlinearity of the input data. Due to the saturating regions of the Tangent Hyperbolic and Sigmoid function, the ReLU function is the defacto-standard in the latest state-of-the-art models. Max pooling and average pooling layers sift out the most prominent values from the output of nonlinearity inducing layers based on maxima or an averaging operations to reduce the dimension of feature vectors and retain useful information while discarding the redundancy. The final layer is the output layer which usually is a classification layer with output nodes equal to the number of desired classes for a given dataset. Each output node produces a score or probability for the associated class. The predicted label is then matched with the ground truth label to calculate accuracy. Figure 3. Open in new tabDownload slide The proposed hybrid system, where ResNet-152 CNN is trained on images that are created by combining the motion-influenced outputs of GMM and optical flow algorithms with raw greyscale video images. This is analogous to three-channel RGB image. ResNet-152 is a modular deep CNN with various hidden layers. The architecture is designed to process images of size 224 × 224 given the fact that this resolution is enough to extract useful features within reasonable computational time. Thus, after applying five pooling layers, the feature map size shrinks to 7 × 7 which can be processed by fully-connected layer of 1000 label prediction nodes, since ResNet-152 was designed to train on a subset of ImageNet dataset with 1000 classes. The complete architecture details of ResNet-152 can be found in He et al. (2016). The arrangement of the above-mentioned layers in this architecture is experimentally determined to yield greater success on visual features from the large-scale ImageNet dataset. Using this network as a pretrained model on our FCS and LCF15 fish datasets, an informative visual representation of fish objects and their motion can be extracted. After applying the pretrained ResNet-152 network on the input which is a concatenation of the raw greyscale video frames and the motion candidates generated by GMM and optical flow, we get the output features. This three-input combination is alternative to the standard three-channel RGB image. The output features extracted by applying ResNet-152 are fed into RPN to generate candidate regions where fish might be present. This is achieved by sliding a small window of size 3 × 3 3×3 on each of the feature maps to produce k proposals, called anchor boxes, of different aspect ratio and scale. We use three different scales (128×128, 256×256, 512×512) each with 3 different aspect rations (2:1, 1:1, and 1:2) to make k=9 proposals. The aim of using different proposals is to capture fish of different sizes that may appear in an image. These proposals are then classified with a binary classification layer of the RPN to detect the ROI. Another sibling layer of RPN outputs coordinate encodings for each classified proposal. This operation is depicted in Figure 4. The ROIs proposed by the RPN are pooled using an ROI pooling layer and passed onto the final classification head which refines and classifies the proposed ROI into the actual number of classes present at hand, namely fish and non-fish. The complete network is trained in an end-to-end fashion using the features generated by ResNet-152 model as the input and the corresponding ground truths provided by the dataset. While training, we employ an error backpropagation algorithm (Hinton et al., 2006). Figure 4 Open in new tabDownload slide Illustration of the functionality of a RPN to detect and localize fish. The proposal with the best fit to the fish instance is selected out of k choices. As mentioned earlier, the parameters of the GMM-based motion detection algorithm are chosen such that it detects even a very small motion by either fish or non-fish objects producing high false alarm or recall rates. The R-CNN architecture, which is a combination of the ResNet-152 based feature extraction and RPN followed by a final classification layer for localizing moving objects, refines the output of the GMM and optical flow motion candidates. Therefore, the information of motion coming from GMM and optical flow is fed into R-CNN to finally detect and localize objects. Apart from motion candidates generated by GMM and optical flow, the use of greyscale raw images in combination with motion candidates as input to the ResNet-152 CNN helps in preserving the textural information of fish appearing in the video frame, which increases the capability of the network to induce separability between fish and non-fish objects. The reason of using greyscale image instead of RGB to fine-tune the R-CNN is the observation that colour information in the employed datasets are not distinct enough to enhance the accuracy of detection as the background is also vibrant in colours. Moreover, doing so increases the computational overhead. In this work, we utilized computer systems equipped with Intel Core-i7 processors and Nvidia Titan X graphical processing units (GPUs). The proposed system is trained and tested using TensorFlow deep learning library (https://www.tensorflow.org) with Tf-faster-rcnn version while GMM and optical flow source codes were taken from publicly available authors’ repositories (https://github.com/andrewssobral/bgslibrary). Fish detection system utility Our software system is available for deployment and ready to be used by marine scientists for automatic fish detection in any dataset. As described in Region-based convolutional neural network Section, the deep network, which is the backbone of our algorithm, is pretrained on a large and generic object image repository called ImageNet and acts as a generic feature extractor. However, for using a pretrained network in such a transfer learning approach, the system must be fine-tuned to the actual datasets in hand; therefore, a complete end-to-end re-training on a new dataset is not required. In our case, we utilized FCS and LCF-15 datasets by updating the weights of the top fully connected layers of ResNet-152 of R-CNN, while keeping the lower layers intact. Furthermore, the GMM and optical flow algorithms can be used as is since they only require the available dataset to generate output. The source code for our proposed hybrid system is available for download from the following repository: https://github.com/ahsan856jalal/Fish-Abundance. Scientists can use this code off-the-shelf for fish/object detection in any dataset, video recordings or even still imagery. Results The underwater video background is modelled by GMM using training data from the initial few frames of the video while the remainder of the video is treated as the test dataset. Since each video in our datasets has a different background, we need to keep the first N frames of each video for background modelling. We take N=50 in our experiments as this value was chosen on a trial basis to get optimum GMM performance on our datasets. Smaller values of N produces an inferior performance, while increasing beyond this value does not bring any improvement and increases GMM training time. Optical flow does not require any training data but simply uses adjacent frames to calculate a motion representation. The R-CNN, on the other hand, requires more data to tune the weight parameters for refined motion detection. The raw video images and the motion candidates generated by the GMM and optical flow are fed to the R-CNN for training. One video from each of the seven categories of FCS dataset is set aside for training the GMM and R-CNN. On top of that, GMM also requires the first 50 frames of each video to make a background model and to generate a blob of moving objects in the test frames. The LCF-15 dataset, on the other hand, is already segmented into training and test sets, 20 videos out of a total of 93 are used in training and the remainder is used for testing. Once again, GMM models for all 93 videos are created using the N initial frames. Table 2 lists the performance measure for the fish detection task as an F-measure (Palazzo and Murabito, 2014) for our proposed hybrid system and its independent constituents of GMM, optical flow and standalone R-CNN, which are trained on raw RGB images from videos. Table 2. Performance analysis of individual components of our proposed hybrid framework in comparison to their joint accuracy. Dataset . GMM . Optical flow . R-CNN . Our hybrid system . FCS Blurred 77.80 45.94 85.62 86.76 Complex background 75.94 49.77 52.74 89.54 Crowded 74.41 67.48 53.23 84.27 Dynamic background 64.30 44.62 62.06 90.36 Luminosity change 59.07 58.67 70.17 81.44 Camouflage foreground object 70.03 67.00 66.25 89.97 Hybrid videos 75.50 59.44 64.90 91.50 Average 71.01 56.13 64.99 87.44 LCF-15 76.21 52.73 77.30 80.02 Dataset . GMM . Optical flow . R-CNN . Our hybrid system . FCS Blurred 77.80 45.94 85.62 86.76 Complex background 75.94 49.77 52.74 89.54 Crowded 74.41 67.48 53.23 84.27 Dynamic background 64.30 44.62 62.06 90.36 Luminosity change 59.07 58.67 70.17 81.44 Camouflage foreground object 70.03 67.00 66.25 89.97 Hybrid videos 75.50 59.44 64.90 91.50 Average 71.01 56.13 64.99 87.44 LCF-15 76.21 52.73 77.30 80.02 F-scores (in percentage) for three different methods i.e. GMM (Stauffer and Grimson, 1999), Optical flow (Warren and Strelow, 1985), and R-CNN (Ren et al., 2017) on FCS and datasets for seven categories of video complexity. Highest scores are highlighted in bold Open in new tab Table 2. Performance analysis of individual components of our proposed hybrid framework in comparison to their joint accuracy. Dataset . GMM . Optical flow . R-CNN . Our hybrid system . FCS Blurred 77.80 45.94 85.62 86.76 Complex background 75.94 49.77 52.74 89.54 Crowded 74.41 67.48 53.23 84.27 Dynamic background 64.30 44.62 62.06 90.36 Luminosity change 59.07 58.67 70.17 81.44 Camouflage foreground object 70.03 67.00 66.25 89.97 Hybrid videos 75.50 59.44 64.90 91.50 Average 71.01 56.13 64.99 87.44 LCF-15 76.21 52.73 77.30 80.02 Dataset . GMM . Optical flow . R-CNN . Our hybrid system . FCS Blurred 77.80 45.94 85.62 86.76 Complex background 75.94 49.77 52.74 89.54 Crowded 74.41 67.48 53.23 84.27 Dynamic background 64.30 44.62 62.06 90.36 Luminosity change 59.07 58.67 70.17 81.44 Camouflage foreground object 70.03 67.00 66.25 89.97 Hybrid videos 75.50 59.44 64.90 91.50 Average 71.01 56.13 64.99 87.44 LCF-15 76.21 52.73 77.30 80.02 F-scores (in percentage) for three different methods i.e. GMM (Stauffer and Grimson, 1999), Optical flow (Warren and Strelow, 1985), and R-CNN (Ren et al., 2017) on FCS and datasets for seven categories of video complexity. Highest scores are highlighted in bold Open in new tab The F-score is calculated as, F=2×Recall×PrecisionRecall+Precision where Precision= True PositivesTrue Positives+False Positives and Recall= True PositivesTrue Positives+False Negatives These scores are computed based on overlap between the areas of bounding boxes related to ground truths and detected fish. An average detection accuracy of 87.44% was achieved by the proposed hybrid system for the FCS dataset for all seven categories of environmental variation. In comparison, the GMM alone yielded an average accuracy of 71.01% exceeding the optical flow and standalone R-CNN with significant margins. We also performed similar experimentation on LCF-15 test dataset of 73 underwater videos. There, our proposed hybrid system outperforms all the other systems, yielding an accuracy of 80.02% as compared with 76.21, 52.73, and 77.30% by the GMM, optical flow and stand-alone R-CNN, respectively. The parameters of the GMM were carefully chosen to produce best possible results by altering the variance for model fitting and the number of frames for training the model on each video. A fewer number of training frames per video results in degraded performance. However, increasing the number of training frames beyond 50 did not improve the overall performance significantly. Similarly, for our proposed hybrid system and also for the stand-alone R-CNN trained on the raw RGB images, various state-of-the-art pretrained CNNs were tried that include Inception-V4 (Szegedy et al., 2016) and DenseNet (Huang et al., 2017). All of these networks are pretrained on the ImageNet dataset with the same experimental settings. Moreover, different numbers of convolution layers for the RPN network were also evaluated and the choice of sliding window size of 2 × 2 2×2 , 4× 4, and 5 × 5 was tested, with the performance maximized at 3 × 3 3×3 . The performance started to deteriorate slightly beyond the 3 × 3 3×3 window size probably due to more overlap between intrinsic size of fish covering the frame of videos in our datasets. The results generated by Inception-V4 and ResNet-152 were comparable without any significant difference but the latter utilizes less processing power in training and testing compared with the former. Our implementation of optical flow on the other hand is a non-trainable processing approach for motion detection and therefore, does not have any trainable parameters. It is worth mentioning here that the GMM chosen for our proposed hybrid system differs with the one listed in Table 2 as its parameters were tuned to produce higher recall rates at the cost of decreased precision to cover maximum possible pixel motion in the video by both fish and non-fish objects. The CNN and RPN subsystems then learn to select the relevant motion candidate through refining the results generated by the GMM and optical flow. Figure 5 shows the performance outcome on a sample video for GMM, optical flow, R-CNN, and the proposed hybrid system for both the FCS and the LCF-15 datasets. It is evident that the optical flow algorithm generates more false alarms and is sensitive to even very slight motion, which can be attributed to disturbances in the scene or luminosity changes. On the other hand, the GMM and stand-alone R-CNN, which is only trained on raw RGB images, also exhibits false alarms and/or missed detection. However, they both yield better scores as compared with the optical flow due to effective background modelling and end-to-end supervised training; capabilities which optical flow lacks and are necessary to reduce the irrelevant motion created by non-fish entities. Our proposed hybrid system, on the other hand is successful in achieving the best performance (see Table 2). Figure 5. Open in new tabDownload slide Example of fish detection outcomes by various algorithms. Left to right, ground truth, optical flow, GMM, stand-alone R-CNN, and proposed hybrid system on all seven categories of FCS dataset category (the first seven rows) and one video of LCF-15 dataset (the last row). To validate the effectiveness of our system, in Table 3 we have drawn a comparison with various published benchmark approaches which are frequently used for motion-based object detection in either still or video imagery. The comparison is made on the FCS dataset for which we can directly tabulate published scores by these techniques with the same experimental settings as ours. It is evident that our proposed hybrid system outperforms all others in most environmental conditions and the overall average F-scores. In another set of experimentation not reported here, we changed the train-test split in the FCS and LCF-15 datasets to calculate the detection scores but observed no significant change. This demonstrates a good generalization capability of our system. Table 3. F-scores (in percentage) for different methods on FCS datasets for fish detection, as given in Spampinato et al. (2014). Video class . KDE . ML-BKG . EIGEN . VIBE . TKDE . Our hybrid system . Blurred 92.56 70.26 81.71 85.13 93.25 86.76 Complex background 87.53 83.67 74.78 74.17 81.79 89.54 Crowded 82.46 79.81 73.87 84.64 84.19 84.27 Dynamic background 59.13 77.51 71.48 67.01 75.59 90.36 Luminosity change 72.06 82.66 70.41 70.37 72.95 81.44 Camouflage foreground object 54.14 73.51 70.20 76.30 82.23 89.97 Hybrid videos 85.69 72.20 80.69 79.75 82.63 91.50 Average 76.22 77.08 74.73 76.76 81.80 87.44 Video class . KDE . ML-BKG . EIGEN . VIBE . TKDE . Our hybrid system . Blurred 92.56 70.26 81.71 85.13 93.25 86.76 Complex background 87.53 83.67 74.78 74.17 81.79 89.54 Crowded 82.46 79.81 73.87 84.64 84.19 84.27 Dynamic background 59.13 77.51 71.48 67.01 75.59 90.36 Luminosity change 72.06 82.66 70.41 70.37 72.95 81.44 Camouflage foreground object 54.14 73.51 70.20 76.30 82.23 89.97 Hybrid videos 85.69 72.20 80.69 79.75 82.63 91.50 Average 76.22 77.08 74.73 76.76 81.80 87.44 The results of our proposed system are copied from Table 2 for easy comparison in this table. Highest scores are highlighted in bold. Open in new tab Table 3. F-scores (in percentage) for different methods on FCS datasets for fish detection, as given in Spampinato et al. (2014). Video class . KDE . ML-BKG . EIGEN . VIBE . TKDE . Our hybrid system . Blurred 92.56 70.26 81.71 85.13 93.25 86.76 Complex background 87.53 83.67 74.78 74.17 81.79 89.54 Crowded 82.46 79.81 73.87 84.64 84.19 84.27 Dynamic background 59.13 77.51 71.48 67.01 75.59 90.36 Luminosity change 72.06 82.66 70.41 70.37 72.95 81.44 Camouflage foreground object 54.14 73.51 70.20 76.30 82.23 89.97 Hybrid videos 85.69 72.20 80.69 79.75 82.63 91.50 Average 76.22 77.08 74.73 76.76 81.80 87.44 Video class . KDE . ML-BKG . EIGEN . VIBE . TKDE . Our hybrid system . Blurred 92.56 70.26 81.71 85.13 93.25 86.76 Complex background 87.53 83.67 74.78 74.17 81.79 89.54 Crowded 82.46 79.81 73.87 84.64 84.19 84.27 Dynamic background 59.13 77.51 71.48 67.01 75.59 90.36 Luminosity change 72.06 82.66 70.41 70.37 72.95 81.44 Camouflage foreground object 54.14 73.51 70.20 76.30 82.23 89.97 Hybrid videos 85.69 72.20 80.69 79.75 82.63 91.50 Average 76.22 77.08 74.73 76.76 81.80 87.44 The results of our proposed system are copied from Table 2 for easy comparison in this table. Highest scores are highlighted in bold. Open in new tab Discussion In this study, we have proposed a R-CNN to detect fish using enhanced features sensitive to natural fish motion in underwater videos in addition to features also representing distinguishable shape and textural information specific to fish in a supervised training hierarchy. The motivation behind using such a deep neural network is to model complex and highly nonlinear attributes in underwater imagery of fish. These attributes are not modelled effectively by conventional machine learning algorithms and image processing techniques (Hinton and Salakhutdinov, 2006; Larochelle et al., 2009). This hybrid approach has resulted in a detection accuracy at reasonable level for use of this technique in fish detection from recorded videos. The most important gain of this research is high detection accuracy of freely swimming fish. With our proposed hybrid system that incorporates motion sensitive features, taken as input to the R-CNN, we are able to achieve 87.44% detection accuracy on the FCS dataset. This performance exceeds the best reported results on this dataset by a significant margin. The second best average accuracy of 81.80% for all seven categories of variability has been produced using KDE to model background and foreground objects by capturing texture information in very low contrast regions of the video frames (Spampinato et al., 2014). An interesting observation can be drawn from Table 3 for video classes Dynamic background, Camouflage foreground object and Hybrid videos that the performance gap between our proposed hybrid system and rest of the techniques is significantly wide. Dynamic background videos exhibit disturbance in water surface and movement of aquatic plants which causes confusion with motion of fish. Therefore, KDE, ML-BKG, and TKDE algorithms, which are based on estimating foreground data distribution by modelling background data, fails in separating motion of fish and non-fish objects. EIGEN and VIBE algorithms also produced poor performance due to similar reasons. Here, our proposed hybrid system utilizes the fish-dependent features captured through the R-CNN component using greyscale images in accurate detection of fish. On the other hand, fish in Camouflage foreground object videos are extremely hard to segregate from the background. Therefore, all the algorithms once again fail to yield better results due to inability in creating difference between foreground and background models. Here, our approach makes use of the motion information from GMM and optical flow to maximize its fish detection potential as shape, texture and colour of fish in this case resemble the background and are difficult to detect by the R-CNN component. Similarly, Hybrid videos combine all the challenges of other six classes and our proposed hybrid system is more effective than all other approaches. To further endorse the effectiveness of our approach, we employed a larger dataset by including LCF-15 with 93 videos. Our solution acquired an average accuracy of 80.02%. Table 2 lists the comparative performance of our proposed hybrid system with three other techniques, namely GMM, optical flow and R-CNN, which are the components of our overall system. The GMM outperforms optical flow and standalone R-CNN, trained on raw images, with a significant margin, for the FCS dataset. On the LCF-15 dataset, the GMM produces better results than optical flow and is comparable with the R-CNN. This signifies effective learning of the background model by the GMM on every new video sequence. The model covers all background variations exhibited by non-fish objects for a static underwater camera configuration, which assists in detecting even subtle movements through non-uniform change in pixel intensities that does not match with the distribution of background pixels. We observe that the training of the GMM background model balances the rate of false alarm and misdetection, which produces a better F-score. The GMM and its variants are considered to give excellent performance in general for motion-based object detection tasks (Yong, 2013; Spampinato et al., 2014). Optical flow, on the other hand, lagged behind all other methods in terms of performance on both datasets. The reason behind this behaviour can be attributed to the non-trainable structure of this algorithm, as the system cannot adapt to the dynamic environment in the videos. There is no learning involved to discriminate background and foreground modelling, like that in the GMM and neural networks. Optical flow involves a direct comparison between adjacent frames of video and any slight disturbance in the pixel intensity, either due to fish or non-fish objects generating luminosity variation, translates into a valid motion. This gives rise to numerous false alarms, which results in a very high recall but consequently a low precision that ends up in producing low F-score. Since the datasets we have chosen involve high environmental variation, especially in the FCS dataset, optical flow fails to perform well as opposed to the other algorithms. On comparative grounds, both GMM and optical flow lags our proposed hybrid system for the FCS and LCF-15 datasets. Another explanation for relatively worse performance of these approaches, as compared with the proposed system shown in Table 2, is an important observation that can be made by watching the videos in the datasets. The fish in each frame may not necessary show motion and sometimes remain dormant for multiple frames, even though for most of the time they are swimming, making the scenes dynamic. The GMM sometimes confuses the fish with the stationary profile as background, especially when the appearance of fish matches the background. Therefore, lack of motion information in video frames results in failure to detect fish by the GMM and optical flow. The R-CNN on the other hand is a tailored neural network used for object localization in the images (Ren et al., 2017) and learns to capture fish-dependent information from stationary images. Underwater fish detection in unconstrained environments is a challenging task as the main aim lies in segregating fish and ignoring non-fish entities in the entire video frame. Conventional machine learning and image processing algorithms are generally designed to detect the objects of interest in the datasets where they exhibit their distinct presence in the imagery, and hence are easier to segment out (Russakovsky et al., 2015). In contrast, a high degree of confusion in separating fish with vibrant, diverse and variable non-fish objects in underwater videos results in a performance compromise for a standalone R-CNN with accuracy of 64.99 and 77.30% on the FCS and LCF-15 datasets, respectively. As mentioned earlier, many videos, especially in the FCS dataset, lacks textural and shape information of fish, a necessary ingredient to yield better performance by systems like standalone R-CNN. This problem is effectively solved by our proposed hybrid system using and learning the information from motion-sensitive and textural features. Figure 6 shows some results from the FCS and LCF-15 datasets where all algorithms including our proposed system failed to detect fish. These are the extreme cases of blurriness, camouflage, water murkiness, and unrecognizable orientation, texture, and shape of fish which either results in generating false alarms or miss detections. In these situations, it is extremely difficult to capture both motion-based and shape/texture-based features. Figure 6. Open in new tabDownload slide Examples of false detection of fish by all algorithms including our proposed hybrid system. Here, bounding boxes signify either miss detections of fish or false alarms. The black and white images are corresponding ground truths. In the future, we aim to employ a unified deep architecture capable of processing the video sequences in real-time through rigorous optimization of our algorithm and better mathematical modelling. Such a setup will be applicable for fish detection as well as their species classification at the same time and, therefore, will be more suitable for effective fish fauna sampling. Furthermore, the accuracy of the system can be improved by tracking the paths of moving fish and having prior information of their movement in several frames. This step can improve the accuracy of detection in the video frames where the proposed approach fails to recognize fish due to extreme blurriness and the camouflage of the background. We plan to incorporate this processing step using recurrent neural networks (Gordon et al., 2018) with temporal processing capability in videos. Conclusions In this article, we have presented an automatic method that employed deep R-CNN networks to detect and localize fish instances in unconstrained underwater videos that exhibit various degrees of scene complexity. The major contribution of this work is that it utilizes a hybrid approach involving GMM and optical flow outputs to combine motion sensitive input features with raw video frames carrying textural and shape information. This mixed data is used as input to a deep R-CNN to fine-tune the categorization of fish in the presence of non-fish entities in the video frame. This assisted in achieving state-of-the-art results for the fish detection task as confirmed by the comparative study. The proposed hybrid system requires relatively more computational resources as compared with the conventional computer vision and machine learning techniques, but comes with the benefit of higher accuracy. However, with an advent of fast microprocessors and GPUs, complex mathematical operation involved in deep neural networks like CNN can be performed quickly, even making them suitable for tasks requiring near real-time processing. Therefore, combining the hybrid fish detection with other fish-related tasks like fish classification even using deep learning (Salman et al., 2016) and tracking can be made possible in the pursuit of realizing fully automated systems for deployment in real world applications of fisheries. We believe that this research will help scientists related to fisheries in adopting automatic approaches for detection, classification and tracking of fish fauna in non-destructive sampling. Moreover, in the future, we aim to employ a unified deep architecture capable of processing the video sequences in real-time through rigorous optimization of our algorithm and better mathematical modelling. Such a setup will be applicable for fish detection as well as their species classification at the same time and therefore, will be more suitable for effective fish fauna sampling. Acknowledgements The authors acknowledge support from the Australian Research Council Grant LP110201008, and German Academic Exchange Service (DAAD) Project Grant 57243488 “FIBEVID”. The authors also acknowledge Nvidia Corporation, USA for their donation of graphics processing units (GPUs) under their GPU Grant Programme. Nvidia GPUs were used to carry out experiments in the work carried out in this article. References Bengio Y. 2009 . Learning deep architectures for AI . Foundations and Trends in Machine Learning , 2 : 1 – 127 . Google Scholar Crossref Search ADS WorldCat Beauchemin S. S. , Barron J. L. 1995 . The computation of optical flow . ACM Computing Surveys , 27 : 433 – 466 . Google Scholar Crossref Search ADS WorldCat Brox T. , Bruhn A., Papenberg N., Weickert J. 2004 . High accuracy optical flow estimation based on theory for warping. In Computer Vision-ECCV 2004. Lecture Notes in Computer Science , 3024 . Ed. by Pajdla T., Matas J.. Springer , Berlin, Heidelberg . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Deng J. , Dong W., Socher R., Li L-J., Li K., Fei-Fei L. 2009 . ImageNet: a large-scale hierarchical image database . IEEE CVPR-2009 , Miami, FL, USA , 248 – 255 pp. Google Scholar OpenURL Placeholder Text WorldCat Fisher R. , Chen-Burger Y-H., Giordano D., Hardman L., Lin F-P. (Eds) 2016 . Fish4Knowledge: Collecting and Analyzing Massive Coral Reef Fish Video Data. Intelligent Systems Reference Library, 104 , Springer International Publishing . DOI: 10.1007/978-3-319-30208-9. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Froese R. 2006 . Cube law, condition factor and weight length relationships: history, meta-analysis and recommendations . Journal of Applied Ichthyology , 22 : 241 – 253 . Google Scholar Crossref Search ADS WorldCat Gordon D. , Farhadi A., Fox D. 2018 . Re3: real-time recurrent regression networks for visual tracking of generic objects . IEEE Robotics and Automation Letters , 3 : 788 – 795 . Google Scholar Crossref Search ADS WorldCat Harvey E. S. , Shortis M. R. 1995 . A system for stereo-video measurement of sub-tidal organisms . Marine Technology Society Journal , 29 : 10 – 22 . Google Scholar OpenURL Placeholder Text WorldCat He K. , Zhang X., Ren S., Sun J. 2016 . Deep residual learning for image recognition . IEEE CVPR-2016 , Las Vegas, NV, USA , 770 – 778 pp. Google Scholar OpenURL Placeholder Text WorldCat Hinton G. , Salakhutdinov R. 2006 . Reducing the dimensionality of data with neural networks . Science , 313 : 504 – 507 . Google Scholar Crossref Search ADS PubMed WorldCat Hinton G. E. , Osindero S., Teh Y-W. 2006 . A fast learning algorithm for deep belief nets . Neural Computation , 18 : 1527 – 1554 . Google Scholar Crossref Search ADS PubMed WorldCat Hsiao Y. , Chen C., Lin S., Lin F. 2014 . Real-world underwater fish recognition and identification using sparse representation . Ecological Informatics , 23 : 13 – 21 . Google Scholar Crossref Search ADS WorldCat Huang G. , Liu Z., Weinberger K. Q. 2017 . Densely connected convolutional networks . IEEE CVPR-2017 , Honolulu, HI, USA , 2261 – 2269 pp. Google Scholar OpenURL Placeholder Text WorldCat Jennings S. , Kaiser M. J. 1998 . The effects of fishing on marine ecosystems . Advances in Marine Biology , 34 : 201 – 352 . Google Scholar Crossref Search ADS WorldCat Larochelle H. , Bengio Y., Louradour J., Lamblin P. 2009 . Exploring strategies for training deep neural networks . Journal of Machine Learning Research , 10 : 1 – 40 . Google Scholar OpenURL Placeholder Text WorldCat LeCun Y. , Boser B., Denker J. S., Henderson D., Howard R. E., Hubbard W., Jackel L. D. 1989 . Backpropagation applied to handwritten zip code recognition . Neural Computing , 1 : 541 – 551 . Google Scholar Crossref Search ADS WorldCat LeCun Y. , Huang F., Bottou L. 2004 . Learning methods for generic object recognition with invariance to pose and lighting . IEEE CVPR-2004 , Washington, DC, USA , 97 – 104 pp. Google Scholar OpenURL Placeholder Text WorldCat LeCun Y. , Bengio Y., Hinton G. 2015 . Deep Learning . Nature , 521 : 436 – 444 . Google Scholar Crossref Search ADS PubMed WorldCat Lin T. Y. , RoyChowdhury A., Maji S. 2015 . Bilinear CNN models for fine-grained visual recognition . IEEE ICCV-2015 , Santiago, Chile , 1449 – 1457 pp. Google Scholar OpenURL Placeholder Text WorldCat McLaren B. W. , Langlois T. J., Harvey E. S., Shortland-Jones H., Stevens R. 2015 . A small no-take marine sanctuary provides consistent protection for small-bodied by-catch species, but not for large-bodied, high-risk species . Journal of Experimental Marine Biology and Ecology , 471 : 153 – 163 . Google Scholar Crossref Search ADS WorldCat Mika S. , Ratsch G., Weston J., Scholkopf B., Mullers K. R. 1999 . Fisher discriminant analysis with kernels . IEEE International Workshop on Neural Networks for Signal Processing , Madison, WI, USA , 41 – 48 pp. Google Scholar OpenURL Placeholder Text WorldCat Moniruzzaman M. , Islam S. M. S., Bennamoun M., Lavery P. 2017 . Deep learning on underwater marine object detection: a survey . In Advanced Concepts for Intelligent Vision Systems. ACIVS 2017. Ed. by J. Blanc-Talon, R. Penne, W. Philips, D. Popescu, and P. Scheunders. Lecture Notes in Computer Science, vol 10617 , Springer, Cham. Google Scholar OpenURL Placeholder Text WorldCat Palazzo S. , Murabito F. 2014 . Fish species identification in real-life underwater images. In 3rd ACM International Workshop on Multimedia Analysis for Ecological Data , Orlando, Florida , pp. 13 – 18 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Ren S. , He K., Girshick R., Sun J. 2017 . Faster R-CNN: towards real-time object detection with region proposal networks . IEEE Transactions on Pattern Analysis and Machine Intelligence , 39 : 1137 – 1149 . Google Scholar Crossref Search ADS PubMed WorldCat Reynolds D. A. , Rose R. C. 1995 . Robust text-independent speaker identification using Gaussian mixture speaker models . IEEE Transactions on Acoustics, Speech and Signal Processing , 3 : 72 – 83 . Google Scholar OpenURL Placeholder Text WorldCat Russakovsky O. , Deng J., Su H., Krause J., Satheesh S., Ma S., Huang Z. et al. 2015 . ImageNet large scale visual recognition challenge . International Journal of Computer Vision , 115 : 211 – 252 . Google Scholar Crossref Search ADS WorldCat Salman A , Jalal A., Shafait F., Mian A., Shortis M., Seager J., Harvey E. 2016 . Fish species classification in unconstrained underwater environments based on deep learning . Limnology and Oceanography: Methods , 14 : 570 – 585 . Google Scholar Crossref Search ADS WorldCat Siddiqui S. A. , Salman A., Malik M. I., Shafait F., Mian A., Shortis M. R., Harvey E. S. 2017 . Automatic fish species classification in underwater videos: exploring pre-trained deep neural network models to compensate for limited labelled data . ICES Journal of Marine Sciences , 75 : 374 – 389 . Google Scholar Crossref Search ADS WorldCat Sheikh Y. , Shah M. 2005 . Bayesian modelling of dynamic scenes for object detection . IEEE Transaction on Pattern Analysis and Machine Intelligence , 27 : 1778 – 1792 . Google Scholar Crossref Search ADS WorldCat Shortis M. , Harvey E. S., Abdo D. 2009 . A review of underwater stereo-image measurement for marine biology. In Oceanography and Marine Biology: An Annual Review . Ed. by Gibson R. N., Atkinson R.J.A., Gordon J.D.M.. CRC Press , USA . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Simonyan K. , Zisserman A. 2014 . Very deep convolutional networks for large-scale image recognition . arXiv : 1409.1556. Google Scholar OpenURL Placeholder Text WorldCat Spampinato C. , Chen-Burger Y., Nadarajan G., Fisher R. B. 2008 . Detecting, tracking and counting fish in low quality unconstrained underwater videos . International Conference on Computer Vision Theory and Applications , Funchal, Madeira, Portugal , 2 : 514 – 519 . Google Scholar OpenURL Placeholder Text WorldCat Spampinato C. , Palazzo S., Kavasidis I. 2014 . A texton-based kernel density estimation approach for background modeling under extreme conditions . International Journal of Computer Vision and Image Understanding , 122 : 74 – 83 . Google Scholar Crossref Search ADS WorldCat Storbeck F. , Daan B. 2001 . Fish species recognition using computer vision and a neural network . Fisheries Research , 51 : 11 – 15 . Google Scholar Crossref Search ADS WorldCat Strachan N. J. C. , Kell L. 1995 . A potential method for the differentiation between haddock fish stocks by computer vision using canonical discriminant analysis . ICES Journal of Marine Science , 52 : 145 – 149 . Google Scholar Crossref Search ADS WorldCat Stauffer C. , Grimson W. E. L. 1999 . Adaptive background mixture models for real-time tracking . IEEE CVPR-1999, Fort Collins, CO, USA , 2 : 246 – 252 . Google Scholar OpenURL Placeholder Text WorldCat Sung M. , Yu S., Girdhar Y. (2017) . Vision based real-time fish detection using convolution neural network . IEEE OCEAN-2017 , Aberdeen, UK , 1 – 6 pp. Google Scholar OpenURL Placeholder Text WorldCat Szegedy C. , Ioffe S., Vanhoucke V. 2016 . Inception-v4, inception-resnet and the impact of residual connections on learning . AAAI Conference on Artificial Intelligence , Phoenix, AZ, USA, 4278 – 4284 pp. Google Scholar OpenURL Placeholder Text WorldCat Tanzer J. , Phua C., Lawrence A., Gonzales A., Roxburgh T., Gamblin P. (Eds) 2015 . Living Blue Planet Report. Species, Habitats and Human Well-Being . WWF , Gland . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Turk M. , Pentland A. 1991 . Eigenfaces for recognition . Journal of Cognitive Neuroscience , 3 : 71 – 86 . Google Scholar Crossref Search ADS PubMed WorldCat Warren D. H. , Strelow E. R. 1985 . Electronic Spatial Sensing for the Blind: Controbutions from Perception . Springer . ISBN 90-247-2689-1. Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Yao J. , Odobez J. M. 2007 . Multi-layer background subtraction based on color and texture . IEEE CVPR-2007, Minneapolis, MN, USA , 1 – 8 pp. Google Scholar OpenURL Placeholder Text WorldCat Yong X. 2013 . Improved Gaussian mixture model in video motion detection . Journal of Multimedia , 8 : 527 – 533 . Google Scholar OpenURL Placeholder Text WorldCat Zivkovic Z. , Heijden F. 2006 . Efficient adaptive density estimation per image pixel for the task of background subtraction . Pattern Recognition Letters , 27 : 773 – 780 . Google Scholar Crossref Search ADS WorldCat © International Council for the Exploration of the Sea 2019. All rights reserved. For permissions, please email: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) © International Council for the Exploration of the Sea 2019. All rights reserved. For permissions, please email: [email protected]
Highlighting growth regulation processes in fish populations by a simplex simulation approach: application to Merluccius hubbsi stocks in the Southwestern AtlanticSemmar,, Nabil;Vaz-dos-Santos, André, M
doi: 10.1093/icesjms/fsz240pmid: N/A
Abstract A new simplex-based simulation approach (Spx) was developed to highlight multidirectional and multi-scale relationships between morphometric variables helping to functionally differentiate biological (fish) groups for better stocks definition and monitoring. Application concerned Merluccius hubbsi sampled in 1968–1972 and 2004 in six Southwestern Atlantic areas. Simulation results highlighted negative trends opposing front to back compartments indicating competition for body biomass distribution. However, top and bottom parts within these compartments were positively correlated indicating cooperative processes in favour of target local growth regulation. Positive and negative trends of growth regulation were also highlighted at lower body scale, notably between smaller components constituting body front compartment. On a geographical scale, average regulation levels of same morphometric variables showed monotonic or alternated variations between successive fish groups. This highlighted target and modulated growth regulations governing biomass distribution in different body parts by geographical-dependent ways. Under dynamical aspect (1968–1972 vs. 2004), growth regulation of mouth tended to increase with time leading to conclude on morphometric responses of M. hubbsi to overfishing pressure. Spx results were confirmed by several traditional approaches which showed less integrative aspect. Introduction Growth represents complex process associated with variable distributions of biomass at several biological scales extending from individual body to intra- and inter-population balance states. Growth states and mechanisms serve as key criteria to define groups with homogeneous vital rates (stock units in fishery management; Kerr et al., 2017). The questions of identity, assessment, and monitoring of unit stocks represent classical aims in fishery management. Complexity of growth system is linked to multifactorial, multidirectional, and multi-scale regulation trends between body components (morphometric variables) leading to specific aspects of different biological groups. Multiple aspects of growth in biological populations are classically approached by means of several complementary quantitative methods (Cadrin, 2000; Cadrin et al., 2014; Secor, 2014): Principal component analysis (PCA) is unsupervised method providing correlation charts between variables helping to highlight different trends associated with different groups (Semmar, 2011). Multiple correspondence analysis (MCA) highlights how different groups are associated with different variation ranges and ways of discrete variables (Escofier and Pagès, 1991). Discriminant analysis (DA) is supervised method used to identify/predict groups from quantitative profiles of measured (morphometric) variables (Seber, 1984). Linear and non-linear modelling methods provide relational aspects of growth governed by a limited number of significant variables (Seber, 1984). Body sizes and weights are generally used to model growth processes by traditional methods including Huxley and von-Bertalanffy models (Huxley, 1924, 1932; Weatherley and Gill, 1987; Weatherley, 1990; Saborido-Rey and Kjesbu, 2005). Alternatively to several decomposition approaches, simulation methods provide integrative structured information on a studied population by exploring potentially its overall variability. Simulation methods helping to treat with multidimensional, multi-scale, and overlapping/mixing aspects of biological groups were encouraged in the recent years under the framework of management evaluation strategy (MES; Kerr et al. 2017). Integrative aspect can be satisfied by Markov Chain Monte Carlo sampling methods (Manly, 2007, Panikian et al., 2015). The current work presents a new simulation approach based on population stratification, simplex mixture design, and bootstrap sampling technique. Simplex rule is appropriate for mass distribution analysis between system constituents and through organization ways obeying to mass conservation law (Semmar and Roux, 2014). Stratification and bootstrap provided robust sampling ways for heteroscedasticity treatment and variability integration, respectively. The new simulation simplex method (Spx) treats mixing aspects of biological groups constituting heterogeneous population and characterized by variability of morphometric profiles. In silico differential mixing and averaging of morphometric profiles helped to highlight growth regulation trends manifesting within and between groups. Mixing-based methodology of Spx is appropriate to treat the natural connections and overlapping aspect of biological groups. Spx allowed fish stock identification, monitoring, and management from morphometric variability by responding to key questions: (i) “how length-expressed biomass was relatively distributed through different body parts conditionally to different fish groups considered a priori?” and (ii) “how lengths of different body parts can reflect specific growth ways of ecological populations in relation to intrinsic and/or external governing factors?” Spx was applied on morphometric variables of the Argentine hake Merluccius hubbsi to analyse growth regulations of different body parts in fish groups associated with several geographical areas. Despite the overlapping aspects in natural populations, Spx was able to extract differentiation mechanisms specific to different groups helping for functional stock management. Merluccius hubbsi Marini 1933 is one of the most important fishing resources in the Southwestern Atlantic Ocean, east coast of South America (Lloris et al., 2005). It is a demersal species continuously distributed from Brazil (21°S) to Argentina (55°S) in the continental shelf and slope, associated with water temperatures up to 23°C. In Brazilian waters (21°–34°S), two stocks (I and II) were defined (Southeastern, SE: 21°–29°S; Southern, S: 29°–34°S) and differentiated by several parameters including spawning’s areas and periods, first maturity, growth rate, otolith ring formation, otolith morphology, and young-of-the-year development (Vaz-dos-Santos et al., 2009, 2017; Costa et al., 2018). Nevertheless, there is a southward gradient favouring dynamical exchanges between fish populations leading to mixing aspects of groups with unknown relative levels (Kerr et al., 2017). This leads to misleading risks that need to be overcome to avoid bias in stock identification, assessment, and management. Spx provided appropriate simulation approach to highlight how different groups can be functionally distinguished within their heterogeneous and naturally interconnected whole biosystem. Spx was applied on two sampling periods of M. hubbsi (A: 1968–1972 vs. B: 2004). Six geographical sites were concerned (all in A vs. four in B) and were parts of the two conventional fish stocks (I and II). Traditional statistical analyses were performed on the two stocks (I and II) using log-transformed length data. Results showed limited differentiation between stocks. More integrative information was provided by Spx which used morphometric variables standardized by their sum (relative levels) to simulate growth regulation mechanisms governing variability and differentiation between and within population groups (stock units). This integrative approach is the first of its kind in fishing management field and it responds well to the need of MES development (Kerr et al., 2017). Material and methods Study area The studied area concerned the Southwestern Atlantic Ocean between 21° and 34°S: Brazilian coast and its Economic Exclusive Zone (Figure 1). It included two defined stocks (I and II) associated with the Southeastern (21°–29°S) and Southern (29°–34°S) areas, respectively. Figure 1. Open in new tabDownload slide Merluccius hubbsi: sample sites in the two study periods 1968–1972 (A) and 2004 (B). Strata Aj and Bj were considered in Spx for simulations of growth regulation ways in different fish groups of periods A, B. I, II, conventionally defined stocks. Figure 1. Open in new tabDownload slide Merluccius hubbsi: sample sites in the two study periods 1968–1972 (A) and 2004 (B). Strata Aj and Bj were considered in Spx for simulations of growth regulation ways in different fish groups of periods A, B. I, II, conventionally defined stocks. In the Southeast region, the cold-water mass “South Atlantic Central Water” (6°–20°C, 34–36 PSU) occupies the upper slope during all the year; during the austral spring to summer (October to March) it also moves to the continental shelf, remaining sub superficial (Piola et al., 2018). This is the main water mass in which M. hubbsi lives, although the species is also found in the Shelf Water, a mixture of water masses close to the coast. The South region presents the same pattern, but during the austral autumn to winter (April to September), it experiences intrusions of the Sub-Antarctic Water (12°–15°C, 33 PSU; Piola et al., 2018). Two important upwelling areas occur at 23°S (Cape Frio) and 28°40′S (Cape of Santa Marta Grande). In this second area, M. hubbsi does not spawn and juveniles do not occur (Vaz-dos-Santos and Schwingel, 2015). Sampling and data acquisition Merluccius hubbsi was sampled in two periods A (1968–1972) and B (2004; Figure 1) by reference to the scientific collection of Zoology Museum of São Paulo University (MZUSP) (A) and commercial landings (B). The specimens A (231 fishes) were obtained in bottom-trawl surveys performed in six areas (A1–A6) between 21° and 36°S along the continental shelf (up to 200 m isobaths; Supplementary material S1). The specimens B (516 fishes) were attained during January to February and July to August of 2004 from bottom-trawlers operating along the Brazilian coast (21°–34°S) in four areas (B1–B3, B5) with depths varying between 40 and 590 m (Supplementary material S1). The specimens were cleaned, fixed in formalin (10%), and stored in alcohol (70%) for ulterior analysis. Morphometric data of M. hubbsi concerned 18 body measurements which were attained with a digital calliper (0.05 mm precision; Figure 2;Supplementary material S2;Inada, 1981; Lloris et al., 2005). The specimens B were measured 2–3 d after sampling, avoiding shrinkage bias. Figure 2. Open in new tabDownload slide Merluccius hubbsi: body measurements for morphometric analyses. Figure 2. Open in new tabDownload slide Merluccius hubbsi: body measurements for morphometric analyses. Traditional data analysis Traditional computational analysis of growth processes based on morphometric data involved five steps detailed in Supplementary material S3 and summarized below (Cadrin, 2000): Comparison of length frequency distributions between the two temporal populations A, B using three tests: Kolmogorov–Smirnov, Scheirer-Ray-Hare (non-parametric two-way analysis; an extension of Kruskal–Wallis test), Permutational Multivariate Analysis of Variance (PERMANOVA) (non-parametric permutation test to improve Multivariate Analysis of Variance (MANOVA); Supplementary material S3 and S4; Sokal and Rohlf, 1995; Anderson et al., 2008). Predictive modelling of the variation of 17 morphometric variables Y in relation to the total lengths (X) by applying 17 Huxley models (Y = aXb) (Supplementary material S5; Huxley, 1924, 1932). Correction of morphometric measures by using calculated parameter b and standard length L0 (fixed at 210 mm): Y ′ = Y (L0/Lt)b (Lombarte and Lleonart, 1993). Effect analysis of both fish stock and total length on the 17 morphometric variables by analysis of covariance (ANCOVA) (Sokal and Rholf, 1995). Application of DA for fish stocks (I and II) recognition from multivariate morphometric patterns (Gotelli and Ellison, 2004). Simplex approach Preliminary data processing Spx is initially based on stratification of whole population into several groups characterized by different relative levels of system constitutive variables. Principle of Spx refers to mass conservation law where a whole resource can be distributed by multiple competing ways under unit sum constraint (sum of all distributed mass parts = 1, 100%) (Cornell, 2002). This conservative principle can be applied to measured lengths of different body parts, illustrated here by M. hubbsi. Initially, size effect was removed by relativizing the different measured lengths by reference to their sum that was associated to a whole unit representing the perimeter of the covered body space. Figure 3a and b represents a generic fish with four body measurements combining front/back and top/bottom with illustrative values. Relativization of different body parts by reference to their perimeter (sum) leads to a morphometric profile with different segments varying relatively the ones at the expense of the others (Figure 3b and c). In the current study, two systems of morphometric variables standardized by their sum were considered: (i) Pdd2, Lpa, Lbsd, and Lba were considered to cover overall body perimeter of M. hubbsi. (ii) At lower morphometric scale, body front compartment was characterized by five constitutive variables: Lbfd, Lpso, Ed, Lpo, and Lppv (Figure 2). Figure 3. Open in new tabDownload slide Different methodological steps of simulation simplex approach (illustrated here by q=3 groups). Figure 3. Open in new tabDownload slide Different methodological steps of simulation simplex approach (illustrated here by q=3 groups). Spx initially required population stratification into q groups representing q modalities or levels of a growth-influencing factor (Figure 3c). The effects of such a factor on growth are highlighted by smoothing group-dependent relationships between morphometric variables (Figure 3i). In this study, fish population was stratified by considering geographical sites and periods. This helped to analyse flexibility and stability of fish developments under spatial–temporal aspect. Initially, fish population was stratified into q = 6 and 4 groups sampled in periods A (A1–A6) and B (B1–B3, B5), respectively (Figure 1; Supplementary material S1). Two Spx were applied on A and B datasets, respectively. Simulation methodological steps Growth regulation processes governing length variation of body parts within and between q fish groups were simulated by Spx using a complete set of N combinations between groups (Figure 3e). The N combinations are given by a mixture design (Scheffé’s matrix) containing N rows and q columns (Scheff-, 1958, 1963; Cornell, 2002, 2011). Each row s is associated with a specific mixture combining the q groups j (j = 1 to q) with q weight values wj (contributions) under the constraint of constant sum w: ∑j=1qwj=w(1) With wj = 0 to w. For instance, for q = 3, (w1 + w2 + w3) = (1 + 2 + 7) means a mixture made by 10%, 20%, and 70% of groups 1, 2, and 3, respectively (Figure 3d). The total number N of mixtures is initially given by a combinatorial formula depending on q and w: N=Cq−1w+q−1=(w+q−1)!(q−1)!w!.(2) With w = 10 contributing individuals per mixture of q groups, the formula gives N = 286 and 3003 mixtures with q = 4 and 6, respectively (associated with 4 and 6 fish groups in 2004 and 1968–1972, respectively). Each mixture was applied by randomly sampling w = 10 individuals (fishes) from the q groups j by respecting the q weights wj given by Scheffé’s matrix. At the output of each mixture, the w contributing individuals representing the q weighted groups were averaged to calculate a barycentric profile containing average relative levels of morphometric variables. In all, a response matrix of N barycentric profiles was calculated (Figure 3f). At this step, each barycentric profile was calculated from only w randomly sampled individuals whereas each population group contained much more individuals which did not contribute to mixtures leading to underestimation of variability within and between groups. This sampling deficit was solved by iterating the mixture design and its response matrix by bootstrap technique (Figure 3g; Manly, 2007): in all, K = 30 iterations were applied to the random sampling of w = 10 individuals for each mixture of Scheffé’s matrix. This led to K = 30 elementary response matrices containing N = 286 and 3003 barycentric profiles in B and A, respectively. Finally, the 30 matrices were averaged to obtain a smoothed matrix integrating potentially the variability between and within groups (Figure 3h). From the smoothed matrix, regulation trends between morphometric variables were plotted (Figure 3i). Geometrical highlighting of regulatory growth trends Group-dependent relationships between paired morphometric variables were highlighted by projecting the weight values wj (0–10) of group j on the corresponding points of simplex space. Points with equal weights were statistically covered by a 95% confidence ellipse. The succession of the w + 1 (=11) weight ellipses from 0 to w = 10 provided a trajectory indicating how the two considered morphometric variables varied the one relatively to the other in favour of the considered fish group (stock unit). Comparison of Spx results to traditional methods Checking, validation, and highlighting advantages of Spx results were performed by comparison with results from other methods (Supplementary material S7–S14). Variation ranges analysis of growth regulation variables Box plots of relative morphometric variables initially helped to characterize spatial–temporal fish groups by different variation ranges and levels of growth regulations. Association analysis between geographical groups and growth regulation profiles DA was applied to separate fish groups by their profiles of relative morphometric variables taken in single or combined forms. Attributions of growth regulation profiles to fish groups were carried out using Bayesian rule allocating a profile to the group showing maximal assignment probability (Seber, 1984). Analysis of growth regulation trends and levels Signs and levels of growth regulation trends in different groups were checked by several statistical methods including non-parametric Spearman correlations (Spm), correlation-based PCA and MCA (Escofier and Pagès, 1991; Semmar 2011, 2013): Preliminary correlation analysis between relative levels of morphometric variables was invested by calculating Spm in different fish groups (Semmar, 2013). PCA was applied on relative levels of morphometric variables to obtain correlation charts highlighting characteristic regulation trends of groups (Semmar, 2011). MCA provides visualized trend analysis by considering variation levels of variables within fish groups. MCA was applied on a complete binary dataset where the variables were organized into four modalities corresponding to values less than first quartile Q1, between Q1 and median (Q2), between Q2 and third quartile (Q3), and more than Q3 (Escofier and Pagès, 1991). Results Traditional data analysis Morphometric measures of total length are summarized in Table 1. Table 1. Merluccius hubbsi: summarized morphometric measures (Lt) in different stocks and years. Years . Stock . Mean ± SD (mm) . [Min, max] (mm) . Sample size . 1968–1972 SE (21°–29°S) 166.34 ± 68.60 74 ≤ Lt ≤ 352 59 S (29°–34°S) 161.48 ± 63.60 66 ≤ Lt ≤ 430 172 2004 SE (21°–29°S) 310.83 ± 84.03 144 ≤ Lt ≤ 618 491 S (29°–34°S) 374.60 ± 73.20 246 ≤ Lt ≤ 501 25 Years . Stock . Mean ± SD (mm) . [Min, max] (mm) . Sample size . 1968–1972 SE (21°–29°S) 166.34 ± 68.60 74 ≤ Lt ≤ 352 59 S (29°–34°S) 161.48 ± 63.60 66 ≤ Lt ≤ 430 172 2004 SE (21°–29°S) 310.83 ± 84.03 144 ≤ Lt ≤ 618 491 S (29°–34°S) 374.60 ± 73.20 246 ≤ Lt ≤ 501 25 Open in new tab Table 1. Merluccius hubbsi: summarized morphometric measures (Lt) in different stocks and years. Years . Stock . Mean ± SD (mm) . [Min, max] (mm) . Sample size . 1968–1972 SE (21°–29°S) 166.34 ± 68.60 74 ≤ Lt ≤ 352 59 S (29°–34°S) 161.48 ± 63.60 66 ≤ Lt ≤ 430 172 2004 SE (21°–29°S) 310.83 ± 84.03 144 ≤ Lt ≤ 618 491 S (29°–34°S) 374.60 ± 73.20 246 ≤ Lt ≤ 501 25 Years . Stock . Mean ± SD (mm) . [Min, max] (mm) . Sample size . 1968–1972 SE (21°–29°S) 166.34 ± 68.60 74 ≤ Lt ≤ 352 59 S (29°–34°S) 161.48 ± 63.60 66 ≤ Lt ≤ 430 172 2004 SE (21°–29°S) 310.83 ± 84.03 144 ≤ Lt ≤ 618 491 S (29°–34°S) 374.60 ± 73.20 246 ≤ Lt ≤ 501 25 Open in new tab The size effect removal was necessary due to significant temporal variation in terms of total length (p < 0.05 in the three tests, Supplementary material S4). After size effect removal, data from both periods were joined (FPERMANOVA = 58.96, p = 0.629) for the subsequent analyses. Huxley models fitted to each stock did not show significant differences (Table 2), except for Pdd2, Lpspe, and Lpa due to the values of a coefficient (ANCOVA, p < 0.05). Table 2. Merluccius hubbsi: summary of models (pooled) from analysis of covariance between log-data of body measurements and total length (Lt) of Southeastern and Southern stocks Variables . ANCOVA . Mann–Whitney test . Constant . Lt . Stock . VIF . a . ±s.e. . b . ±s.e. . p-Value . Coefficient . ±s.e. . p-Value . U . p-Value . Hc 1.47260 0.02250 −0.002710 0.00969 0.780 0.001040 0.002090 0.618 1.48 186531.0 0.8029 Hpc 0.87750 0.01920 0.005060 0.00821 0.538 −0.002180 0.001770 0.219 1.34 201298.0 0.2306 Lbfd 1.35480 0.01470 −0.000380 0.00626 0.952 0.000170 0.001350 0.901 1.35 204250.5 0.9529 Lbsd 1.91500 0.00794 −0.001770 0.00339 0.601 0.000755 0.000730 0.301 1.34 206671.0 0.4004 Lba 1.91773 0.00755 −0.000440 0.00322 0.892 0.000179 0.000695 0.796 1.34 207094.0 0.5918 Lh 1.75498 0.00765 0.000330 0.00326 0.920 −0.000134 0.000705 0.849 1.35 208294.5 0.2681 Pdd1 1.77879 0.00643 0.002090 0.00274 0.447 −0.000911 0.000593 0.125 1.35 202483.0 0.2562 Pdd2 1.94936 0.00605 −0.005700 0.00258 0.027 0.002415 0.000554 <0.001 1.34 208052.0 0.0005 Lpo 1.29500 0.01020 0.000540 0.00434 0.900 −0.000227 0.000938 0.809 1.35 206634.0 0.6407 Lm 1.45569 0.00906 −0.001760 0.00387 0.650 0.000762 0.000836 0.362 1.35 207861.5 0.3470 Lppv 1.71126 0.00899 −0.002990 0.00383 0.436 0.001258 0.000827 0.129 1.35 209239.0 0.0845 Lppe 1.75729 0.00741 −0.002830 0.00316 0.371 0.001197 0.000683 0.080 1.35 209325.5 0.1322 Lpspe 1.95080 0.00748 −0.003200 0.00319 0.317 0.001377 0.000690 0.046 1.34 206592.0 0.1883 Lpa 1.95330 0.00736 −0.004740 0.00314 0.132 0.002017 0.000679 0.003 1.35 201406.0 0.0210 Ed 1.07250 0.01510 −0.001080 0.00645 0.867 0.000470 0.001390 0.736 1.34 205779.0 0.8077 Lpso 1.42701 0.00995 −0.000480 0.00424 0.909 0.000198 0.000917 0.829 1.35 205605.0 0.7490 Lpv 1.45800 0.01160 −0.000490 0.00497 0.922 0.000220 0.001070 0.837 1.34 203838.5 0.6403 Variables . ANCOVA . Mann–Whitney test . Constant . Lt . Stock . VIF . a . ±s.e. . b . ±s.e. . p-Value . Coefficient . ±s.e. . p-Value . U . p-Value . Hc 1.47260 0.02250 −0.002710 0.00969 0.780 0.001040 0.002090 0.618 1.48 186531.0 0.8029 Hpc 0.87750 0.01920 0.005060 0.00821 0.538 −0.002180 0.001770 0.219 1.34 201298.0 0.2306 Lbfd 1.35480 0.01470 −0.000380 0.00626 0.952 0.000170 0.001350 0.901 1.35 204250.5 0.9529 Lbsd 1.91500 0.00794 −0.001770 0.00339 0.601 0.000755 0.000730 0.301 1.34 206671.0 0.4004 Lba 1.91773 0.00755 −0.000440 0.00322 0.892 0.000179 0.000695 0.796 1.34 207094.0 0.5918 Lh 1.75498 0.00765 0.000330 0.00326 0.920 −0.000134 0.000705 0.849 1.35 208294.5 0.2681 Pdd1 1.77879 0.00643 0.002090 0.00274 0.447 −0.000911 0.000593 0.125 1.35 202483.0 0.2562 Pdd2 1.94936 0.00605 −0.005700 0.00258 0.027 0.002415 0.000554 <0.001 1.34 208052.0 0.0005 Lpo 1.29500 0.01020 0.000540 0.00434 0.900 −0.000227 0.000938 0.809 1.35 206634.0 0.6407 Lm 1.45569 0.00906 −0.001760 0.00387 0.650 0.000762 0.000836 0.362 1.35 207861.5 0.3470 Lppv 1.71126 0.00899 −0.002990 0.00383 0.436 0.001258 0.000827 0.129 1.35 209239.0 0.0845 Lppe 1.75729 0.00741 −0.002830 0.00316 0.371 0.001197 0.000683 0.080 1.35 209325.5 0.1322 Lpspe 1.95080 0.00748 −0.003200 0.00319 0.317 0.001377 0.000690 0.046 1.34 206592.0 0.1883 Lpa 1.95330 0.00736 −0.004740 0.00314 0.132 0.002017 0.000679 0.003 1.35 201406.0 0.0210 Ed 1.07250 0.01510 −0.001080 0.00645 0.867 0.000470 0.001390 0.736 1.34 205779.0 0.8077 Lpso 1.42701 0.00995 −0.000480 0.00424 0.909 0.000198 0.000917 0.829 1.35 205605.0 0.7490 Lpv 1.45800 0.01160 −0.000490 0.00497 0.922 0.000220 0.001070 0.837 1.34 203838.5 0.6403 The cut-off probability for significance of model parameters was p ≤ 0.05. In bold: significant values. s.e., standard error of the coefficient; VIF, variance inflation factor; U, Mann–Whitney statistics, p-value, probability. Open in new tab Table 2. Merluccius hubbsi: summary of models (pooled) from analysis of covariance between log-data of body measurements and total length (Lt) of Southeastern and Southern stocks Variables . ANCOVA . Mann–Whitney test . Constant . Lt . Stock . VIF . a . ±s.e. . b . ±s.e. . p-Value . Coefficient . ±s.e. . p-Value . U . p-Value . Hc 1.47260 0.02250 −0.002710 0.00969 0.780 0.001040 0.002090 0.618 1.48 186531.0 0.8029 Hpc 0.87750 0.01920 0.005060 0.00821 0.538 −0.002180 0.001770 0.219 1.34 201298.0 0.2306 Lbfd 1.35480 0.01470 −0.000380 0.00626 0.952 0.000170 0.001350 0.901 1.35 204250.5 0.9529 Lbsd 1.91500 0.00794 −0.001770 0.00339 0.601 0.000755 0.000730 0.301 1.34 206671.0 0.4004 Lba 1.91773 0.00755 −0.000440 0.00322 0.892 0.000179 0.000695 0.796 1.34 207094.0 0.5918 Lh 1.75498 0.00765 0.000330 0.00326 0.920 −0.000134 0.000705 0.849 1.35 208294.5 0.2681 Pdd1 1.77879 0.00643 0.002090 0.00274 0.447 −0.000911 0.000593 0.125 1.35 202483.0 0.2562 Pdd2 1.94936 0.00605 −0.005700 0.00258 0.027 0.002415 0.000554 <0.001 1.34 208052.0 0.0005 Lpo 1.29500 0.01020 0.000540 0.00434 0.900 −0.000227 0.000938 0.809 1.35 206634.0 0.6407 Lm 1.45569 0.00906 −0.001760 0.00387 0.650 0.000762 0.000836 0.362 1.35 207861.5 0.3470 Lppv 1.71126 0.00899 −0.002990 0.00383 0.436 0.001258 0.000827 0.129 1.35 209239.0 0.0845 Lppe 1.75729 0.00741 −0.002830 0.00316 0.371 0.001197 0.000683 0.080 1.35 209325.5 0.1322 Lpspe 1.95080 0.00748 −0.003200 0.00319 0.317 0.001377 0.000690 0.046 1.34 206592.0 0.1883 Lpa 1.95330 0.00736 −0.004740 0.00314 0.132 0.002017 0.000679 0.003 1.35 201406.0 0.0210 Ed 1.07250 0.01510 −0.001080 0.00645 0.867 0.000470 0.001390 0.736 1.34 205779.0 0.8077 Lpso 1.42701 0.00995 −0.000480 0.00424 0.909 0.000198 0.000917 0.829 1.35 205605.0 0.7490 Lpv 1.45800 0.01160 −0.000490 0.00497 0.922 0.000220 0.001070 0.837 1.34 203838.5 0.6403 Variables . ANCOVA . Mann–Whitney test . Constant . Lt . Stock . VIF . a . ±s.e. . b . ±s.e. . p-Value . Coefficient . ±s.e. . p-Value . U . p-Value . Hc 1.47260 0.02250 −0.002710 0.00969 0.780 0.001040 0.002090 0.618 1.48 186531.0 0.8029 Hpc 0.87750 0.01920 0.005060 0.00821 0.538 −0.002180 0.001770 0.219 1.34 201298.0 0.2306 Lbfd 1.35480 0.01470 −0.000380 0.00626 0.952 0.000170 0.001350 0.901 1.35 204250.5 0.9529 Lbsd 1.91500 0.00794 −0.001770 0.00339 0.601 0.000755 0.000730 0.301 1.34 206671.0 0.4004 Lba 1.91773 0.00755 −0.000440 0.00322 0.892 0.000179 0.000695 0.796 1.34 207094.0 0.5918 Lh 1.75498 0.00765 0.000330 0.00326 0.920 −0.000134 0.000705 0.849 1.35 208294.5 0.2681 Pdd1 1.77879 0.00643 0.002090 0.00274 0.447 −0.000911 0.000593 0.125 1.35 202483.0 0.2562 Pdd2 1.94936 0.00605 −0.005700 0.00258 0.027 0.002415 0.000554 <0.001 1.34 208052.0 0.0005 Lpo 1.29500 0.01020 0.000540 0.00434 0.900 −0.000227 0.000938 0.809 1.35 206634.0 0.6407 Lm 1.45569 0.00906 −0.001760 0.00387 0.650 0.000762 0.000836 0.362 1.35 207861.5 0.3470 Lppv 1.71126 0.00899 −0.002990 0.00383 0.436 0.001258 0.000827 0.129 1.35 209239.0 0.0845 Lppe 1.75729 0.00741 −0.002830 0.00316 0.371 0.001197 0.000683 0.080 1.35 209325.5 0.1322 Lpspe 1.95080 0.00748 −0.003200 0.00319 0.317 0.001377 0.000690 0.046 1.34 206592.0 0.1883 Lpa 1.95330 0.00736 −0.004740 0.00314 0.132 0.002017 0.000679 0.003 1.35 201406.0 0.0210 Ed 1.07250 0.01510 −0.001080 0.00645 0.867 0.000470 0.001390 0.736 1.34 205779.0 0.8077 Lpso 1.42701 0.00995 −0.000480 0.00424 0.909 0.000198 0.000917 0.829 1.35 205605.0 0.7490 Lpv 1.45800 0.01160 −0.000490 0.00497 0.922 0.000220 0.001070 0.837 1.34 203838.5 0.6403 The cut-off probability for significance of model parameters was p ≤ 0.05. In bold: significant values. s.e., standard error of the coefficient; VIF, variance inflation factor; U, Mann–Whitney statistics, p-value, probability. Open in new tab For Pdd2, the values for each stock were aSE = 1.95178 and aS = 1.94695; for Lpspe aSE = 1.95217 and aS = 1.94942; for Lpa aSE = 1.95532 and aS = 1.95128, revealing longer morphometric measures of M. hubbsi in Southeastern stock (SE) than Southern stock (S). Concerning b-coefficient, negative values revealed a decrease in the proportion of the measurements when Lt increased (Table 2). The models showed low level of multicollinearity (Variance Inflation Factors close to 1, the value of no correlation) and the residue analysis did not show any trend. Mann–Whitney test highlighted only Pdd2 and Lpa as measurements capable to distinguish stocks after size effect removal. The joint analysis of the measurements reinforced the difference between the SE and S stocks (FPERMANOVA = 3.843, p = 0.007), although it was difficult to distinguish clear patterns. In the DA, the correct assignments after cross-validation presented an overall coincidence of 58.96%, with 62% for SE stock and 50% for S stock, revealing a limited power for stock identification based on body measurements. Spx-simulation results Compared growth regulations between the two periods at overall population scale Spx provided smoothed regulation trends between morphometric variables (Figure 4a and b). Initially, these trends were analysed at overall population scale including all the fish groups of period A vs. period B. Figure 4. Open in new tabDownload slide Merluccius hubbsi: smoothed trends between morphometric variables of body perimeter (a) and front part (b). Dark and light greys: periods A and B, respectively. Numbers of points in scatter plots = 30030 in A (10 × 3003) and 2860 in B (10 × 286) due to 10 replications of smoothed matrix. Figure 4. Open in new tabDownload slide Merluccius hubbsi: smoothed trends between morphometric variables of body perimeter (a) and front part (b). Dark and light greys: periods A and B, respectively. Numbers of points in scatter plots = 30030 in A (10 × 3003) and 2860 in B (10 × 286) due to 10 replications of smoothed matrix. Period A showed larger variation ranges (diversification) than period B concerning growth regulations of perimeter body variables (Pdd2, Lpa, Lbsd, Lba; Figure 4a). This could indicate some perturbations in B due to overfishing activities. Periods A and B showed distinct variation spaces in (Pdd2 vs. Lpa) and (Pdd2 vs. Lba) leading to strong differentiation between the two temporal populations. Positive trends concerned the pairs (Pdd2, Lpa) and (Lbsd, Lba) whereas negative trends opposed (Pdd2, Lpa) (front body variables) to (Lba, Lbsd) (back body variabes), especially in A. Moreover, in B, Lbsd showed a loss of monotonicity towards the other variables. At lower body scale, regulation trends between front body variables (Lbfd, Lpso, Ed, Lpo, Lppv) showed distinct variation spaces separating periods A and B (Figure 4b): higher regulation levels concerned Lbfd, Ed in A vs. Lpo, Lpso, Lppv in B. This highlighted significant increase of growth regulations of mouth part in B. Functional differentiations between fish groups Spx highlighted group-dependent trends making morphometric variables to reach different balance levels through different variation ways and relational shapes according to fish groups (Figure 5). The pairs (Pdd2, Lpa) and (Lpso, Lppv) were considered to illustrate such growth flexibility through body perimeter and within front compartment, respectively. Figure 5. Open in new tabDownload slide Merluccius hubbsi: smoothed plots showing different trends (bold trajectories) and local variations (weight ellipses’ inclinations) governing regulations of morphometric variables Lpa and Pdd2 in different fish groups associated with different geographical sites: A1 (a), A2 (b), A3 (c), A4 (d), A5 (e), A6 (f) in 1968–1972, B1 (g), B2 (h), B3 (i), and B5 (j) in 2004. The numbers of smoothed points in the scatter plots were equal to 30 030 in A (10 × 3003) and 2860 in B (10 × 286) due to 10 replications of smoothed matrix. Figure 5. Open in new tabDownload slide Merluccius hubbsi: smoothed plots showing different trends (bold trajectories) and local variations (weight ellipses’ inclinations) governing regulations of morphometric variables Lpa and Pdd2 in different fish groups associated with different geographical sites: A1 (a), A2 (b), A3 (c), A4 (d), A5 (e), A6 (f) in 1968–1972, B1 (g), B2 (h), B3 (i), and B5 (j) in 2004. The numbers of smoothed points in the scatter plots were equal to 30 030 in A (10 × 3003) and 2860 in B (10 × 286) due to 10 replications of smoothed matrix. Figure 5. Open in new tabDownload slide continued Figure 5. Open in new tabDownload slide continued Pdd2 vs. Lpa showed positive global trends resulting from monotonic succession of different full-weight ellipses (w = 10) associated with different fish groups (Figure 5). This indicated that perimeter front variables (Pdd2, Lpa) were similarly influenced by environmental conditions associated with different geographical sites: growth regulation levels of both Pdd2 and Lpa were high in groups A2, A4, B3, B5, low in A1, A6, B1, intermediate in A3, A5, B2. These different levels highlighted diversified/heterogeneous aspects within conventional stocks (I and II) due to different growth regulation mechanisms between neighbour fish groups (A1 vs. A2 vs. A3), (B1 vs. B2 vs. B3) (for stock I) and (A4 vs. A5 vs. A6) (for stock II). Such intra-stock diversity showed alternation of high and intermediate growth regulation levels in middle sites vs. minimal states in extreme sites (A1, B1, A6). Positive global trends were also highlighted for back variables (Lbsd, Lba) which were opposite to (Pdd2, Lpa) Supplementary material S6a, b, e, and f). This could be directly due to two different mechanistic needs in fish: eating vs moving; it could also indicate opposite effects of environment on growth regulatory ways in body front and back compartments. At local (within-group) scale, Pdd2 vs. Lpa showed different inclinations of full-weight ellipses indicating flexible correlations between same variables (Figure 5): correlations were negative in A3, A4, A6, B5, slightly positive in A1, B1, and not significant in A2, A5, B2, B3. Negative local trends could be indicative of competition between top (Pdd2) and bottom (Lpa) body parts for biomass distribution in front compartment (A3, A4, A6, B5). However, positive trends indicated some cooperation mechanism between Pdd2 and Lpa in favour of their common front compartment (in A1, B1). Local competing and local cooperative mechanisms were also highlighted for the pair (Lbsd, Lba) of back compartment with negative trends in A4, A6 vs. positive in A2, A5 (Supplementary material S6d and f vs S6b and e). Further differentiations of fish groups were provided by analysing growth regulation trends between lower constituents of body front compartment (illustrated by Lppv vs. Lpso; Figure 6): in A, Lpso and Lppv showed strong negative trends at both global and local scales (Figure 6a–f). This indicated growth differentiation processes based on opposite regulations between top and bottom body segments. Extreme cases concerned A2 and A3 where (Lpso, Lppv) showed (maximal, minimal) and (minimal, maximal) states, respectively (Figure 6b and c). In B-fish groups, extreme global states concerned B3 and B5 with (minimal, maximal) and (maximal, maximal) regulation levels of (Lpso, Lppv), respectively (Figure 6i and j). Double maximum in B5 indicated a trend for bigger mouth development favoured by some environment conditions. However, locally, negative inclinations of full-weight ellipses indicated well-conserved intra-compartment competition for biomass distributions between top (Lpso) and bottom (Lppv) segments within the body front compartment (Figure 6j). Figure 6. Open in new tabDownload slide Merluccius hubbsi: smoothed plots showing different trends (bold trajectories) and local variations (weight ellipses’ inclinations) governing regulations of morphometric variables Lppv and Lpso in different fish groups associated with different geographical sites A1 (a), A2 (b), A3 (c), A4 (d), A5 (e), A6 (f) in 1968–1972, B1 (g), B2 (h), B3 (i), and B5 (j) in 2004. The numbers of smoothed points in the scatter plots were equal to 30 030 in A (10 × 3003) and 2860 in B (10 × 286) due to 10 replications of smoothed matrix. Figure 6. Open in new tabDownload slide Merluccius hubbsi: smoothed plots showing different trends (bold trajectories) and local variations (weight ellipses’ inclinations) governing regulations of morphometric variables Lppv and Lpso in different fish groups associated with different geographical sites A1 (a), A2 (b), A3 (c), A4 (d), A5 (e), A6 (f) in 1968–1972, B1 (g), B2 (h), B3 (i), and B5 (j) in 2004. The numbers of smoothed points in the scatter plots were equal to 30 030 in A (10 × 3003) and 2860 in B (10 × 286) due to 10 replications of smoothed matrix. Spatial-temporal monitoring of average growth regulation states of stock units Smoothed average states of body perimeter variables showed higher variability in period A than in B (Figure 7a). This could indicate higher resiliency of fish groups in A. Stock units were more differentiated by growth regulation levels than variation ranges: A4–A6 (stock II in 1968–1972) showed high Lba vs. low Lpa whereas B1–B3 (stock I in 2004) showed opposite aspect in favour of bigger mouth trend. A1–A3 (stock I in 1968–1972) had intermediate aspect with particularly extreme state in A2 (low Lbsd, Lba vs. high Pdd2, Lpa) making this fish group to be well distinguished. Opposition between A4–A6 and B1–B3 concerned also morphometric variables of front compartment (Figure 7b): high Lbfd, Lpso, Lpo vs. low Ed in B1–B3 with opposite trends in A4–A6. This highlighted clear effects of both geographical location and time to differentiate stock units. Figure 7. Open in new tabDownload slide Merluccius hubbsi: variation of average states of morphometric variables corresponding to full-weight ellipses in different smoothed relationships by Spx. (a) Body perimeter variables and (b) front compartment variables. A1–A6, fish groups of period A (1968–1972); B1–B3 and B5, fish groups of period B (2004). Figure 7. Open in new tabDownload slide Merluccius hubbsi: variation of average states of morphometric variables corresponding to full-weight ellipses in different smoothed relationships by Spx. (a) Body perimeter variables and (b) front compartment variables. A1–A6, fish groups of period A (1968–1972); B1–B3 and B5, fish groups of period B (2004). Successive geographical areas resulted in target or modulated regulations of growths: for instance, in Ed vs. Lbfd, Ed increased from A1 to A6 whereas Lbfd showed alternated (cyclic) variations (Figure 7b). Target and modulated aspects were also highlighted for Lba vs. Lbsd, respectively (A2–A6; Figure 7a). Comparison between Spx results and other methods Variation ranges analysis of growth regulation variables Box plots showed different average relative levels (growth regulations) of morphometric variables in different fish groups (Supplementary material S7; Table 3). Mann–Whitney test highlighted significant differences between the two periods A and B concerning Lba, Ed (higher in A), Lpa, Lpso, Lpo, Lppv (higher in B; p < 10−4; Supplementary material S7). This agreed with Spx results showing clear variations of relative growth investments in fish groups with time (Figures 4 and 7). Also, box plots showed characteristic growth regulation levels of fish groups that agreed with Spx results (Table 3; Figure 7; Supplementary material S7): low Lpa (A1), low Lbsd vs. maximal Lpa (A2), maximal Lbfd vs. low Lpo (A3), maximal Pdd2 vs. minimal Lbsd (A4), maximal Ed vs. minimal Lpo (A5), minimal Pdd2, Lpa, Lpo vs. maximal Lbsd, Lba (A6). In B, extreme growth regulation levels essentially concerned site B5 showing minimal Lba, Lbfd, Ed vs. maximal Lpa, Lpso, Lpo, Lppv. Table 3. Merluccius hubbsi: means and standard deviations of regulation levels of different morphometric variables of body perimeter (Lbsd, Lba, Pdd2, Lpa) and front compartment (Lbfd, Lpso, Ed, Lpo, Lppv) in fish groups of periods A (A1–A6) and B (B1–B3, B5). Strata . Lbsd . Lba . Pdd2 . Lpa . Lbfd . Lpso . Ed . Lpo . Lppv . A 0.241 ± 0.009 0.246 ± 0.008 0.256 ± 0.007 0.257 ± 0.009 0.172 ± 0.011 0.198 ± 0.009 0.098 ± 0.011 0.149 ± 0.008 0.383 ± 0.012 B 0.241 ± 0.006 0.243 ± 0.006 0.255 ± 0.005 0.261 ± 0.007 0.172 ± 0.009 0.208 ± 0.008 0.082 ± 0.008 0.152 ± 0.005 0.386 ± 0.01 A1 0.242 ± 0.008 0.244 ± 0.006 0.256 ± 0.006 0.258 ± 0.007 0.17 ± 0.012 0.198 ± 0.011 0.099 ± 0.012 0.149 ± 0.005 0.383 ± 0.011 A2 0.237 ± 0.007 0.243 ± 0.008 0.258 ± 0.006 0.262 ± 0.012 0.174 ± 0.011 0.202 ± 0.009 0.094 ± 0.012 0.149 ± 0.009 0.381 ± 0.013 A3 0.240 ± 0.009 0.244 ± 0.008 0.258 ± 0.01 0.258 ± 0.009 0.176 ± 0.014 0.195 ± 0.01 0.097 ± 0.009 0.148 ± 0.006 0.384 ± 0.011 A4 0.236 ± 0.01 0.244 ± 0.008 0.259 ± 0.006 0.261 ± 0.009 0.172 ± 0.012 0.198 ± 0.008 0.096 ± 0.011 0.150 ± 0.015 0.383 ± 0.016 A5 0.240 ± 0.009 0.244 ± 0.008 0.257 ± 0.006 0.259 ± 0.009 0.171 ± 0.011 0.197 ± 0.012 0.101 ± 0.012 0.148 ± 0.006 0.383 ± 0.012 A6 0.244 ± 0.009 0.248 ± 0.007 0.253 ± 0.007 0.255 ± 0.008 0.172 ± 0.01 0.198 ± 0.008 0.098 ± 0.011 0.148 ± 0.005 0.383 ± 0.01 B1 0.241 ± 0.006 0.244 ± 0.006 0.254 ± 0.005 0.26 ± 0.006 0.173 ± 0.009 0.209 ± 0.008 0.08 ± 0.008 0.152 ± 0.005 0.385 ± 0.009 B2 0.240 ± 0.006 0.243 ± 0.006 0.256 ± 0.006 0.261 ± 0.007 0.173 ± 0.01 0.209 ± 0.007 0.082 ± 0.008 0.152 ± 0.005 0.384 ± 0.009 B3 0.241 ± 0.006 0.241 ± 0.006 0.256 ± 0.005 0.262 ± 0.008 0.170 ± 0.009 0.206 ± 0.009 0.084 ± 0.008 0.151 ± 0.006 0.388 ± 0.01 B5 0.241 ± 0.004 0.240 ± 0.005 0.256 ± 0.004 0.264 ± 0.006 0.168 ± 0.009 0.211 ± 0.007 0.077 ± 0.005 0.155 ± 0.005 0.389 ± 0.01 Strata . Lbsd . Lba . Pdd2 . Lpa . Lbfd . Lpso . Ed . Lpo . Lppv . A 0.241 ± 0.009 0.246 ± 0.008 0.256 ± 0.007 0.257 ± 0.009 0.172 ± 0.011 0.198 ± 0.009 0.098 ± 0.011 0.149 ± 0.008 0.383 ± 0.012 B 0.241 ± 0.006 0.243 ± 0.006 0.255 ± 0.005 0.261 ± 0.007 0.172 ± 0.009 0.208 ± 0.008 0.082 ± 0.008 0.152 ± 0.005 0.386 ± 0.01 A1 0.242 ± 0.008 0.244 ± 0.006 0.256 ± 0.006 0.258 ± 0.007 0.17 ± 0.012 0.198 ± 0.011 0.099 ± 0.012 0.149 ± 0.005 0.383 ± 0.011 A2 0.237 ± 0.007 0.243 ± 0.008 0.258 ± 0.006 0.262 ± 0.012 0.174 ± 0.011 0.202 ± 0.009 0.094 ± 0.012 0.149 ± 0.009 0.381 ± 0.013 A3 0.240 ± 0.009 0.244 ± 0.008 0.258 ± 0.01 0.258 ± 0.009 0.176 ± 0.014 0.195 ± 0.01 0.097 ± 0.009 0.148 ± 0.006 0.384 ± 0.011 A4 0.236 ± 0.01 0.244 ± 0.008 0.259 ± 0.006 0.261 ± 0.009 0.172 ± 0.012 0.198 ± 0.008 0.096 ± 0.011 0.150 ± 0.015 0.383 ± 0.016 A5 0.240 ± 0.009 0.244 ± 0.008 0.257 ± 0.006 0.259 ± 0.009 0.171 ± 0.011 0.197 ± 0.012 0.101 ± 0.012 0.148 ± 0.006 0.383 ± 0.012 A6 0.244 ± 0.009 0.248 ± 0.007 0.253 ± 0.007 0.255 ± 0.008 0.172 ± 0.01 0.198 ± 0.008 0.098 ± 0.011 0.148 ± 0.005 0.383 ± 0.01 B1 0.241 ± 0.006 0.244 ± 0.006 0.254 ± 0.005 0.26 ± 0.006 0.173 ± 0.009 0.209 ± 0.008 0.08 ± 0.008 0.152 ± 0.005 0.385 ± 0.009 B2 0.240 ± 0.006 0.243 ± 0.006 0.256 ± 0.006 0.261 ± 0.007 0.173 ± 0.01 0.209 ± 0.007 0.082 ± 0.008 0.152 ± 0.005 0.384 ± 0.009 B3 0.241 ± 0.006 0.241 ± 0.006 0.256 ± 0.005 0.262 ± 0.008 0.170 ± 0.009 0.206 ± 0.009 0.084 ± 0.008 0.151 ± 0.006 0.388 ± 0.01 B5 0.241 ± 0.004 0.240 ± 0.005 0.256 ± 0.004 0.264 ± 0.006 0.168 ± 0.009 0.211 ± 0.007 0.077 ± 0.005 0.155 ± 0.005 0.389 ± 0.01 Black and grey highlighting indicated maximal and minimal levels of different systems. Open in new tab Table 3. Merluccius hubbsi: means and standard deviations of regulation levels of different morphometric variables of body perimeter (Lbsd, Lba, Pdd2, Lpa) and front compartment (Lbfd, Lpso, Ed, Lpo, Lppv) in fish groups of periods A (A1–A6) and B (B1–B3, B5). Strata . Lbsd . Lba . Pdd2 . Lpa . Lbfd . Lpso . Ed . Lpo . Lppv . A 0.241 ± 0.009 0.246 ± 0.008 0.256 ± 0.007 0.257 ± 0.009 0.172 ± 0.011 0.198 ± 0.009 0.098 ± 0.011 0.149 ± 0.008 0.383 ± 0.012 B 0.241 ± 0.006 0.243 ± 0.006 0.255 ± 0.005 0.261 ± 0.007 0.172 ± 0.009 0.208 ± 0.008 0.082 ± 0.008 0.152 ± 0.005 0.386 ± 0.01 A1 0.242 ± 0.008 0.244 ± 0.006 0.256 ± 0.006 0.258 ± 0.007 0.17 ± 0.012 0.198 ± 0.011 0.099 ± 0.012 0.149 ± 0.005 0.383 ± 0.011 A2 0.237 ± 0.007 0.243 ± 0.008 0.258 ± 0.006 0.262 ± 0.012 0.174 ± 0.011 0.202 ± 0.009 0.094 ± 0.012 0.149 ± 0.009 0.381 ± 0.013 A3 0.240 ± 0.009 0.244 ± 0.008 0.258 ± 0.01 0.258 ± 0.009 0.176 ± 0.014 0.195 ± 0.01 0.097 ± 0.009 0.148 ± 0.006 0.384 ± 0.011 A4 0.236 ± 0.01 0.244 ± 0.008 0.259 ± 0.006 0.261 ± 0.009 0.172 ± 0.012 0.198 ± 0.008 0.096 ± 0.011 0.150 ± 0.015 0.383 ± 0.016 A5 0.240 ± 0.009 0.244 ± 0.008 0.257 ± 0.006 0.259 ± 0.009 0.171 ± 0.011 0.197 ± 0.012 0.101 ± 0.012 0.148 ± 0.006 0.383 ± 0.012 A6 0.244 ± 0.009 0.248 ± 0.007 0.253 ± 0.007 0.255 ± 0.008 0.172 ± 0.01 0.198 ± 0.008 0.098 ± 0.011 0.148 ± 0.005 0.383 ± 0.01 B1 0.241 ± 0.006 0.244 ± 0.006 0.254 ± 0.005 0.26 ± 0.006 0.173 ± 0.009 0.209 ± 0.008 0.08 ± 0.008 0.152 ± 0.005 0.385 ± 0.009 B2 0.240 ± 0.006 0.243 ± 0.006 0.256 ± 0.006 0.261 ± 0.007 0.173 ± 0.01 0.209 ± 0.007 0.082 ± 0.008 0.152 ± 0.005 0.384 ± 0.009 B3 0.241 ± 0.006 0.241 ± 0.006 0.256 ± 0.005 0.262 ± 0.008 0.170 ± 0.009 0.206 ± 0.009 0.084 ± 0.008 0.151 ± 0.006 0.388 ± 0.01 B5 0.241 ± 0.004 0.240 ± 0.005 0.256 ± 0.004 0.264 ± 0.006 0.168 ± 0.009 0.211 ± 0.007 0.077 ± 0.005 0.155 ± 0.005 0.389 ± 0.01 Strata . Lbsd . Lba . Pdd2 . Lpa . Lbfd . Lpso . Ed . Lpo . Lppv . A 0.241 ± 0.009 0.246 ± 0.008 0.256 ± 0.007 0.257 ± 0.009 0.172 ± 0.011 0.198 ± 0.009 0.098 ± 0.011 0.149 ± 0.008 0.383 ± 0.012 B 0.241 ± 0.006 0.243 ± 0.006 0.255 ± 0.005 0.261 ± 0.007 0.172 ± 0.009 0.208 ± 0.008 0.082 ± 0.008 0.152 ± 0.005 0.386 ± 0.01 A1 0.242 ± 0.008 0.244 ± 0.006 0.256 ± 0.006 0.258 ± 0.007 0.17 ± 0.012 0.198 ± 0.011 0.099 ± 0.012 0.149 ± 0.005 0.383 ± 0.011 A2 0.237 ± 0.007 0.243 ± 0.008 0.258 ± 0.006 0.262 ± 0.012 0.174 ± 0.011 0.202 ± 0.009 0.094 ± 0.012 0.149 ± 0.009 0.381 ± 0.013 A3 0.240 ± 0.009 0.244 ± 0.008 0.258 ± 0.01 0.258 ± 0.009 0.176 ± 0.014 0.195 ± 0.01 0.097 ± 0.009 0.148 ± 0.006 0.384 ± 0.011 A4 0.236 ± 0.01 0.244 ± 0.008 0.259 ± 0.006 0.261 ± 0.009 0.172 ± 0.012 0.198 ± 0.008 0.096 ± 0.011 0.150 ± 0.015 0.383 ± 0.016 A5 0.240 ± 0.009 0.244 ± 0.008 0.257 ± 0.006 0.259 ± 0.009 0.171 ± 0.011 0.197 ± 0.012 0.101 ± 0.012 0.148 ± 0.006 0.383 ± 0.012 A6 0.244 ± 0.009 0.248 ± 0.007 0.253 ± 0.007 0.255 ± 0.008 0.172 ± 0.01 0.198 ± 0.008 0.098 ± 0.011 0.148 ± 0.005 0.383 ± 0.01 B1 0.241 ± 0.006 0.244 ± 0.006 0.254 ± 0.005 0.26 ± 0.006 0.173 ± 0.009 0.209 ± 0.008 0.08 ± 0.008 0.152 ± 0.005 0.385 ± 0.009 B2 0.240 ± 0.006 0.243 ± 0.006 0.256 ± 0.006 0.261 ± 0.007 0.173 ± 0.01 0.209 ± 0.007 0.082 ± 0.008 0.152 ± 0.005 0.384 ± 0.009 B3 0.241 ± 0.006 0.241 ± 0.006 0.256 ± 0.005 0.262 ± 0.008 0.170 ± 0.009 0.206 ± 0.009 0.084 ± 0.008 0.151 ± 0.006 0.388 ± 0.01 B5 0.241 ± 0.004 0.240 ± 0.005 0.256 ± 0.004 0.264 ± 0.006 0.168 ± 0.009 0.211 ± 0.007 0.077 ± 0.005 0.155 ± 0.005 0.389 ± 0.01 Black and grey highlighting indicated maximal and minimal levels of different systems. Open in new tab Association analysis between growth regulation profiles and geographical groups DA used morphometric variables with single and/or combined (interactive) forms to provide predictive models of the ten fish groups (A1–A6, B1–B3, B5; Supplementary material S8). The percentages of corrected predictions varied between 65% (A5) and 83% (B5) (Supplementary material S9). Discriminant variables showed extreme average levels in associated groups (Supplementary material S7 and S9; Table 3): A1, Lbfd (minimum average level); A2, Lppv (Min), Lbfd (high); A3, Pdd2 (Max), Lpso (Min); A4, Lbsd (Min), Lpso (Low); A5, Ed (Max), Lpso (Min); A6, Lba, Lbsd (Max), Lpa, Pdd2 (Min); B3, Lpso, Lppv (Max), Lpa (high), Ed, Lba (low); B5, Lpo (Max), Ed (Min). These variables showed extreme locations of full-weight ellipses of corresponding fish groups in Spx (Figure 7). Although DA models highlighted differentiation roles of morphometric variables, the development ways of groups remain unanswered question. Analysis of trends between growth regulation variables Spm negative correlations corresponded to negatively inclined full-weight ellipses in Spx (Supplementary material S10). Negative trends mainly concerned (Lpa, Lbsd), (Pdd2, Lba), (Lbfd, Lppv), (Ed, Lpso) in periods A and B followed by (Pdd2, Lbsd), (Lpa, Lba), (Lpo, Lppv), (Lpso, Lppv) in B (Spm) and both periods (Spx). These negative trends were also highlighted by PCA through opposite projections of variable points (Supplementary material S11): (Lpa, Lbsd) in A1, A3, A5, B3; (Lba, Lpa) in A2, B1, B2, B5; (Pdd2, Lba) in A1, A5; (Ed, Lpso) in A2, A5, A6, B1, B2, B3; (Lppv, Lbfd) in A2, A3, A6, B2, B3, B5; (Pdd2, Lpa) in A4. Spm positive trends corresponded to positively inclined full-weight ellipses in Spx (Supplementary material S10). The most concerned pairwise was (Lpo, Lpso) in A2, A3, B1, B3, B5 (Spx) vs. B1, B2, B3 (Spm). It was followed by (Ed, Lppv) in A3 (Spx) vs. A3, B1, B3 (Spm), (Lpso, Lbfd) in A1, A2, A4 (Spx) vs. A4 (Spm), (Pdd2, Lpa) in A1, B1 (Spx). These positive correlations were confirmed in PCA by variable points showing close projections in factorial plots (Supplementary material S11). Also, Spx and Spm agreed concerning not detected (not significant) trends indicated by high p-values in Spm and by circular or not inclined full-weight ellipses in Spx (Supplementary material S10): e.g. Lba vs. Lbsd in A1, A3, B1–B3, B5 (Supplementary material S6a, c, and g–j). This could reveal non-linear dependences between the variables. Compared with Spm and PCA, Spx advantageously provided visualization of variation spaces of group-dependent trends. Such group-depending associations between correlation signs and variation ranges of variables were checked by MCA leading to further validation of Spx results. Multiple correspondence analysis MCA highlighted strong separation between fish groups of periods A and B due to differences in growth regulation of front body variables (Supplementary material S12a). Period B showed high regulations of Lpso and Lpo (Lpso4, Lpo4) contrary to Ed (Ed1) vs. opposite aspect in A (Lpso1, Lpo1, Ed4; Supplementary material S12b). These results agreed with those of Spx which highlighted strong separation between groups A1–A6 and B1–B3, B5 within the variation space of balance average regulations Lpso vs. Ed (Supplementary material S12c and d). MCA applied on body perimeter variables (Pdd2, Lpa, Lbsd, Lba) highlighted clear differentiations between fish groups A1–A6 (Supplementary material S13a and b). Results highlighted extreme states of A6, A4, and A2: high regulations of Lba, Lbsd vs. low levels of Lpa, Pdd2 in A6 (Lbsd4, Lba4, Lpa1, Pdd2.1); high regulations of Lpa, Pdd2 vs. low levels of Lbsd in A4, A2 (Lpa4, Pdd2.4, Lbsd1) and Lba in A2 (Lba1). A1 showed intermediate states between A6 and A2, A4 (Pdd2.2, Pdd2.3, Lbsd2, Lbsd3, Lba3, Lpa3). All these results agreed with those of Spx (Supplementary material S13c). Concerning intermediate states of A3, A5 highlighted by Spx, they were revealed by MCA by widely dispersed points along F1 for A5 and F2 for A3. ACM results highlighted some heterogeneity aspects of A3, A5 compared with other groups. MCA applied on constitutive parts of front compartment (Lbfd, Lpso, Lppv, Lpo, Ed) highlighted differentiations between the four fish groups (B1–B3, B5) in 2004 (Supplementary material S14a and b). The four modalities of different variables were projected in different factorial subspaces that were associated with different fish groups: Lpso and Lpo increased from B3 to B5 via B2 and B1, respectively; Ed showed opposite trend. Regulation levels of Lbfd increased from B5 (Lbfd1) to B1 (Lbfd4) via B3 and B2, respectively; Lppv showed opposite trend. These results were highlighted by Spx which also provided the balance levels of variables the ones relatively to the others (Supplementary material S14c). Discussion Spx provided an original simulation method of growth regulation trends between body constitutive variables helping to understand better mechanistic origins of population diversification and differentiation ways of biological groups. Based on population stratification and groups’ weighting (combinations), Spx has the advantage to treat heterogeneous systems containing several groups characterized by unlimited number of variables whatever their variability levels could be. Spx highlighted different scales mechanisms governing complex growth variability between and within several fish population groups: global trends governed growth regulation levels and differentiations between groups (under different overall/environmental conditions). However, local trends highlighted intrinsic growth regulation ways within groups leading to competing or cooperating body parts for local biomass distribution. By this way, Spx provides mechanistic information helping to define/identify stocks with specific growth regulation ways of different body parts. This provides highlighting of characteristic processes responsible for self-sustainability and resiliency of population groups in responses to environmental conditions and fishing activities. Six fish groups coming from six geographical areas in 1968–1972 were well discriminated by well distinct regulation levels of body perimeter variables including Pdd2, Lpa in front compartment and Lbsd, Lba in back compartment (Figure 5; Supplementary material S6a and c). This indicated a strong order of biomass distribution in fish groups at whole body scale due to intense growth rates of juveniles (Vaz-dos-Santos and Rossi-Wongtschowski, 2007; Costa et al., 2018). However, in 2004, regulation ratios among most hierarchical variables (Pdd2, Lpa, Lbsd, Lba) showed narrower variation ranges indicating more homogeneous growth of adults (Figure 4a). This could result from lower intraspecific competition for foraging due to overfishing: overfishing leads to a reduction of stock abundance favouring relatively high prey abundance and subsequently higher capture per individual predator (M. hubbsi) (Perez et al., 2009; Vaz-dos-Santos et al., 2010). Moreover, higher and lower variability in A and B, respectively, seemed to be associated with two development strategies in fish populations (Minto et al., 2008; Panikian et al., 2015): higher resiliency in A vs. higher viability in B due to overfishing resulting in lower intraspecific competition for feeding. At low body scale, Lpso, Lppv, and Lpo showed frankly higher regulation levels in 2004 than 1968–1972 indicating strong positive trend within body front compartment during the three decades (Figure 7b). Considering the relatively higher food availability (prey abundances) in 2004 (Muto and Soares, 2011), bigger jaws could be developed for higher capture-yielding in M. hubbsi vs. a decrease of relative growth of back body part (Lba; less required moving for prey capture; Figure 7a and b). This provides mechanistic argument on the fact that overexploited fishes tend to grow faster with target biomass distribution in body (Hart and Reynolds, 2002). Spx results were compatible with regional hydrological phenomena: in the period of 1968–1972, A3 was distinct from the Southeast A1 and South extremes (A5, A6) by extreme growth regulations of front body parts (more particularly Lbfd; Figure 7b). This result reinforced the key role of the upwelling systems and the Patos-Mirim and La Plata systems (Piola et al., 2018). Moreover, the sequence and closeness of A4, A5, and A6 confirms that M. hubbsi is shared with Uruguay and Argentina (Vaz-dos-Santos et al., 2017). Spx shows several methodological strengths and application perspectives for stock identification, monitoring, and management: Fish groups can be openly characterized by regulation trends of several types of constitutive body variables including length, weight, feeding, body composition, otoliths, etc. Moreover, population stratification admits flexible criteria that can be intrinsic (physiological states, maturity degrees, ages, body length ranges, population sizes, stock biomasses) or extrinsic (environmental conditions, fishing activities). Therefore, simulated growth regulation trends could provide helpful information on gradual, transitory, or chained mechanisms implied in variability, persistency, adaptive strategies, and/or stock recruitment relationships (Minto et al., 2008; Panikian et al., 2015). However, higher number of stratification groups q can be limiting because of resulting drastic increase of total number of combinations (N). Although Spx does not provide predictive formulation, it has the ultimate advantage to highlight multidirectional, multi-shape, and multi-scale trends (group-dependent trends) preparing for next predictive modelling between regulation variables. Application perspectives of Spx extend beyond fishpopulations. Among modern approaches, DNA barcoding provides reliable background reference libraries from which unidentified states or specimens can be delineated, classified, or precisely recognized (Hubert and Hanner, 2015). DNA barcodes can be used as categorical variables for population stratification in Spx. This helps to highlight regulation trends associating genomic categories of fish populations with different body constitution patterns (length, weight, feeding, metabolism, etc.). Critics were given about working with preserved and fixed specimens leading to shrinkage bias on morphometric measures. Such a bias is strongly reduced by Spx because of three methodological factors: (i) initial working on ratio variables followed by (ii) pattern averaging then (iii) iteration of complete set of average patterns. Working on ratio variables (X/Y) leads to variation space with considerably lower variance and coefficient of variation compared to measured variable X (Supplementary material S15). Using morphometric data initially increased with random errors of 5–10%, simulation revealed stable regulation trends indicating strong attenuation of shrinkage bias with relative variables (Supplementary Material S16). Finally, mathematical demonstration based on error theory was given on bias reduction associated with ratio variable (Supplementary Material S17). Supplementary data Supplementary material is available at the ICESJMS online version of the manuscript. Acknowledgements NS and AMVS thank the Brazilian National Council for the Scientific and Technological Development (CNPq) for the financial support (Process 453871/2016-0) and research grant (Process 310451/2018-3). AMVS is grateful to Professors Carmen Lúcia Del Bianco Rossi-Wongtschowski, José Lima de Figueiredo, Naércio Aquino Menezes, and Francisco Manoel de Souza Braga, whose discussions on the theme in the middle of the 2000s were valuable. References Anderson M. J. , Gorley R. N., Clarke K. R. 2008 . PERMANOVA+ for PRIMER: Guide to Software and Statistical Methods. PRIMER-E Ltd, Plymouth. Cadrin S. X. 2000 . Advances in morphometric identification of fishery stocks . Reviews in Fish Biology and Fisheries , 10 : 91 – 112 . Google Scholar Crossref Search ADS WorldCat Cadrin S. X. , Kerr L. A., Mariani S. (Eds). 2014 . Stock Identification Methods: Applications in Fishery Science , 2nd edn. Elsevier , Oxford . 588 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Cornell J. A. 2002 . Experiments with Mixtures: Designs, Models and the Analysis of Mixture Data , 2nd edn. John Willey & Son , New York . 649 pp. Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Cornell J. A. 2011 . A retrospective view of mixture experiments . Quality Engineering , 23 : 315 – 331 . Google Scholar Crossref Search ADS WorldCat Costa P. A. S. , Braga A., da C., Vieira J. M., da S., Martins R. R. M., de São-Clemente R. R. B., Couto B. R. 2018 . Age estimation, growth and maturity of the Argentine hake (Merluccius hubbsi Marini, 1933) along the northernmost limit of its distribution in the south-western Atlantic. Marine Biology Research , 14 : 728 – 738 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Escofier B. , Pagès J. 1991 . Presentation of correspondence analysis and multiple correspondence analysis with the help of examples. In Applied Multivariate Analysis in SAR and Environmental Studies , pp. 1 – 32 . Ed. by Devillers J., Karcher W.. Kluwer Academic Publishers , Dordrecht . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Gotelli N. J. , Ellison A. M. 2004 . A Primer of Ecological Statistics . Sinauer Associates , Sunderland . 510 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hart P. J. B. , Reynolds J. D. 2002 . Handbook of Fish Biology and Fisheries, Volume 1: Fish Biology. Blackwell Publishing , Malden . 413 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hubert N. , Hanner R. 2015 . DNA barcoding, species delineation and taxonomy: a historical perspective . DNA Barcodes , 3 : 44 – 58 . OpenURL Placeholder Text WorldCat Huxley J. S. 1924 . Constant differential growth-ratios and their significance . Nature , 114 : 895 – 896 . Google Scholar Crossref Search ADS WorldCat Huxley J. S. 1932 . Problems of Relative Growth. Metheun , London . 319 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Inada T. 1981 . Studies on the merlucciid fishes . Bulletin of the Far Seas Fisheries Research Laboratory , 18 : 1 – 172 . OpenURL Placeholder Text WorldCat Kerr L. A. , Hintzen N. T., Cadrin S. X., Clausen L. W., Dickey-Collas M., Goethel D. R., Hatfield E. M. C. et al. 2017 . Lessons learned from practical approaches to reconcile mismatches between biological population structure and stock units of marine fish . ICES Journal of Marine Science , 74 : 1708 – 1722 . Google Scholar Crossref Search ADS WorldCat Lloris D. , Matallanas J., Oliver P. 2005 . Hakes of the World (Family Merlucciidae). An Annotated and Illustrated Catalogue of Hake Species Know to Date. FAO , Rome . 57 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lombarte A. , Lleonart J. 1993 . Otolith size changes related with body growth, habitat depth and temperature . Environmental Biology of Fishes , 37 : 297 – 306 . Google Scholar Crossref Search ADS WorldCat Manly B. F. J. 2007 . Randomization, Bootstrap and Monte Carlo Methods in Biology . Chapman & Hall , Boca Raton . 447 pp . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Minto C. , Myers R. A., Blanchard W. 2008 . Survival variability and population density in fish populations . Nature , 452 : 344 – 348 . Google Scholar Crossref Search ADS PubMed WorldCat Muto E. Y. , Soares L. S. H. 2011 . Spatio-temporal variations in the diet and stable isotope composition of the Argentine hake Merluccius hubbsi Marini, 1933 of the continental shelf of Southeastern Brazil . Marine Biology , 158 : 1619 – 1630 . Google Scholar Crossref Search ADS WorldCat Panikian G. , Cussens J., Pitchford W. 2015 . Identification and quantification of heteroscedasticity in stock-recruitment relationships . Canadian Journal of Fisheries and Aquatic Sciences , 72 : 1259 – 1271 . Google Scholar Crossref Search ADS WorldCat Perez J. A. , Pezzuto P. R., Wahrlich R., Soares A. L. S. 2009 . Deep-water fisheries in Brazil: history, status and perspectives . Latin American Journal of Aquatic Research , 37 : 513 – 542 . Google Scholar Crossref Search ADS WorldCat Piola A. R. , Palma E. D., Bianchi A. A., Castro B. M., Dottori M., Guerrero R. A., Marrari M. et al. 2018 . Physical oceanography of the SW Atlantic Shelf: a review. In Plankton Ecology of the Southwestern Atlantic , pp. 37 – 56 . Ed. by Hoffmeyer M. S., Sabatini M. E., Brandini F. P., Calliari D. L., Santinelli N. H., Springer , Cham . 586 pp. Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Saborido-Rey F. , Kjesbu O. 2005 . Growth and maturation dynamics . Digital CSIC Paper , 1 – 26 . http://hdl.handle.net/10261/47150 (last accessed 4 December 2019). OpenURL Placeholder Text WorldCat Scheffé H. 1958 . Experiments with mixtures . Journal of the Royal Statistical Society: Series B , 20 : 344 – 360 . OpenURL Placeholder Text WorldCat Scheffé H. 1963 . The simplex-centroid design for experiments with mixtures . Journal of the Royal Statistical Society: Series B , 25 : 235 – 263 . OpenURL Placeholder Text WorldCat Seber G. A. F. 1984 . Multivariate Observations . Wiley , New York . 686 pp. Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Secor D. H. 2014 . The unit stock concept: bounded fish and fisheries. In Stock Identification Methods: Applications in Fishery Science , 2nd edn, pp. 7 – 28 . Ed. by Cadrin S. X., Kerr L. A., Mariani S., Elsevier , Oxford . 588 pp. Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Semmar N. 2011 . Computational Metabolomics . Nova Science Publishers , New York . 238 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Semmar N. 2013 . Native Statistics for Natural Sciences . Nova Science Publishers , New York . 515 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Semmar N. , Roux M. 2014 . A new simplex approach to highlight multi-scale feeding behaviors in forager species from stomach contents: application to insectivore lizard population . Biosystems , 118 : 60 – 75 . Google Scholar Crossref Search ADS PubMed WorldCat Sokal R. R. , Rohlf F. J. 1995 . Biometry: The Principles and Practice of Statistics in Biological Research , 3rd edn. W.H. Freeman and Company , New York . 887 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Vaz-dos-Santos A. M. , Rossi-Wongtschowski C. L. D. B. 2007 . Age and growth of the Argentine hake Merluccius hubbsi Marini, 1933 in the Brazilian South-Southeast Region during 1996-2001 . Neotropical Ichthyology , 5 : 375 – 386 . Google Scholar Crossref Search ADS WorldCat Vaz-dos-Santos A. M. , Rossi-Wongtschowski C. L. D. B., Figueiredo J. L. 2009 . Merluccius hubbsi (Teleostei: Merlucciidae): stock identification based on reproductive biology in the South-Southeast Brazilian region . Brazilian Journal of Oceanography , 57 : 17 – 31 . Google Scholar Crossref Search ADS WorldCat Vaz-dos-Santos A. M. , Rossi-Wongtschowski C. L. D. B., Figueiredo J. L., Ávila-da-Silva A. O. 2010 . Threatened fishes of the world: Merluccius hubbsi Marini, 1933 (Merlucciidae) . Environmental Biology of Fishes , 87 : 349 – 350 . Google Scholar Crossref Search ADS WorldCat Vaz-dos-Santos A. M. , Schwingel P. R. 2015 . Biology and fisheries of hake (Merluccius hubbsi) in Brazilian waters, Southwest Atlantic Ocean. In Hakes: Biology and Exploitation , pp. 211 – 233 . Ed. by Arancibia H.. John Wiley and Sons, Chichester, West Sussex . 348 pp. Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Vaz-dos-Santos A. M. , Santos-Cruz N. N., Souza D., Silva A. G., Gris B., Rossi-Wongtschowski C. L. D. B. 2017 . Otoliths sagittae of Merluccius hubbsi: an efficient tool for the differentiation of stocks in the Southwestern Atlantic . Brazilian Journal of Oceanography , 65 : 520 – 525 . Google Scholar Crossref Search ADS WorldCat Weatherley A. H. 1990 . Approaches to understanding fish growth . Transactions of the American Fisheries Society , 119 : 662 – 672 . Google Scholar Crossref Search ADS WorldCat Weatherley A. H. , Gill H. S. 1987 . The Biology of Fish Growth . Academic Press , London . 443 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC © International Council for the Exploration of the Sea 2019. All rights reserved. For permissions, please email: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
Non-linearity in stock–recruitment relationships of Atlantic cod: insights from a multi-model approachSguotti, Camilla; Otto, Saskia A; Cormon, Xochitl; Werner, Karl M; Deyle, Ethan; Sugihara, George; Möllmann, Christian
doi: 10.1093/icesjms/fsz113pmid: N/A
Abstract The stock–recruitment relationship is the basis of any stock prediction and thus fundamental for fishery management. Traditional parametric stock–recruitment models often poorly fit empirical data, nevertheless they are still the rule in fish stock assessment procedures. We here apply a multi-model approach to predict recruitment of 20 Atlantic cod (Gadus morhua) stocks as a function of adult biomass and environmental variables. We compare the traditional Ricker model with two non-parametric approaches: (i) the stochastic cusp model from catastrophe theory and (ii) multivariate simplex projections, based on attractor state-space reconstruction. We show that the performance of each model is contingent on the historical dynamics of individual stocks, and that stocks which experienced abrupt and state-dependent dynamics are best modelled using non-parametric approaches. These dynamics are pervasive in Western stocks highlighting a geographical distinction between cod stocks, which have implications for their recovery potential. Furthermore, the addition of environmental variables always improved the models’ predictive power indicating that they should be considered in stock assessment and management routines. Using our multi-model approach, we demonstrate that we should be more flexible when modelling recruitment and tailor our approaches to the dynamical properties of each individual stock. Introduction Forecasting complex trajectories of marine resources is essential to fishery management and one of the major challenges of our time (Schindler and Hilborn, 2015; Ye et al., 2015). An important factor to be considered in fishery management is the stock–recruitment relationship (SRR), which serves as a basis for any stock assessment procedure to ultimately calculate reference points (Hilborn, 2002; ICES, 2017, 2018). SRRs are based on the assumption that recruitment (the number of fish that enter the adult population) is directly related to adult stock size (Kraus et al., 2000; Jennings et al., 2001). Parametric approaches, such as the Ricker model, were developed around the 1950s (Ricker, 1954) and in some cases still represent the method of choice in stock assessments (ICES, 2017). These models are very specific in the type of functional response curve to describe the SRR, and are linear, in the sense that, the relationship between recruitment and biomass can be linearized through log-transformation (Ye et al., 2015). However, they often fail to capture the high variability in recruitment data and this has led to questioning the existence of the relationship itself (Szuwalski et al., 2015; Britten et al., 2016; Perlala et al., 2017). The fit of the SRR is sometimes so poor, that short-term predictions of spawning-stock biomass (SSB) are conducted using an average of recruitment over a particular number of years, instead of a SRR model (Deyle et al., 2018). Both approaches, using average recruitment or a parametric model, assume that natural systems behave in a linear way, which may lead to biased fishery management decisions when stocks show complex dynamics such as aperiodic chaos, non-linearity, or non-stationarity (Ye and Sugihara, 2016; Perlala et al., 2017; Deyle et al., 2018). Chaos and non-stationary dynamics are pervasive in natural systems and characterize many marine ecosystems and populations (May and Oster, 1976; Scheffer et al., 2001; Möllmann et al., 2015). These dynamics emerge from the inherent complexity of nature, governed by a multitude of factors (Ye et al., 2015; Deyle et al., 2016; Tu et al., 2018). Assuming linearity and stability in recruitment models can, thus, result in wrong stock predictions (Glaser et al., 2014; Ye et al., 2015). As a consequence, new non-parametric modelling frameworks were developed to predict stock trajectories accounting for state-dependent and chaotic behaviour, such as the empirical dynamic modelling (EDM) framework (Sugihara et al., 2012; Ye et al., 2015; Deyle et al., 2018). EDM is a minimal assumptive approach based on time-series observations, which reconstructs the temporal dynamics of a system by constructing a so-called attractor manifold (Sugihara et al., 2012; Ye et al., 2015). EDM is able to predict the future system trajectory based on its past dynamics (Ye et al., 2015; Deyle et al., 2018), thus accounting for state-dependent dynamics (Sugihara, 1994). This approach, and in particular multivariate simplex projection (MSP) has been applied to predict non-linear fish recruitment dynamics in a range of studies, and has also been applied directly to management, e.g. for the menhaden stocks along the East Coast of the United States (Perretti et al., 2015; Ye et al., 2015; Deyle et al., 2018). Another non-parametric approach suitable for modelling state-dependent and discontinuous recruitment dynamics is the stochastic cusp model (SCM), which is based on catastrophe theory (Zeeman, 1976; Thom, 1977; Grasman et al., 2009; Petraitis and Dudgeon, 2016; Sguotti et al., 2019). Here, a state variable z (for instance recruitment), depends on two control variables alpha and beta. The model allows z to move from a state A (e.g. high recruitment) to a state B (e.g. low recruitment) following either a continuous or discontinuous path (Diks and Wang, 2016). SCM has been widely applied to economic and behavioural studies (van der Maas et al., 2003; Diks and Wang, 2016), but to a lesser degree to marine ecological studies (Jones and Walters, 1976; Jones, 1977; Petraitis and Dudgeon, 2015; Sguotti et al., 2019). Another point often neglected in recruitment prediction is the effect of multiple external drivers and potential interactions such as predation, competition, and environmental variables (Myers et al., 1995; Brander, 2005; Ottersen et al., 2006; Stiasny et al., 2016). However, in multiple cases the relationship between recruitment and environment can be spurious, non-linear, or non-stationary, and therefore is often not considered in stock assessments (Myers, 1998; Perlala et al., 2017). Traditional stock–recruitment models, which often are parametric models, assuming fixed parameters, usually fail to correctly incorporate the environmental information, because they often just consider additive effects of SSB and climate variables. Instead, non-parametric models such as MSP and SCM, can model interactions between the different drivers (i.e. biomass and climate variables) and thus may be able to integrate the climate information correctly (Ye et al., 2015; Deyle et al., 2018; Sguotti et al., 2019). This is important because for effectively predicting the status of living marine resources the integration of environmental variables is becoming crucial given the widespread impacts of climate change on ecosystems and marine resources such as commercially important fish (Britten et al., 2016; Gaines et al., 2018). Atlantic cod (Gadus morhua) is an iconic species from ecological, cultural, and economic points of view (Myers et al., 1996). In recent decades, most North Atlantic cod stocks have collapsed, followed by prolonged periods of no recovery even after the application of strict management measures (e.g. fishing moratoria; Myers et al., 1996; Hutchings, 2000; Hutchings and Rangeley, 2011; Frank et al., 2016; Sguotti et al., 2019). This failed recovery of Atlantic cod stocks suggests the presence of discontinuous dynamics and hysteresis (Frank et al., 2011; Steneck et al., 2011; Sguotti et al., 2019). Eastern and western Atlantic stocks differ in life-history traits, exploitation trajectories and recovery potential (Pörtner et al., 2008; Wang et al., 2014; Frank et al., 2016). Indeed, stocks in the West collapsed more abruptly compared to stocks in the East, which on average show more gradual declines (Frank et al., 2016). Cod recruitment is highly state-dependent, depending on the dimension of the stock and environment conditions. Recruitment is fundamental to Atlantic cod recovery (Myers and Barrowman, 1996; Brander, 2005) and influenced by climate change (Myers and Drinkwater, 1989; Planque and Frédou, 1999; Stige et al., 2006; Pörtner et al., 2008; Pershing et al., 2015). We here used stock assessment data from 20 Atlantic cod stocks to (i) investigate whether cod recruitment can be best described by the parametric Ricker model, by the non-parametric, “discontinuous” SCM, or by the non-parametric, state-dependent MSP approach, and (ii) test whether the model's predictive power can be improved when including environmental variables. We show that the adoption of a multi-model approach should be considered when modelling stocks presenting different dynamics. Material and methods Data We used recruitment (i.e. number of fish for a particular age and stock that recruit to the adult biomass in thousands, R) and SSB (i.e. biomass of mature fish in tonnes) data derived from stock assessments of 20 Atlantic cod stocks (Figure 1, Supplementary Figure S1). Data were provided by the International Council for the Exploration of the Sea (ICES), the National Oceanic and Atmospheric Administration of the United States (NOAA), the northwest Atlantic Fisheries Organization (NAFO), the Department of Fisheries and Ocean in Canada (DFO), and by personal communication (Supplementary Table S1). Recent assessments for cod stocks in the Kattegat, the western Baltic as well as the Norwegian coast have been conducted only for reduced periods. Therefore, we combined recent and older stock assessments after consistency checks of SSB and R time-series, by simply replacing the newer part of the time-series of the older assessments with the newer time-series assessment (see Supplementary Figure S2). Figure 1. Open in new tabDownload slide Map of cod stocks in the North Atlantic. Each circle corresponds to the centre of distribution of an Atlantic cod stock. The colour code corresponds to the division between western (orange, 12–20) and eastern stocks (pink, 1–11). Figure 1. Open in new tabDownload slide Map of cod stocks in the North Atlantic. Each circle corresponds to the centre of distribution of an Atlantic cod stock. The colour code corresponds to the division between western (orange, 12–20) and eastern stocks (pink, 1–11). We selected sea surface temperature (SST) and the indices of the North Atlantic oscillation (NAO) and Atlantic multidecadal oscillation (AMO) as climate predictors in our models. SST data were collected from the NOAA Extended Reconstructed Sea Surface Temperature dataset (ERSST, www.ncdc.noaa.gov) version 4. The dataset represents a reconstruction of SST from 1854 to the present and comprises monthly anomalies computed with respect to the period 1971–2000, resolved in a 2° × 2° grid of spatial resolution. The data were averaged per year and per management unit. SST was chosen because of its importance for recruitment of Atlantic cod and is also a proxy for climate change at a local scale (Planque and Frédou, 1999). NAO and AMO were used as indices of climate variability at the supraregional scale. In particular NAO has been shown to highly correlate with Atlantic cod recruitment (Stige et al., 2006), whereas AMO is a good proxy for climate change at longer timescales in this area. The NAO is a large-scale, high frequency (7–25 years) climatic index depending on the different atmospheric pressure at sea level between Iceland and Azores. The AMO is instead a large-scale, low frequency (60 years) multidecadal index representing climate-related SST changes in the Atlantic Ocean. The data for both indices were collected from the Earth System Research Laboratory of NOAA (www.esrl.noaa.gov), and the AMO was averaged to annual values, whereas the NAO was averaged annually but just between December and March. Modelling strategy We compared multiple stock–recruitment models, the traditional Ricker model, the SCM, and MSPs (from the EDM framework). Recruitment models include either SSB alone or SSB in combination with one of the climate variables (i.e. SST, NAO, and AMO) as predictors. Because recruitment can be influenced by climatic factors at different life-stages (i.e. eggs, larvae, and juveniles), we applied multiple lags on the climate variables depending on recruitment age (Supplementary Table S1). We assessed the predictive power of the different models (three modelling approaches and explanatory variables and corresponding lags selection) on the test data using fivefold cross-validation, which randomly splits the time-series in five parts, fitting the model on four (training data), and using the results to predict the last one (test data). In each of the five iterations, we compared the predicted with the observed test values using Pearson correlation coefficients (ρ; Ye et al., 2015; Deyle et al., 2018). We repeated this procedure 100 times to increase the robustness and eventually used the median of the 500 values of ρ for model comparison (Figure 2). Figure 2. Open in new tabDownload slide Schematic summary of the modelling approach. Box 1 shows the three model types applied. In the Ricker model (left) a curve depending on parameters (a) and (b) is fitted to show the relationship between recruitment (R) and SSB. A third parameter (c), allows the introduction of a climate variable and thus allows the curve to change between different climate states (i.e. above the mean in red, below the mean in blue). The stochastic cusp model (SCM; center) is shown both in a three-dimensional landscape and its projection in two-dimensional. In SCM the state variable R depends on two control variables (i.e. SSB and temperature). Although SSB controls the dimension of R, i.e. if R is found at the upper or lower shield, temperature controls the type of path that R will follow, either linear or discontinuous (s-haped, i.e. the two red paths). The three-dimensional landscape can be projected on the two-dimensional plane in which the folded area, i.e. the area of discontinuous dynamics, is shown in grey. In the attractor reconstruction of R depending on SSB and climate made with multivariate simplex projection (MSP; right) every point in the attractor correspond to a time-step of the system. MSP allows for the prediction of the future step of the system based on Euclidean Distance Calculations and thus is a state-dependent approach. All methods were fit as baseline models using just SSB as control variable then adding the environmental variable. Box 2 shows the model evaluation procedure using fivefold cross-validation. From this procedure, the predictive power of all the models was calculated and finally compared (Box 3) to assess performance among models of the same type and between the three different methods. Figure 2. Open in new tabDownload slide Schematic summary of the modelling approach. Box 1 shows the three model types applied. In the Ricker model (left) a curve depending on parameters (a) and (b) is fitted to show the relationship between recruitment (R) and SSB. A third parameter (c), allows the introduction of a climate variable and thus allows the curve to change between different climate states (i.e. above the mean in red, below the mean in blue). The stochastic cusp model (SCM; center) is shown both in a three-dimensional landscape and its projection in two-dimensional. In SCM the state variable R depends on two control variables (i.e. SSB and temperature). Although SSB controls the dimension of R, i.e. if R is found at the upper or lower shield, temperature controls the type of path that R will follow, either linear or discontinuous (s-haped, i.e. the two red paths). The three-dimensional landscape can be projected on the two-dimensional plane in which the folded area, i.e. the area of discontinuous dynamics, is shown in grey. In the attractor reconstruction of R depending on SSB and climate made with multivariate simplex projection (MSP; right) every point in the attractor correspond to a time-step of the system. MSP allows for the prediction of the future step of the system based on Euclidean Distance Calculations and thus is a state-dependent approach. All methods were fit as baseline models using just SSB as control variable then adding the environmental variable. Box 2 shows the model evaluation procedure using fivefold cross-validation. From this procedure, the predictive power of all the models was calculated and finally compared (Box 3) to assess performance among models of the same type and between the three different methods. The recruitment models The Ricker Model fits a curve between recruitment and SSB depending on parameters a and b (Ricker, 1954). These parameters allow the curve to bend in the middle, so that at very high SSB values recruitment will be low because of density-dependent effects. However, this model is log-linear, i.e. the relationship between recruitment and biomass can be linearized through log-transformation, thus, we will refer to it as a linear model throughout the text. Climate effects can be added through a new parameter (c; Figure 2, box1): Rt= SSBtexp (a-b*SSBt-ageR)(1a) Rt= SSBtexp a-b*SSBt-ageR+c*climatet-lags,(1b) where ageR is the age at recruitment, and lags the offset between the effect of a climate variable and R depending on the age of recruitment [i.e. for each stocks the climate variables were lagged from Rt to Rt-ageR depending on the age at recruitment (Supplementary Table S1)]. The starting values for the parameters were estimated from the linearized version of the function using the FSA (Ogle, 2016) package [Equation (2)]: log RtSSBt=SSBtexp a-b*SSBt-ageR+c*climatet-lags.(2) Subsequently the Stock–Recruitment function was fitted to the data using a non-linear model with as response variable the log-transformed recruitment. SCM is based on the cusp, one of the seven canonical forms of catastrophe theory that describe sudden changes in a system because of small changes of external drivers (Thom, 1977; van der Maas et al., 2003; Petraitis and Dudgeon, 2016; Sguotti et al., 2019). The cusp model is based on a cubic differential Equation (3) and describes discontinuous transitions in a state variable zt, controlled by two control variables α and β, and thus can be used to describe discontinuous dynamics in recruitment (Figure 2, box 1). -Vzt;α,β=-14zt4+12βzt2+αzt,(3) where Vzt;α,β is a potential function representing the rate of change of the system (zt), depending on the two control variables ( α,β). Because natural processes and empirical data often include stochasticity, Equation (3) was reformulated as a stochastic differential equation, adding the Wiener process ( σzdWt ) with variance σ2: -∂Vz;α,β∂z=-zt3+βzt+αdt + σzdWt=0,(4) where the first part of the equation is the drift term, σz is the diffusion parameter, and Wt represents the Wiener process. The state variable, zt, and the parameters, α and β [Equations (3) and (4)] are estimated as a linear function of one or more exogenous variables using a likelihood approach [Equations (5a)–(5c)]. z= w0+w1y1+w2y2+…wyyy(5a) α= α0+α1x1+α2x2+…αyxy(5b) β= β0+β1x1+β2x2+…βyxy(5c) where α0,β0,and z0 are the intercepts and α1,β1,and z1 the slopes of the models. In our study zt , the state variable, is a linear function of recruitment. α is the so-called asymmetry parameter and controls the size of zt, thus in our study is a function of SSB. β is called the bifurcation parameter because it controls whether the state variable follows a continuous or discontinuous path (Petraitis and Dudgeon, 2015; Diks and Wang, 2016; Sguotti et al., 2019), and in our study is a function of the environmental variables [climate, see below Equations (7a)–(7c)]. The system presents multiple equilibria if it follows a discontinuous path (i.e. two stable and one unstable) and just one if it follows a continuous path. The number of equilibria of the system depends on the solution of Equation (4), from which the Cardan’s discriminant (δ) is derived: δ=27α2-4β3.(6) If δ > 0, the system has one equilibrium, indicating a continuous path. Whereas if δ ≤0 the system has three equilibria, indicating a discontinuous path (Diks and Wang, 2016). Therefore, SCM allows the detection of interactive effects of the two control variables on the state variable. Any changes in the bifurcation parameter β, will lead to changes in the relationship between α and zt and consequently dramatic changes of the state variable (Figure 2, box1, Supplementary Figure S3). The model is represented by: zt= w0+w1Rt(7a) α= α0+α1SSBt-ageR(7b) β= β0+β1SSBt−ageR or β= β0+β1climatet−lags(7c) In order to test the predictive power of the model, we first produced the linear predictors of the parameters and the state variable. These were then fit into Equation (8) to predict the new points on the surface. V´(z^)=α+βz^-z^3 =0.(8) MSP is based on the EDM framework. The cornerstone of this framework is the Simplex Projection method. The principle of EDM is to reconstruct the dynamics of one or multiple time-series in a multidimensional space, i.e. an attractor manifold, and predict the future trajectory of the system based on these past dynamics (Figure 2, box1; Sugihara et al., 2012; Ye et al., 2015; Chang et al., 2017). Reconstructing the past dynamics of a system (in our case recruitment) is possible either using multiple variables (i.e. SSB or climate indices) or just time-lags of the system itself (i.e. recruitment; Sugihara et al., 2012). We here used differentiated recruitment time-series to build the attractor for each cod stock, and Simplex Projection [Equations (9) and (10)] to approximate the attractor dynamics of the system (Sugihara et al., 2012; Ye et al., 2015; Deyle et al., 2018). The time-series is transformed in a set of time-delayed coordinate vectors: x={xt, xt−τ, xt−2τ, xt−3τ,…xt−(E−1)τ},(9) where x is the system, in our study recruitment, t is time, τ is the time-lag, and E the embedding dimension. E represents the dimensionality of the attractor (Ye et al., 2015). E is selected by predicting the attractor manifold one step ahead into the future (using leave-one-out cross-validation) then comparing the predictive power of models with a varying E. In order to predict the system into the future, x^t+1 , Euclidean distance is used and the system is predicted using nearest neighbourhood estimations x^t+1= (∑i=1E+1wi,txi,t+1)∑i=1E+1wi,t,(10) where wi represents the weights [Equation (11)], which are the Euclidean distance to the neighbour vector i relative to the nearest neighbour d¯ . wi=exp (-dxt,xid¯).(11) MSP uses Equation (9) but with multiple variables. In our study, the attractor reconstruction of recruitment was based on SSB alone or together with climate variables [climate; Equations (12a) and (12b)]: Rt=SSBt-ageR(12a) Rt=SSBt-ageR, climatet-lags(12b) Performing MSP requires two preliminary tests, the S-Map and the convergent cross mapping (CCM), to unravel recruitment dynamics and the relationship between recruitment and explanatory variables, respectively. EDM-specific preliminary tests S-Map and CCM The S-Map, was performed after the attractor reconstruction with Simplex Projection. This test includes a tuning parameter θ that controls the weights wifrom Equation (10), and, if bigger than 0 indicates non-linearity (Sugihara, 1994; Klein et al., 2016; Dakos et al., 2017). Significance of non-linearity was assessed using a null distribution generated from 500 surrogate time-series for each S-Map model. The surrogate time-series were created following Deyle et al. (2018) and were phase-randomized, which preserves the basic statistical properties of the original time-series (Ebisuzaki, 1997). We averaged the S-Map results for all Atlantic cod recruitment time-series to understand the overall dynamics. We performed CCM between R and SSB and the climate variables (SST, NAO, and AMO), a technique to understand causality between time-series without assuming any distribution (Sugihara et al., 2012; Deyle et al., 2016; Pierre et al., 2018). CCM is based on the principle that, if SSB or climate variables have an influence, then the R time-series will contain information about the past state of these variables. CCM is performed using Equation (10; see Deyle et al., 2018). Software All analyses were performed in the programming environment R (R Core Team, 2016, version 3.3.1) using the packages FSA (Ogle, 2016), cusp (Grasman et al., 2009), rEDM (Ye et al., 2016). Results In our multi-model approach, we compared the parametric, linear Ricker model with two non-parametric, state-dependent approaches, i.e. the SCM and the state-dependent MSP, with or without environmental variables as predictors (Figure 2). The two preliminary tests of the EDM, necessary to perform the MSP, revealed on average significantly non-linear dynamics in recruitment of Atlantic cod stocks, and an appropriate choice of explanatory variables (Supplementary Figures S4 and S5), thus allowing us to proceed with the analyses. For most of the Atlantic cod stocks, the best performing models produced high correlations between observed and predicted values (0.7 < ρ < 0.8). An exception were northeast Arctic, Iceland, and Gulf of Maine cod stocks where the predictive power was lower compared to the other stocks (about ρ = 0.4). Differences between the three model types were in general low (Figure 3). The Ricker model performed best for seven stocks, the SCM for eight stocks, and the MSP for five stocks (Figure 3, Supplementary Table S2). For stocks where SCM was the best, the MSP generally showed also a high predictive power, indicating that both models can well describe abrupt dynamics [e.g. Figure 3 (12, 13, 14, 17)]. The addition of climate variables as explanatory variables to the baseline SSB models generally increased the predictive power, independently of the model type, although SSB was often the most correlated explanatory variable (Figure 3, Supplementary Table S2 and as shown in CCM, Supplementary Figure S5). SST and AMO were selected, based on the predictive power of the model, in respectively eight stocks and NAO in the remaining four stocks, generally agreeing with CCM results (Figure 4, Supplementary Figure S5). However, adding a climate variable had only a weak or even no additional effect when the baseline SSB model performed already poorly [e.g. Figure 3 (8), Ricker model]. Figure 3. Open in new tabDownload slide Stock–recruitment model comparison. The comparison between the predictive power of the best models resulting from the model selection between the Ricker, stochastic cusp model (SCM), and multivariate simplex projection (MSP) models. The median of the predictive power, derived from the cross-validation is shown for the three models without (blue, darker) and with (green, lighter) the inclusion of climate variables. The best model among the three, i.e. the model presenting the highest Pearson ρ between observed and predicted values of the test dataset, is indicated by a yellow (lighter) star for each stock. Black stars indicate the best models which however had a poor fit to recruitment and thus were substituted by the second-best model. The environmental variable that resulted in the best predictions can be found in Figure 4 and Supplementary Table S2. The number of years in the time-series is indicated for each stock. The colours underlying the names of the stocks correspond to the geographical location of the stock in pink (or from 1 to 11) in the East Atlantic, and orange (from 12 to 20) in the West. The numbers correspond to the stocks number in Figure1. Figure 3. Open in new tabDownload slide Stock–recruitment model comparison. The comparison between the predictive power of the best models resulting from the model selection between the Ricker, stochastic cusp model (SCM), and multivariate simplex projection (MSP) models. The median of the predictive power, derived from the cross-validation is shown for the three models without (blue, darker) and with (green, lighter) the inclusion of climate variables. The best model among the three, i.e. the model presenting the highest Pearson ρ between observed and predicted values of the test dataset, is indicated by a yellow (lighter) star for each stock. Black stars indicate the best models which however had a poor fit to recruitment and thus were substituted by the second-best model. The environmental variable that resulted in the best predictions can be found in Figure 4 and Supplementary Table S2. The number of years in the time-series is indicated for each stock. The colours underlying the names of the stocks correspond to the geographical location of the stock in pink (or from 1 to 11) in the East Atlantic, and orange (from 12 to 20) in the West. The numbers correspond to the stocks number in Figure1. Figure 4. Open in new tabDownload slide Visualization of stock–recruitment relationships (SRR) for Atlantic cod stocks. Every panel contains a representation of the best model for each of the stocks. Abbreviations used: spawning-stock biomass (SSB), recruitment (R), sea surface temperature (SST), North Atlantic oscillation (NAO), and Atlantic multidecadal oscillation (AMO). (4, 5, 6, 10, 9, 11, 19). Results for cod stocks best represented by the Ricker model. The colour of the dots corresponds to the state of the climate variable, red (lighter) above the mean, blue (darker) below the mean. The two lines indicate separate SRRs for the two climate states (above and below the mean. (12–18, 1). Results for cod stocks best represented by the stochastic cusp model (SCM), dots are scaled to the size of R. The blue area (or grey in black and white) corresponds to the instability area, thus the fold in the three-dimensional visualization (Figure 2, Supplementary Figure S3) where three equilibria are possible. (3, 2, 7, 8, 20). Results for cod stocks best represented by multivariate simplex projection (MSP); dots are scaled to the size of SSB. Colours correspond to the state of the climate variable indicated, red (lighter) above the mean, blue below the mean. The lines show the predicted trends of R over time. The colours in the boxes correspond to the geographical location of the stock in pink (or from 1 to 11) in the East Atlantic, and orange (from 12 to 20) in the West. The numbers correspond to stock numbers in Figure 1. The comparison between observed and predicted values for each model can be seen in Supplementary Figure S6. Figure 4. Open in new tabDownload slide Visualization of stock–recruitment relationships (SRR) for Atlantic cod stocks. Every panel contains a representation of the best model for each of the stocks. Abbreviations used: spawning-stock biomass (SSB), recruitment (R), sea surface temperature (SST), North Atlantic oscillation (NAO), and Atlantic multidecadal oscillation (AMO). (4, 5, 6, 10, 9, 11, 19). Results for cod stocks best represented by the Ricker model. The colour of the dots corresponds to the state of the climate variable, red (lighter) above the mean, blue (darker) below the mean. The two lines indicate separate SRRs for the two climate states (above and below the mean. (12–18, 1). Results for cod stocks best represented by the stochastic cusp model (SCM), dots are scaled to the size of R. The blue area (or grey in black and white) corresponds to the instability area, thus the fold in the three-dimensional visualization (Figure 2, Supplementary Figure S3) where three equilibria are possible. (3, 2, 7, 8, 20). Results for cod stocks best represented by multivariate simplex projection (MSP); dots are scaled to the size of SSB. Colours correspond to the state of the climate variable indicated, red (lighter) above the mean, blue below the mean. The lines show the predicted trends of R over time. The colours in the boxes correspond to the geographical location of the stock in pink (or from 1 to 11) in the East Atlantic, and orange (from 12 to 20) in the West. The numbers correspond to stock numbers in Figure 1. The comparison between observed and predicted values for each model can be seen in Supplementary Figure S6. The Ricker model best represented more gradual declines in recruitment, typical for cod stocks around the British Isles (i.e. North Sea, West of Scotland, and Irish Sea), those closer to the Arctic (i.e. Faroe Plateau, northeast Arctic, and Iceland cod) and Georges Bank cod [Figure 4 (4, 5, 6, 10, 9, 11, 19), Supplementary Table S3], as illustrated by their individual time-series (Supplementary Figure S1). All of these stocks, except Georges Bank, displayed strong density-dependence in recruitment (Supplementary Table S3), which is characteristic for the Ricker model. Furthermore, Ricker models clearly revealed that recruitment in warmer years is usually lower for the same level of SSB when compared to colder conditions (as indicated by low SST, NAO, or AMO in Figure 4 (4, 5, 6). The only exception with the reverse pattern of higher recruitment values at warmer conditions was northeast Arctic, hence the only cod stock that really profited from climate warming [Figure 4 (9)]. SCM instead is an approach from catastrophe theory, which models best discontinuous dynamics characterized by abrupt sudden shifts and hysteresis (i.e. in this case delayed recovery). The recruitment and SSB time-series of Canadian stocks on the western Atlantic side, but also Greenland and eastern Baltic cod (Supplementary Figure S1) show this type of dynamics, and hence SCM was the best approach for these stocks. SCM identified discontinuous stock–recruitment dynamics caused by the interaction of SSB and the climate variable. Moreover, SCM can identify catastrophic collapse, which occurs when SSB is found in the “folded” area, or area of instability (see blue shaded areas in Figure 4 (12, 13, 14, 15, 16, 17, 18, 1), Supplementary Figure S3). Recruitment collapsed in these stocks, when in the instability area, in response to only small reductions in stock size [Figure 4 (12–18, 1)]. Consequently, SSB was a significant predictor in all SCMs, controlling the stocks dimension, while the climate variables modified the relationship between recruitment and SSB rendering it discontinuous, and thus inducing hysteresis (Supplementary Table S4). These two factors lead to the presence of stable low recruitment levels towards the end of the time-series. Low SSB coupled with warming (as indicated by climate variables SST, NAO, and AMO, Supplementary Table S4) had the potential to stabilize low recruitment. This is indicated by values outside the bifurcation area as best demonstrated by northern and Grand Banks cod [Figure 4 (13, 17)]. Other cod stocks such as those from the Gulf of St Lawrence, on the eastern Scotian Shelf and off Greenland were at the boarder of stable low recruitment levels [Figure 4 (12, 14–16)]. Eventually, we found MSP to be the best model for recruitment of stocks that did not show collapses, but mostly fluctuating dynamics such as cod in the western Baltic, the Kattegat (because, even if the SCM was the best the model, the fit was invalid), the Celtic Sea, the Norwegian coast and in the Gulf of Maine [Figure 4 (3, 2, 7, 8, 20), Supplementary Figure S1]. The MSP however, seemed also appropriate to model catastrophic dynamics, but less effectively than the SCM. In contrast to the stocks best modelled with SCM and Ricker, stocks best modelled with MSP showed a mixed response to recent warming with a clear negative effect on recruitment in the western Baltic only [Figure 3 (2)]. Discussion Short-term predictions of the size of an incoming year-class is essential to modern assessments of commercial fish species, but often suffers from the accuracy of available models predicting recruitment based on continuous, linear relationship with SSB. In our study, we investigated (i) whether recruitment dynamics in Atlantic cod stocks are better predicted by non-parametric, state-dependent, or catastrophic statistical methodology compared to traditional parametric, linear approaches such as the Ricker stock–recruitment model and (ii) whether using climate variables as predictors in addition to SSB improves the predictive performance of such models. The main result of our study is that predicting fish stock–recruitment can be improved by tailoring the modelling approach to the dynamical properties of each individual stock. We found cod stocks with more gradual and mostly linear dynamics to be best predicted by the traditional linear Ricker model, whereas stocks that experienced sudden abrupt changes in recruitment and stock size are best described by the SCM. SCM, based on catastrophe theory, is well suited to represent such discontinuous regime shift dynamics (Thom, 1972; Grasman et al., 2009; Diks and Wang, 2016; Sguotti et al., 2019). SCM allows for the identification of drivers and how their interactions result in unstable recruitment dynamics and hence provides a form of vulnerability assessment that can be instrumental in management (Petraitis and Dudgeon, 2015; Diks and Wang, 2016; Sguotti et al., 2019). Eventually, MSP was most appropriate for stocks that displayed more chaotic and fluctuating behaviours (Sugihara et al., 2012; Ye et al., 2015). Indeed, being a minimally assumptive model the most complex dynamics are better captured by it. MSP as part of the EDM suite of methods is based on attractor reconstruction and accounts for state-dependent dynamics (Ye et al., 2015), which makes it a suitable approach to model also discontinuous dynamics (Ye et al., 2015; Deyle et al., 2018). Mostly, both SCM and MSP models performed similarly in our analysis and their relatively high predictive power indicated the importance of using state-dependent and/or discontinuous approaches to model recruitment (Ye et al., 2015; Deyle et al., 2018; Munch et al., 2018). Our study highlights that important differences exist between cod stocks in the eastern and western areas of the North Atlantic (Frank et al., 2016). Stocks from the western Atlantic and in particular off Canada and Greenland often experienced pronounced catastrophic dynamics, i.e. abrupt and sudden changes in stock size and recruitment. Eastern Atlantic stocks instead showed more continuous dynamics and thus a higher degree of stability. In general western Atlantic cod stocks seemed to be less resilient to abrupt collapses because of more fragile life-history traits, an overall more extreme and difficult environment, and different exploitation histories (Rätz and Lloret, 2003; Pörtner et al., 2008; Wang et al., 2014; Frank et al., 2016). Moreover, SST was selected in eastern Atlantic cod stocks models, whereas for Western stocks the climate indices explained better the recruitment variability. This difference might indicate that the eastern cod stocks are more influenced by local processes, whereas in the western Atlantic large-scale climatic fluctuations are more important. Nevertheless, the addition of the climate factors in the best stock–recruitment models almost always increased its predictive power and thus highlights the importance of using environmental information also in stock assessment and management considerations to consider broader ecosystem dynamics (Punt et al., 2013; Skern-Mauritzen et al., 2016). These results highlight the presence of multiple dynamics in cod stocks, which are also supported by the results of the preliminary S-Maps tests revealing a significant level of non-linearity in recruitment time-series of Atlantic cod stock. However, the non-linearity signal is lower than expected, which we assume is because of the nature of the stock assessment data we used, and thus could be an underestimation (Brooks and Deroba, 2015). Such model output tends to be smoother and more linear than survey data (Storch et al., 2017), which are unfortunately not available for all cod stocks and longer periods needed for our study. Finally, the different models allow us to draw conclusions about the recovery potential for collapsed Atlantic cod stocks. Most of the stocks are negatively influenced by warming and climate variability, because the lowest recruitment and SSB coincide with the highest temperature (Brander, 2005; Drinkwater, 2005; Pörtner et al., 2008). The only exception is northeast Arctic cod where a warming environment positively influences recruitment, because this stocks resides at the northern distribution limit of the species (Pörtner et al., 2008). Apart from northeast Arctic and Iceland cod where SSB has recently reached high levels, the stocks for which the traditional Ricker model performed best, such as the ones from the North Sea and around the British Isles, show continuously low recruitment and SSB in recent years and a continuous relationship between these parameters. These imply that, with low exploitation pressure these stocks have a higher recovery potential, but with climate change the productivity will likely remain low (Drinkwater, 2005). The situation is even worse for stocks that are best described by the SCM such as the western Atlantic stocks where the relationship between recruitment and SSB is discontinuous and thus the stocks display a strong hysteresis effect. Most of them are at present in a stable low state, suggesting that recovery might be even further delayed and productivity will remain low. Conclusions We demonstrated that discontinuous, state-dependent dynamics are pervasive in at least half of Atlantic cod stocks and need to be considered when predicting year-class strength. Indeed, even if our study does not necessarily reflect the goodness of the models to predict future recruitment, because the cross-validation included years after those predicted, we highlight the presence of different dynamics between stocks. Furthermore, we show the importance of accounting for environmental factors in recruitment predictions. Our findings indicate the need for more flexibility in the stock assessment process and highlight the importance for an adaptive multi-model approach that accounts for the inherent dynamics of living marine resource populations (Punt et al., 2016). Flexible models and adaptive management are fundamental to move towards an ecosystem-based management approach, especially in the face of climate change. To achieve this, we need to move away from fixed and established model procedures and explore other options, to be ready to adapt to the new challenges that climate change will pose (King et al., 2015). Supplementary data Supplementary material is available at the ICESJMS online version of the manuscript. Acknowledgements CS was supported through MARmaED (MARine MAnagement and Ecosystem Dynamics under climate change). The MARmaED project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 675997. The results of this publication reflect only the author's view and the commission is not responsible for any use that may be made of the information it contains. Further financial support was received (to CM, SAO, and XC) by the Federal Ministry of Education and Research of Germany in the framework of marEEshift (project No. 01LC17058). K-MW was supported by the CLIMA project, reference RER 15/0008, Ministry of Foreign Affairs Norway. GS was supported by the Department of Defence SERDP 15 RC-2509, National Science Foundation NSF-DEB-1020372, and NSF-ABI-Innovation DBI-1667584, McQuown Chair, UCSD. References Brander K. M. 2005 . Cod recruitment is strongly affected by climate when stock biomass is low . ICES Journal of Marine Science , 62 : 339 – 343 . Google Scholar Crossref Search ADS WorldCat Britten G. L. , Dowd M., Worm B. 2016 . Changing recruitment capacity in global fish stocks . Proceedings of the National Academy of Sciences of the United States of America , 113 : 134 – 139 . Google Scholar Crossref Search ADS PubMed WorldCat Brooks E. N. , Deroba J. J. 2015 . When “data” are not data: the pitfalls of post hoc analyses that use stock assessment model output . Canadian Journal of Fisheries and Aquatic Sciences , 72 : 634 – 641 . Google Scholar Crossref Search ADS WorldCat Chang C-W. , Ushio M., Hsieh C. 2017 . Empirical dynamic modeling for beginners . Ecological Research , 32 : 785 . Google Scholar Crossref Search ADS WorldCat Dakos V. , Glaser S. M., Hsieh C-H., Sugihara G. 2017 . Elevated nonlinearity as an indicator of shifts in the dynamics of populations under stress . Journal of Royal Society Interface , 14 : 20160845. Google Scholar Crossref Search ADS WorldCat Deyle E. , Schueller A. M., Ye H., Pao G. M., Sugihara G. 2018 . Ecosystem-based forecasts of recruitment in two menhaden species . Fish and Fisheries , 19 : 769 – 781 . Google Scholar Crossref Search ADS WorldCat Deyle E. R. , May R. M., Munch S. B., Sugihara G. 2016 . Tracking and forecasting ecosystem interactions in real time . Proceedings of the Royal Society B: Biological Sciences , 283 : 20152258. Google Scholar Crossref Search ADS WorldCat Diks C. , Wang J. 2016 . Can a stochastic cusp catastrophe model explain housing market crashes? Journal of Economic Dynamics and Control , 69 : 68 – 88 . Google Scholar Crossref Search ADS WorldCat Drinkwater K. F. 2005 . The response of Atlantic cod (Gadus morhua) to future climate change . ICES Journal of Marine Science , 62 : 1327 – 1337 . Google Scholar Crossref Search ADS WorldCat Ebisuzaki W. 1997 . A method to estimate the statistical significance of a correlation when the data are serially correlated . Journal of Climate , 10 : 2147 – 2153 . Google Scholar Crossref Search ADS WorldCat Frank K. T. , Petrie B., Fisher J. A. D., Leggett W. C. 2011 . Transient dynamics of an altered large marine ecosystem . Nature , 477 : 86 – 89 . Google Scholar Crossref Search ADS PubMed WorldCat Frank K. T. , Petrie B., Leggett W. C., Boyce D. G. 2016 . Large scale, synchronous variability of marine fish populations driven by commercial exploitation . Proceedings of the National Academy of Sciences , 113 : 8248 – 8253 . Google Scholar Crossref Search ADS WorldCat Gaines S. D. , Costello C., Owashi B., Mangin T., Bone J., García Molinos J., Burden M. 2018 . Improved fisheries management could offset many negative effects of climate change . Science Advances , 4 : eaao1378 . Google Scholar Crossref Search ADS PubMed WorldCat Glaser S. M. , Fogarty M. J., Liu H., Altman I., Hsieh C. H., Kaufman L., Maccall A. D., et al. 2014 . Complex dynamics may limit prediction in marine fisheries . Fish and Fisheries , 15 : 616 – 633 . Google Scholar Crossref Search ADS WorldCat Grasman R. P. P. P. , Van Der Maas H. L. J., Wagenmakers E-J. 2009 . Fitting the cusp catastrophe in R: a cusp package primer . Journal of Statistical Software , 32 : 1 – 27 . Google Scholar Crossref Search ADS WorldCat Hilborn R. 2002 . The dark side of reference points . Bulletin of Marine Science , 70 : 403 – 408 . Google Scholar OpenURL Placeholder Text WorldCat Hutchings J. A. 2000 . Collapse and recovery of marine fishes . Nature , 406 : 882 – 885 . Google Scholar Crossref Search ADS PubMed WorldCat Hutchings J. A. , Rangeley R. W. 2011 . Correlates of recovery for Canadian Atlantic cod (Gadus morhua) . Canadian Journal of Zoology , 89 : 386 – 400 . Google Scholar Crossref Search ADS WorldCat ICES. 2017 . ICES Advice Technical Guidelines. http://ices.dk/sites/pub/PublicationReports/Advice/2017/2017/12.04.03.01_Reference_points_for_category_1_and_2.pdf (last accessed 6 March 2019). ICES. 2018 . Advice Basis. http://www.ices.dk/explore-us/Documents/ICES and EBM.pdf (last accessed 6 March 2019). Jennings S. , Kaiser M. J., Reynolds J. D. 2001 . Marine Fisheries Ecology . Blackwell Science , Oxford, England . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Jones D. D. 1977 . Catastrophe theory applied to ecological systems . Simulation , 29 : 1 – 15 . Google Scholar Crossref Search ADS WorldCat Jones D. D. , Walters C. J. 1976 . Catastrophe theory and fisheries regulation . Journal of the Fisheries Research Board of Canada , 33 : 2829 – 2833 . Google Scholar Crossref Search ADS WorldCat King J. R. , Mcfarlane G. A., Punt A. E. 2015 . Shifts in fisheries management: adapting to regime shifts . Philosophical Transactions B , 370 : 20130277. Google Scholar Crossref Search ADS WorldCat Klein E. , Glaser S., Jordaan A., Kaufman L., Rosenberg A. 2016 . A complex past: historical and contemporary fisheries demonstrate nonlinear dynamics and a loss of determinism . Marine Ecology Progress Series , 557 : 237 – 246 . Google Scholar Crossref Search ADS WorldCat Kraus G. , Muller A., Trella K., Kouster F. W. 2000 . Fecundity of Baltic cod: temporal and spatial variation . Journal of Fish Biology , 56 : 1327 – 1341 . Google Scholar Crossref Search ADS WorldCat May R. M. , Oster G. F. 1976 . Bifurcations and dynamic complexity in simple ecological models . Source: The American Naturalist , 110 : 573 – 599 . Google Scholar Crossref Search ADS WorldCat Möllmann C. , Folke C., Edwards M., Conversi A., Mollmann C., Folke C., Edwards M. 2015 . Marine regime shifts around the globe: theory, drivers and impacts . Philosophical Transactions of the Royal Society B , 370 : 20130260. Google Scholar Crossref Search ADS WorldCat Munch S. B. , Giron-Nava A., Sugihara G. 2018 . Nonlinear dynamics and noise in fisheries recruitment: a global meta-analysis . Fish and Fisheries , 19 : 964 – 973 . Google Scholar Crossref Search ADS WorldCat Myers R. , Hutchings J., Barrowman N. 1996 . Hypotheses for the decline of cod in the North Atlantic . Marine Ecology Progress Series , 138 : 293 – 308 . Google Scholar Crossref Search ADS WorldCat Myers R. A. 1998 . When do environment-recruitment correlations work? Reviews in Fish Biology and Fisheries , 8 : 285 – 305 . Google Scholar Crossref Search ADS WorldCat Myers R. A. , Barrowman N. J. 1996 . Is fish recruitment related to spawner abundance? Fishery Bulletin , 94 : 707 – 724 . Google Scholar OpenURL Placeholder Text WorldCat Myers R. A. , Barrowman N. J., Hutchings J. A., Rosenberg A. A. 1995 . Population dynamics of exploited fish stocks at low population levels . Science , 269 : 1106 – 1108 . Google Scholar Crossref Search ADS PubMed WorldCat Myers R. A. , Drinkwater K. 1989 . The influence of Gulf Stream warm core rings on recruitment of fish in the northwest Atlantic . Journal of Marine Research , 47 : 635 – 656 . Google Scholar Crossref Search ADS WorldCat Ogle D. H. 2016 . FSA: Simple Fisheries Stock Assessment Methods. https://cran.r-project.org/web/packages/FSA/FSA.pdf (last accessed 4 September 2018). Ottersen G. , Hjermann D. O., Stenseth N. C. 2006 . Changes in spawning stock structure strengthen the link between climate and recruitment in a heavily fished cod (Gadus morhua) stock . Fisheries Oceanography , 15 : 230 – 243 . Google Scholar Crossref Search ADS WorldCat Perlala T. A. , Swain D. P., Kuparinen A. 2017 . Examining nonstationarity in the recruitment dynamics of fishes using Bayesian change point analysis . Canadian Journal of Fisheries and Aquatic Sciences , 74 : 751 – 765 . Google Scholar Crossref Search ADS WorldCat Perretti C. T. , Munch S. B., Fogarty M. J., Sugihara G. 2015 . Global evidence for non-random dynamics in fish recruitment. 1 – 22 . arXiv:1509.01434v1. Pershing A. J. , Alexander M. A., Hernandez C. M., Kerr L. A., Le Bris A., Mills K. E., Nye J. A., et al. 2015 . Slow adaptation in the face of rapid warming leads to collapse of the Gulf of Maine cod fishery . Science , 350 : 809 – 812 . Google Scholar Crossref Search ADS PubMed WorldCat Petraitis P. S. , Dudgeon S. R. 2015 . Variation in recruitment and the establishment of alternative community states . Ecology , 96 : 3186 – 3196 . Google Scholar Crossref Search ADS PubMed WorldCat Petraitis P. S. , Dudgeon S. R. 2016 . Cusps and butterflies: multiple stable states in marine systems as catastrophes . Marine and Freshwater Research , 67 : 37 – 46 . Google Scholar Crossref Search ADS WorldCat Pierre M. , Rouyer T., Bonhommeau S., Fromentin J. M. 2018 . Assessing causal links in fish stock–recruitment relationships . ICES Journal of Marine Science , 75 : 903 – 911 . Google Scholar Crossref Search ADS WorldCat Planque B. , Frédou T. 1999 . Temperature and the recruitment of Atlantic cod (Gadus morhua) . Canadian Journal of Fisheries and Aquatic Sciences , 56 : 2069 – 2077 . Google Scholar Crossref Search ADS WorldCat Pörtner H. O. , Bock C., Knust R., Lannig G., Lucassen M., Mark F. C., Sartoris F. J. 2008 . Cod and climate in a latitudinal cline: physiological analyses of climate effects in marine fishes . Climate Research , 37 : 253 – 270 . Google Scholar Crossref Search ADS WorldCat Punt A. E. , A'mar T., Bond N. A., Butterworth D. S., de Moor C. L., De Oliveira J. A. A., Haltuch M. A., et al. 2013 . Fisheries management under climate and environmental uncertainty: control rules and performance simulation . ICES Journal of Marine Science , 71 : 2208 – 2220 . Google Scholar Crossref Search ADS WorldCat Punt A. E. , Butterworth D. S., de Moor C. L., De Oliveira J. A. A., Haddon M. 2016 . Management strategy evaluation: best practices . Fish and Fisheries , 17 : 303 – 334 . Google Scholar Crossref Search ADS WorldCat Rätz H. J. , Lloret J. 2003 . Variation in fish condition between Atlantic cod (Gadus morhua) stocks, the effect on their productivity and management implications . Fisheries Research , 60 : 369 – 380 . Google Scholar Crossref Search ADS WorldCat R Core Team. 2016 . R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Ricker W. E. 1954 . Stock and recruitment . Journal of the Fisheries Research Board of Canada , 11 : 559 – 623 . Google Scholar Crossref Search ADS WorldCat Scheffer M. , Carpenter S. R., Foley J. A., Folke C., Walker B. 2001 . Catastrophic shifts in ecosystems . Nature , 413 : 591 – 596 . Google Scholar Crossref Search ADS PubMed WorldCat Schindler D. E. , Hilborn R. 2015 . Sustainability. Prediction, precaution, and policy under global change . Science (New York, N.Y.) , 347 : 953 – 954 . Google Scholar Crossref Search ADS PubMed WorldCat Sguotti C. , Otto S. A., Frelat R., Langbehn T. J., Ryberg M. P., Lindegren M., Durant J. M., et al. 2019 . Catastrophic dynamics limit Atlantic cod recovery . Proceedings of the Royal Society B , 286 : 20182877 . Google Scholar Crossref Search ADS PubMed WorldCat Skern-Mauritzen M. , Ottersen G., Handegard N. O., Huse G., Dingsør G. E., Stenseth N. C., Kjesbu O. S. 2016 . Ecosystem processes are rarely included in tactical fisheries management . Fish and Fisheries , 17 : 165 – 175 . Google Scholar Crossref Search ADS WorldCat Steneck R. S. , Hughes T. P., Cinner J. E., Adger W. N., Arnold S. N., Berkes F., Boudreau S. A., et al. 2011 . Creation of a gilded trap by the high economic value of the Maine lobster fishery . Conservation Biology , 25 : 904 – 912 . Google Scholar Crossref Search ADS PubMed WorldCat Stiasny M. H. , Mittermayer F. H., Sswat M., Voss R., Jutfelt F., Chierici M., Puvanendran V., et al. 2016 . Ocean acidification effects on Atlantic cod larval survival and recruitment to the fished population . PLoS One , 11 : 1 – 11 . Google Scholar Crossref Search ADS WorldCat Stige L. C. , Ottersen G., Brander K., Chan K. S., Stenseth N. C. 2006 . Cod and climate: effect of the North Atlantic oscillation on recruitment in the North Atlantic . Marine Ecology Progress Series , 325 : 227 – 241 . Google Scholar Crossref Search ADS WorldCat Storch L. S. , Glaser S. M., Ye H., Rosenberg A. A. 2017 . Stock assessment and end-to-end ecosystem models alter dynamics of fisheries data . PLoS One , 12 : e0171644. Google Scholar Crossref Search ADS PubMed WorldCat Sugihara G. 1994 . Nonlinear forecasting for the classification of natural time series. Philosophical transactions of the Royal Society of London . Series A Biological Sciences , 348 : 477 – 495 . Google Scholar OpenURL Placeholder Text WorldCat Sugihara G. , May R., Ye H., Hsieh C., Deyle E., Fogarty M., Munch S. 2012 . Detecting causality in complex ecosystems . Science , 338 : 496 – 500 . Google Scholar Crossref Search ADS PubMed WorldCat Szuwalski C. S. , Vert-Pre K. A., Punt A. E., Branch T. A., Hilborn R. 2015 . Examining common assumptions about recruitment: a meta-analysis of recruitment dynamics for worldwide marine fisheries . Fish and Fisheries , 16 : 633 – 648 . Google Scholar Crossref Search ADS WorldCat Thom R. 1972 . Structural Stability and Morphogenesis: An Outline of a Theory of Models . Benjamin , New York . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Thom R. 1977 . Structural stability, catastrophe theory and applied mathematics . SIAM Review , 19 : 189 – 201 . Google Scholar Crossref Search ADS WorldCat Tu C-Y. , Chen K-T., Hsieh C-H. 2018 . Fishing and temperature effects on the size structure of exploited fish stocks . Scientific Reports , 8 : 7132 . Google Scholar Crossref Search ADS PubMed WorldCat van der Maas H. L. J. , Kolstein R., van der Pligt J. 2003 . Sudden transitions in attitudes . Sociological Methods & Research , 32 : 125 – 152 . Google Scholar Crossref Search ADS WorldCat Wang H. Y. , Botsford L. W., White J. W., Fogarty M. J., Juanes F., Hastings A., Holland M. D., et al. 2014 . Effects of temperature on life history set the sensitivity to fishing in Atlantic cod Gadus morhua . Marine Ecology Progress Series , 514 : 217 – 229 . Google Scholar Crossref Search ADS WorldCat Ye H. , Beamish R. J., Glaser S. M., Grant S. C. H., Hsieh C., Richards L. J., Schnute J. T., et al. 2015 . Equation-free mechanistic ecosystem forecasting using empirical dynamic modeling . Proceedings of the National Academy of Sciences , 112 : E1569 – E1576 . Google Scholar Crossref Search ADS WorldCat Ye H. , Clark A., Deyle E. R., Sugihara G. 2016 . rEDM: an R package for empirical dynamic modelling and convergent cross-mapping. Ye H. , Sugihara G. 2016 . Information leverage in interconnected ecosystems: overcoming the curse of dimensionality . Science , 353 : 922 – 925 . Google Scholar Crossref Search ADS PubMed WorldCat Zeeman E. 1976 . Catastrophe theory . Scientific American , 234 : 65 – 83 . Google Scholar Crossref Search ADS WorldCat © International Council for the Exploration of the Sea 2019. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © International Council for the Exploration of the Sea 2019.
Responses of ecological indicators to fishing pressure under environmental change: exploring non-linearity and thresholdsFu, Caihong; Xu, Yi; Grüss, Arnaud; Bundy, Alida; Shannon, Lynne; Heymans, Johanna J; Halouani, Ghassen; Akoglu, Ekin; Lynam, Christopher P; Coll, Marta; Fulton, Elizabeth A; Velez, Laure; Shin, Yunne-Jai
doi: 10.1093/icesjms/fsz182pmid: N/A
Abstract Marine ecosystems are influenced by multiple stressors in both linear and non-linear ways. Using generalized additive models (GAMs) fitted to outputs from a multi-ecosystem, multi-model simulation experiment, we investigated 14 major ecological indicators across ten marine ecosystems about their responses to fishing pressure under: (i) three different fishing strategies (focusing on low-, high-, or all-trophic-level taxa); and (ii) four different scenarios of directional or random primary productivity change, a proxy for environmental change. From this work, we draw four major conclusions: (i) responses of indicators to fishing mortality in shapes, directions, and thresholds depend on the fishing strategies considered; (ii) most of the indicators demonstrate decreasing trends with increasing fishing mortality, with a few exceptions depending on the type of fishing strategy; (iii) most of the indicators respond to fishing mortality in a linear way, particularly for community and biomass-based indicators; and (iv) occurrence of threshold for non-linear-mixed type (i.e. non-linear with inflection points) is not prevalent within the fishing mortality rates explored. The conclusions drawn from the present study provide a knowledge base in indicators’ dynamics under different fishing and primary productivity levels, thereby facilitating the application of ecosystem-based fisheries management worldwide. Introduction Marine ecosystems are influenced by multiple stressors including climatic, oceanographic, ecological, and anthropogenic (e.g. fishing) changes. At the species level, these stressors may affect survival, growth, and reproduction, which can influence the structure and functioning of the entire ecosystem through complex species interactions. Ecosystem responses to stressors can be either linear or non-linear, depending, inter alia, on the properties of individual stressors and the way stressors interact. Particular attention needs to be paid to non-linear responses to stressors, in which thresholds (i.e. tipping points), may exist, beyond which large changes in ecosystem state or properties can occur abruptly (also described as regime shifts; Scheffer et al., 2001; Scheffer and Carpenter, 2003; Groffman et al., 2006; Andersen et al., 2009; Kelly et al., 2015; Moore, 2018). Identifying thresholds in non-linear responses to stressors can facilitate the implementation of risk-based management targets, the avoidance of detrimental ecosystem shifts and the optimization of socio-economic returns from ecosystem exploitation (Foley et al., 2015; Kelly et al., 2015). In the present study, two ubiquitous stressors are considered: fishing pressure and primary productivity change, where the latter serves as a proxy for environmental change. Fishing has been a major stressor in exploited marine ecosystems around the globe for the last several decades (Boldt et al., 2014), leading to substantial changes in the state of ecosystems (e.g. Jackson et al., 2001). In parallel, changes in ocean temperature, circulation patterns, vertical mixing, and availability of nutrients have resulted in changes in the primary productivity of marine ecosystems (e.g. Brown et al., 2011; Hunt et al., 2011; Poloczanska et al., 2016; Capuzzo et al., 2018). Fishing pressure and primary productivity may work synergistically, affecting the structure and function of marine ecosystems (Kirby et al., 2009; Möllmann et al., 2009) and how they should be managed under future climate change conditions (Serpetti et al., 2017; Baudron et al., 2019). To explore how an ecosystem responds (e.g. linearly vs. non-linearly) to fishing pressure and primary productivity change, various ecological indicators that represent ecosystem state and properties can be employed. Ecological indicators have been extensively studied using empirical data about their trends (e.g. Blanchard et al., 2010), relationships with fishing pressure and environmental changes (e.g. Link et al., 2010; Shannon et al., 2010; Shannon et al., 2014; Fu et al., 2015), and potential thresholds (e.g. Large et al., 2013, 2015; Tam et al., 2017). Empirical studies have suggested that past and present exploitation strategies, mechanisms of productivity dynamics, and dominant ecological and environmental traits are all likely to influence how ecological indicators respond to various stressors (Shannon et al., 2014). Consequently, it has long been a major challenge to draw general conclusions about the response patterns of indicators to fishing pressure and environmental changes. To overcome this major challenge, we used a comparative multi-ecosystem, multi-model simulation experiment to explore the differential responses of ecological indicators to fishing mortality and primary productivity change. We took advantage of the work conducted by the IndiSeas working group (Shin et al., 2012) in which researchers across different continents carried out unified simulation experiments for ten marine ecosystems around the globe. Prior to the simulation experiments, each of the ten ecosystems already had a working ecosystem simulation model developed (Shin et al., 2018; Fu et al., 2018, 2019), thus having overcome the time constraint and resource limitation of developing brand new ecosystem models. Relying on different ecosystem modelling approaches applied to different marine ecosystems allowed us to draw conclusions independent of the specific structure of ecosystems or the specific structure and assumptions of ecosystem simulation models. The present study primarily addresses the following questions: (i) How do ecological indicators respond to fishing mortality in trends and directions under different types of primary productivity change? (ii) How frequently do non-linear responses to fishing mortality occur in different ecological indicators? (3) How frequently do thresholds occur in non-linear indicators’ responses to fishing mortality? Results from the present study will provide insights on specific responses of indicators to fishing mortality; the thresholds of non-linear responses that we identify can inform the development of management strategies integrating environmental considerations, thereby facilitating the move towards ecosystem-based fisheries management worldwide. Material and methods Ecosystem models In the present study, ten marine ecosystems are considered. The structure and functioning of these ecosystems were simulated using one of four different ecosystem simulation modelling approaches: Ecopath with Ecosim (EwE; Christensen and Walters, 2004), OSMOSE (Shin and Cury, 2004), Atlantis (Fulton et al., 2004), or a multispecies size spectrum model (Blanchard et al., 2014). Five of the ten ecosystems were modelled with EwE: the Black Sea (Akoglu, 2013; Akoglu et al., 2014), the southern Benguela (Shannon et al., 2004, 2008), the southern Catalan Sea (Coll et al., 2006, 2013), the western Scotian Shelf (Araújo and Bundy, 2011, 2012), and western Scotland (Alexander et al., 2015). OSMOSE was employed to model three of ten ecosystems: the West Coast of Canada (Fu et al., 2013), the West Florida Shelf (Grüss et al., 2016a, 2016b), and the Gulf of Gabes, Tunisia (Halouani et al., 2016). Finally, the Southeastern Australian ecosystem was modelled with Atlantis (Fulton et al., 2014), whereas the North Sea was modelled with the size spectrum model (Blanchard et al., 2014). All the ecosystem models used in the present study have been published and validated against observed or estimated abundance, biomass and/or catch data. Details on the structure, assumptions, and species/taxa represented in the ten ecosystem models are provided in Supplementary Tables S1 and S2. Fishing mortality For each ecosystem, various levels of fishing mortality rate relative to FMSY (the fishing mortality rate corresponding to the maximum sustainable yield of a specific species/taxon within a specific ecosystem) were implemented to capture the impacts of different fishing levels. The same design was carried out consistently across all ecosystems for valid comparisons. Prior to running the experimental simulation within each ecosystem, FMSY was estimated for each exploited species/taxon, by constructing the yield to fishing mortality curve (fisheries yield as a function of fishing mortality rate) of that species/taxon while holding the fishing mortality rate of all other species/taxa constant at their respective current fishing mortality rates (Fcurr). Then, during the experimental simulation runs, fishing mortality rates were incrementally increased relative to FMSY for each species/taxon of interest through the use of a multiplier λ. In other words, fishing mortality rates Fi=λ×FMSYi (year−1) were applied to each species/taxon of interest i, where λ ∈ {0.25, 0.5, 0.75, 1, 1.25, 1.5}. This range of fishing mortality rates covers representative values for the yield to fishing mortality curves (Fu et al., 2018; Shin et al., 2018). For each ecosystem, we consider three different fishing strategies: a “low-trophic-level” (F_ltl), a “high-trophic-level” (F_htl), and an “all-trophic-level” (F_all) strategy. With the F_ltl fishing strategy, fishing experiments focus on all forage species/taxa retained in commercial or subsistence fisheries. Here, forage species/taxa are defined as species/taxa whose adults mainly feed on plankton (phyto-, zoo-, or ichthyoplankton) or similar low-trophic level organisms (e.g. meiobenthos). With the F_htl fishing strategy, fishing experiments focus on predatory taxa, which include large demersal and large pelagic taxa. Finally, the F_all fishing strategy represents broad-scale exploitation, where the focus taxa are all taxa retained in commercial or subsistence fisheries, including both high-trophic-level (HTL) and low-trophic-level (LTL) taxa. Any pre-recruit stages that are represented in the ecosystem models are excluded from the fishing scenarios, as are marine mammals, marine turtles, and seabirds. The HTL and LTL species/taxa represented in each of the ten ecosystem models are detailed in Supplementary Table S1. Primary productivity For all ten ecosystems, changes in phytoplankton biomass, comparable across models and ecosystems, were used to represent changes in primary productivity. For the Atlantis ecosystem model, however, it is not possible to directly scale phytoplankton biomass, because of the role of phytoplankton in the biogeochemical cycles represented in Atlantis. Therefore, in this case, time-series of nutrient inputs were scaled so that the resulting changes in phytoplankton biomass were in line with those implemented in other models. Two types of changes in primary productivity were simulated: directional and random. For directional changes in primary productivity, a multiplier γ ∈ {0.85, 0.9, 0.95, 1, 1.05, 1.1} was directly applied to modelled phytoplankton biomass. This range of variability encompasses the range of changes observed in the ten ecosystems in the last decade (Boyce et al., 2014). With respect to random changes in primary productivity, the modelled phytoplankton biomass was multiplied by a random multiplier drawn from a lognormal distribution with a mean μ of 1 and a range of standard deviations σ = {0.1, 0.2, 0.3}. These standard deviation values are consistent with the observed annual satellite-derived chlorophyll-a concentration levels from the MODIS (Moderate Resolution Imaging Spectroradiometer) Aqua spectral data (Shin et al., 2018). A set of 30 random multipliers was generated for each value of σ, so as to adequately sample the random distribution. Consequently, within each ecosystem, 12 combinations of fishing strategies (F_ltl, F_htl, and F_all) and primary productivity change (directional change, and random change while assuming σ = 0.1, 0.2, or 0.3) were simulated to explore interactions of these two factors. Ecological indicators Following Fu et al. (2019), 14 ecological indicators were explored (Table 1), including biomass to fisheries catch ratio (B/C), proportion of predatory fish (Pred), mean intrinsic vulnerability (IVI), mean lifespan (Lifesp), mean trophic level of the community (TLco), and marine trophic index (MTI), the trophic level (TL) of catch calculated assuming a constant TL for taxa (TLc), the TL of catch calculated assuming a variable TL for taxa (TLcVar; Reed et al., 2017), and six community-level biomass indicators [biomass of all-trophic-level taxa (B_all), the biomass of HTL taxa (B_htl), the biomass of LTL taxa (B_ltl), the ratio of B_htl to B_all (B_H2A), the ratio of B_ltl to B_all (B_L2A), and the ratio of B_ltl to B_htl (B_L2H)]. Aside from the six community-level biomass indicators, the indicators Pred, Lifesp, TLco were also derived from survey biomass data (i.e. biomass-based); in contrast, the indicators IVI, MTI, TLc, and TLcVar were calculated from catch data (i.e. catch-based). Within each of the 10 ecosystem models, these 14 ecological indicators were averaged over the last 10 years of simulation period under each scenario of varying fishing mortality rate and primary productivity, which implies that these indicator data points corresponding to different fishing mortality rates and primary productivity levels are independent of each other. Table 1. List of the indicators considered in this study. Indicator . Definition . Abbreviation . Biomass to catch ratio B/C B/C Proportion of predatory fish B (predatory fish)/ B (surveyed) Pred Mean intrinsic vulnerability ∑sIVIsCs∑sCs IVI Mean lifespan ∑sAgemax, sBs∑sBs Lifesp Trophic level of catch ∑sTLsCs∑sCs TLc Trophic level of catch with variable TL ∑sTL'sCs∑sCs TLcVar Marine trophic index ∑s(TL>3.25)TLsCs∑s(TL>3.25)Cs MTI Mean trophic level of community ∑sTLsBs∑sBs TLco Biomass of all surveyed species ∑s(all)Bs B_all Biomass of high-trophic-level (htl) species ∑s(htl)Bs B_htl Biomass of low-trophic-level (ltl) species ∑s(ltl)Bs B_ltl Ratio of htl biomass to total biomass B_htlB_all B_H2A Ratio of ltl biomass to total biomass B_ltlB_all B_L2A Ratio of ltl biomass to htl biomass BltlBhtl B_L2H Indicator . Definition . Abbreviation . Biomass to catch ratio B/C B/C Proportion of predatory fish B (predatory fish)/ B (surveyed) Pred Mean intrinsic vulnerability ∑sIVIsCs∑sCs IVI Mean lifespan ∑sAgemax, sBs∑sBs Lifesp Trophic level of catch ∑sTLsCs∑sCs TLc Trophic level of catch with variable TL ∑sTL'sCs∑sCs TLcVar Marine trophic index ∑s(TL>3.25)TLsCs∑s(TL>3.25)Cs MTI Mean trophic level of community ∑sTLsBs∑sBs TLco Biomass of all surveyed species ∑s(all)Bs B_all Biomass of high-trophic-level (htl) species ∑s(htl)Bs B_htl Biomass of low-trophic-level (ltl) species ∑s(ltl)Bs B_ltl Ratio of htl biomass to total biomass B_htlB_all B_H2A Ratio of ltl biomass to total biomass B_ltlB_all B_L2A Ratio of ltl biomass to htl biomass BltlBhtl B_L2H High-trophic-level and low-trophic-level taxa for all the ecosystems are listed in Supplementary Table S1. B, biomass (in tonnes); C, catch (in tonnes); s, species; TL, trophic level; TL', variable TL; IVI, intrinsic vulnerability index. Open in new tab Table 1. List of the indicators considered in this study. Indicator . Definition . Abbreviation . Biomass to catch ratio B/C B/C Proportion of predatory fish B (predatory fish)/ B (surveyed) Pred Mean intrinsic vulnerability ∑sIVIsCs∑sCs IVI Mean lifespan ∑sAgemax, sBs∑sBs Lifesp Trophic level of catch ∑sTLsCs∑sCs TLc Trophic level of catch with variable TL ∑sTL'sCs∑sCs TLcVar Marine trophic index ∑s(TL>3.25)TLsCs∑s(TL>3.25)Cs MTI Mean trophic level of community ∑sTLsBs∑sBs TLco Biomass of all surveyed species ∑s(all)Bs B_all Biomass of high-trophic-level (htl) species ∑s(htl)Bs B_htl Biomass of low-trophic-level (ltl) species ∑s(ltl)Bs B_ltl Ratio of htl biomass to total biomass B_htlB_all B_H2A Ratio of ltl biomass to total biomass B_ltlB_all B_L2A Ratio of ltl biomass to htl biomass BltlBhtl B_L2H Indicator . Definition . Abbreviation . Biomass to catch ratio B/C B/C Proportion of predatory fish B (predatory fish)/ B (surveyed) Pred Mean intrinsic vulnerability ∑sIVIsCs∑sCs IVI Mean lifespan ∑sAgemax, sBs∑sBs Lifesp Trophic level of catch ∑sTLsCs∑sCs TLc Trophic level of catch with variable TL ∑sTL'sCs∑sCs TLcVar Marine trophic index ∑s(TL>3.25)TLsCs∑s(TL>3.25)Cs MTI Mean trophic level of community ∑sTLsBs∑sBs TLco Biomass of all surveyed species ∑s(all)Bs B_all Biomass of high-trophic-level (htl) species ∑s(htl)Bs B_htl Biomass of low-trophic-level (ltl) species ∑s(ltl)Bs B_ltl Ratio of htl biomass to total biomass B_htlB_all B_H2A Ratio of ltl biomass to total biomass B_ltlB_all B_L2A Ratio of ltl biomass to htl biomass BltlBhtl B_L2H High-trophic-level and low-trophic-level taxa for all the ecosystems are listed in Supplementary Table S1. B, biomass (in tonnes); C, catch (in tonnes); s, species; TL, trophic level; TL', variable TL; IVI, intrinsic vulnerability index. Open in new tab Generalized additive models: linearity/non-linearity and directions Nonparametric regression methods such as generalized additive models (GAMs; Hastie and Tibshirani, 1990) are powerful tools for exploring linear or non-linear responses of a variable to predictors without being constrained to an underlying parametric model of a specific form, which is particularly useful when ecological thresholds of non-linear responses are of interest (Lintz et al., 2011; Nicolaou and Constandinou, 2016). Therefore, GAMs were employed in this study to explore the responses of ecological indicators to fishing mortality and primary productivity change. Each of the two stressors was modelled individually in order to avoid convergence issues because of the small number of fishing mortality and primary productivity levels considered in this study, as well as to avoid potential correlation between fishing mortality and primary productivity. The GAM model is formulated as (Large et al., 2013): Y=α+sX+ɛ,(1) where Y is the ecological indicator; α is the intercept held constant; X is the predictor, representing either fishing mortality or primary productivity level; s() is the smoothing function; and ɛ is the error term of a normal distribution with zero mean and finite variance. The smoothing parameter was estimated using the generalized cross-validation criterion (GCV criterion; Wood et al., 2016) implemented in R package “mgcv” (Version 1.8-24; Wood, 2017). Thin plate regression spline smoothing terms with a ridge penalty were employed, such that the smoothing term could be eliminated and the predicted curve could be reduced to a linear function when the smoothing term did not improve fit. Confidence intervals were calculated for each GAM with a naive bootstrap using random sampling and replacement [see Large et al. (2013) for more details on the naive bootstrap approach]. For the purpose of exploring linear or non-linear responses of ecological indicators to fishing mortality, the present study focuses on the development and exploration of 10 × 14 × 3 × 4 = 1680 GAMs, covering 10 ecosystems, 14 ecological indicators, 3 fishing strategies (F_ltl, F_htl, and F_all), and 4 scenarios of primary productivity change (directional change, and random change while assuming σ = 0.1, 0.2, or 0.3). The performance of a GAM model is often gauged by examining the amount of variance it explains (i.e. the proportion of the deviance in the data explained). Comparisons of GAM predictions across the different scenarios of primary productivity change provide perspectives on how the nature and degree of environmental variability affects the capacity of GAMs to capture the detrimental impacts of fishing, as well as whether indicators’ responses to fishing pressure are robust to varying levels of environmental variability. Predicted GAM curves based on empirical time-series have often been explored for non-linear responses to drivers using first and second derivatives (e.g. Large et al., 2013; Burthe et al., 2016). However, empirical time-series are typically faced with statistical problems such as temporal autocorrelation, confounding among unknown sources of variability and inadequacy in replication, resulting in difficulty in identifying linearity/non-linearity (Litzow and Hunsicker, 2016). In contrast to empirical time-series data explored in Large et al. (2013) and Burthe et al. (2016), the ecological indicator data in the present study were obtained through individual simulations with different fishing mortality and primary productivity and are thus independent data points and not time-series, as stated above. Thus, the linearity/non-linearity was explored in a more literal sense using specifically designed methods for these datasets, as described below. For each of the 1680 GAMs with respect to fishing mortality, we first determine whether the indicator’s response to fishing mortality is linear by fitting a linear regression to the GAM predicted values. If the correlation coefficient between the GAM predicted values and the linear regression fitted values is >0.95, the indicator’s response to fishing mortality (derived from the GAM) is considered linear. Otherwise, the indicator’s response to fishing mortality is assumed to be non-linear. If the relationship between the indicator and fishing mortality is linear and if the slope of the linear regression is statistically significant (p < 0.05), the direction of indicator’s change in response to fishing mortality is defined as a “linear increase” when the slope is >0 and a “linear decrease” when the slope is <0. When the relationship between the indicator and fishing mortality is non-linear, the GAM curve is further investigated for identifying monotonicity and determining the direction of indicator’s change in response to fishing mortality by examining the first derivative of the curve. A “non-linear-monotonic increase” is associated with all first derivatives being non-negative; whereas a “non-linear-monotonic decrease” is associated with all first derivatives being non-positive. However, if the first derivatives have both positive and negative values, the GAM curve is defined as non-linear-mixed. For this type of non-linearity, the direction of an indicator’s change in response to fishing mortality is determined by fitting a linear regression to GAM predicted values. When a significant linear relationship exists (i.e. the correlation coefficient between the GAM predicted values and the linear regression fitted values is >0.95), the non-linear-mixed curve is defined as a “non-linear-mixed increase” when the linear regression fitted slope is >0 and as a “non-linear-mixed decrease” when the linear regression fitted slope is <0. Thresholds For non-linear-mixed GAM curves, there are regions along the fishing mortality gradient where the first derivative of the response function changes sign at an inflection point. Inflection points on GAM curves have been considered as ecological thresholds (e.g. Large et al., 2013; Tam et al., 2017), and such ecological thresholds, also known as tipping points (e.g. Moore, 2018), may have implication for fisheries management (e.g. indicate the occurrence of a regime shift) and, therefore, deserve special attention. Because of the fact that there are only six discrete fishing mortality rates (FMSY multiplier = 0.25, 0.5, 0.75, 1.0, 1.25, 1.5), we do not use the threshold-determination approach developed by Large et al. (2013), rather, we fit a fifth order polynomial model [degree of freedom = 6, Equation (1)] in all the cases where the indicator’s response to fishing mortality is identified as non-linear-mixed: fx=a0+a1x+a2x2+a3x3+a4x4+a5x5,(2) where x is the predictor, representing the levels of fishing mortality; fx is the ecological indicator as a function of x ; and a0, ⋯, a5 are coefficients of the polynomial model. The analytical root(s) of the second derivatives, where ecological thresholds reside, are derived from: f''(x)=2a2+6a3x+12a4x2+20a5x3=0.(3) The analytical root(s) (x = xi) are solved using the “rootSolve” package in R (R Core Team, 2019). Only root(s) within the range of FMSY multiplier (0.25∼1.25) are retained for consideration. In addition, root(s) with their predicted values f(xi) being outside the range defined by the indicator value on either side of the xi value are not considered as thresholds, as they are inflection points resulting from local maxima or minima artificially introduced by high-order polynomial fits. Results Variance explained The amount of variance explained by fishing mortality (i.e. RF2 ), though varying among the ten ecosystems under each of the three fishing strategies (Figures 1–3), has some common themes. First, the RF2 values for all indicators tend to be the highest under the F_all fishing strategy with all four scenarios of primary productivity change across all ecosystems, whereas those under the F_ltl fishing strategy tend to be the lowest (the higher the RF2 values, the stronger the indication of fishing impacts). Within a given ecosystem, an individual indicator can have drastically different RF2 values among the three fishing strategies, implying that fishing patterns are important when interpreting the dynamics of ecological indicators. Second, the RF2 value of a specific indicator under any fishing strategy is generally smaller under directional than under random primary productivity change, and the RF2 value consistently decreases with increasing primary productivity variability (i.e. when σ increases from 0.1 to 0.3), suggesting that an indicator’s response to fishing mortality can be increasingly obscured when environmental variability increases. Third, comparing across the 14 indicators, the biomass indicators B_all, B_htl, and B_ltl tend to have the lowest RF2 values under each of the three fishing strategies, which can be more clearly seen from their values averaged over the ten ecosystems (Figure 4). These ecosystem-averaged RF2 values indicate the ability of the 14 indicators to reflect fishing impacts. Fourth, the RF2 value for the indicator B/C under the F_all fishing strategy is close to one under all four scenarios of primary productivity change for all ecosystems except the southern Catalan Sea, indicating that the trends for the indicator B/C are usually well explained by fishing mortality when both HTL and LTL taxa are targeted by fishing. Figure 1. Open in new tabDownload slide Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators for the 10 study marine ecosystems under the low-trophic-level (F_ltl) fishing strategy under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1. Figure 1. Open in new tabDownload slide Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators for the 10 study marine ecosystems under the low-trophic-level (F_ltl) fishing strategy under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1. Figure 2. Open in new tabDownload slide Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators for the 10 study marine ecosystems under the high-trophic-level (F_htl) fishing strategy under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1. Figure 2. Open in new tabDownload slide Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators for the 10 study marine ecosystems under the high-trophic-level (F_htl) fishing strategy under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1. Figure 3. Open in new tabDownload slide Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators in the 10 study marine ecosystems under the all-trophic-level (F_all) fishing strategy under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1. Figure 3. Open in new tabDownload slide Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators in the 10 study marine ecosystems under the all-trophic-level (F_all) fishing strategy under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1. Figure 4. Open in new tabDownload slide Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators for the 10 study marine ecosystems that were averaged over all the 10 ecosystems under each of the 3 fishing strategies: (top) F_ltl (the low-trophic-level fishing strategy), (middle) F_htl (the high-trophic-level fishing strategy), and (bottom) F_all (the all-trophic-level fishing strategy) under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1. Figure 4. Open in new tabDownload slide Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators for the 10 study marine ecosystems that were averaged over all the 10 ecosystems under each of the 3 fishing strategies: (top) F_ltl (the low-trophic-level fishing strategy), (middle) F_htl (the high-trophic-level fishing strategy), and (bottom) F_all (the all-trophic-level fishing strategy) under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1. Responses of indicators to fishing mortality: linearity/non-linearity Overall, the shapes of the GAM curves (linear, non-linear-monotonic, or non-linear-mixed) for most indicators are consistent among the four scenarios of primary productivity change (Figures 5–7), suggesting that most indicators’ responses to fishing mortality are largely insensitive to primary productivity variability. In contrast, the shapes of the GAM curves for most indicators across the ten ecosystems tend to vary with the three different fishing strategies, indicating the importance of fishing patterns when considering the responses of ecological indicators to fishing mortality. Figure 5. Open in new tabDownload slide Schematic plots of the shapes, including linear (shown as line), non-linear-monotonic (solid curve), and non-linear-mixed (dashed curve), and directions, including increase (upward arrowhead), decrease (downwards arrowhead), and no significant trend (no arrowhead), of the generalized additive models (GAMs) for each of the 14 indicators for the 10 marine ecosystems under the low-trophic-level (F_ltl) fishing strategy using fishing mortality as the predictor under 4 different scenarios of primary productivity (directional, random changes with σ = 0.1, σ = 0.2, and σ = 0.3, respectively). The locations of thresholds are shown as dark dots and those for inflection points that do not constitute thresholds are shown as light dots. The full name and meaning of ecological indicators are provided in Table 1. Figure 5. Open in new tabDownload slide Schematic plots of the shapes, including linear (shown as line), non-linear-monotonic (solid curve), and non-linear-mixed (dashed curve), and directions, including increase (upward arrowhead), decrease (downwards arrowhead), and no significant trend (no arrowhead), of the generalized additive models (GAMs) for each of the 14 indicators for the 10 marine ecosystems under the low-trophic-level (F_ltl) fishing strategy using fishing mortality as the predictor under 4 different scenarios of primary productivity (directional, random changes with σ = 0.1, σ = 0.2, and σ = 0.3, respectively). The locations of thresholds are shown as dark dots and those for inflection points that do not constitute thresholds are shown as light dots. The full name and meaning of ecological indicators are provided in Table 1. Figure 6. Open in new tabDownload slide Schematic plots of the shapes, including linear (shown as line), non-linear-monotonic (solid curve), and non-linear-mixed (dashed curve), and directions, including increase (upward arrowhead), decrease (downwards arrowhead), and no significant trend (no arrowhead), of the generalized additive models (GAMs) for each of the 14 indicators for the 10 marine ecosystems under the high-trophic-level (F_htl) fishing strategy using fishing mortality as the predictor variable under 4 different scenarios of primary productivity (directional, random changes with σ = 0.1, σ = 0.2, and σ = 0.3, respectively). The locations of thresholds are shown as dark dots and those for inflection points that do not constitute thresholds are shown as light dots. The full name and meaning of ecological indicators are provided in Table 1. Figure 6. Open in new tabDownload slide Schematic plots of the shapes, including linear (shown as line), non-linear-monotonic (solid curve), and non-linear-mixed (dashed curve), and directions, including increase (upward arrowhead), decrease (downwards arrowhead), and no significant trend (no arrowhead), of the generalized additive models (GAMs) for each of the 14 indicators for the 10 marine ecosystems under the high-trophic-level (F_htl) fishing strategy using fishing mortality as the predictor variable under 4 different scenarios of primary productivity (directional, random changes with σ = 0.1, σ = 0.2, and σ = 0.3, respectively). The locations of thresholds are shown as dark dots and those for inflection points that do not constitute thresholds are shown as light dots. The full name and meaning of ecological indicators are provided in Table 1. Figure 7. Open in new tabDownload slide Schematic plots of the shapes, including linear (shown as line), non-linear-monotonic (solid curve), and non-linear-mixed (dashed curve), and directions, including increase (upward arrowhead), decrease (downwards arrowhead), and no significant trend (no arrowhead), of the generalized additive models (GAMs) for each of the 14 indicators for the 10 marine ecosystems under the all-trophic-level (F_all) fishing strategy using fishing mortality as the predictor variable under 4 different scenarios of primary productivity (directional, random changes with σ = 0.1, σ = 0.2, and σ = 0.3, respectively). The locations of thresholds are shown as dark dots and those for inflection points that do not constitute thresholds are shown as light dots. The full name and meaning of ecological indicators are provided in Table 1. Figure 7. Open in new tabDownload slide Schematic plots of the shapes, including linear (shown as line), non-linear-monotonic (solid curve), and non-linear-mixed (dashed curve), and directions, including increase (upward arrowhead), decrease (downwards arrowhead), and no significant trend (no arrowhead), of the generalized additive models (GAMs) for each of the 14 indicators for the 10 marine ecosystems under the all-trophic-level (F_all) fishing strategy using fishing mortality as the predictor variable under 4 different scenarios of primary productivity (directional, random changes with σ = 0.1, σ = 0.2, and σ = 0.3, respectively). The locations of thresholds are shown as dark dots and those for inflection points that do not constitute thresholds are shown as light dots. The full name and meaning of ecological indicators are provided in Table 1. Across the ten ecosystems, the response to fishing mortality is largely linear (over 50% for most indicators; Figure 8). Within a given ecosystem, the shapes of GAM curves tend to differ among the 14 indicators. The catch-based indicators (e.g. B/C, IVI, TLc, MTI) tend to have larger proportion of non-linear responses. In contrast, the community and biomass-based indicators (e.g. TLco, B_all, B_htl, B_ltl) are largely linear in their response to fishing mortality (Figure 8), which may suggest that the responses of the community and biomass-based indicators to fishing mortality are more predictable. Figure 8. Open in new tabDownload slide Proportion of non-linear GAM (generalized additive model) types (including monotonic and mixed) for each of the 14 indicators averaged over the 10 marine ecosystems under each of the 3 fishing strategies (F_ltl for the low-trophic-level fishing strategy, F_htl for the high-trophic-level fishing strategy, and F_all for the all-trophic-level fishing strategy) under 2 types of primary productivity changes: (top) directional, and (bottom) random (the results for random were averaged over 3 levels of variability with σ = 0.1, 0.2, and 0.3). The full name and meaning of ecological indicators are provided in Table 1. Figure 8. Open in new tabDownload slide Proportion of non-linear GAM (generalized additive model) types (including monotonic and mixed) for each of the 14 indicators averaged over the 10 marine ecosystems under each of the 3 fishing strategies (F_ltl for the low-trophic-level fishing strategy, F_htl for the high-trophic-level fishing strategy, and F_all for the all-trophic-level fishing strategy) under 2 types of primary productivity changes: (top) directional, and (bottom) random (the results for random were averaged over 3 levels of variability with σ = 0.1, 0.2, and 0.3). The full name and meaning of ecological indicators are provided in Table 1. Responses of indicators to fishing mortality: direction Similar to the shapes of GAM curves, one common theme across the ten ecosystems is that the response direction of GAM curves (“decrease” or “increase”) for a specific indicator depends on the type of fishing strategy (Figures 5–7 and 9). Under the F_ltl fishing strategy, the indicators Pred, Lifesp, and B_H2A show increasing GAM trends in more than half of the study ecosystems under both directional and random primary productivity changes. However, these three indicators all show decreasing GAM trends in most ecosystems under the F_htl and F_all fishing strategies. Under the F_htl fishing strategy, there are six indicators (IVI, TLc, MTI, B_ltl, B_L2A, and B_L2H) that show an increasing trend with increasing fishing mortality in over 50% of the ecosystems under both directional and random primary productivity changes. However, these six indicators predominantly decrease with increasing fishing mortality under both the F_ltl and F_all fishing strategies. Under the F_all fishing strategy, the indicator B_L2H is the only one that predominantly shows an increasing GAM trend with increasing fishing mortality, demonstrating that the response direction of all other indicators becomes more consistent (i.e. always decreasing with increasing fishing mortality) when fishing targets both HTL and LTL taxa (i.e. under the F_all fishing strategy). It is worth noting that, among all the indicators, B/C is the only one that predominantly exhibits a decreasing trend with increasing fishing mortality under both directional and random primary productivity changes and under all three fishing strategies, which indicates that B/C is the most informative indicator of fishing impacts. Figure 9. Open in new tabDownload slide Proportion of GAM (generalized additive model) directions (increase and decrease) for each of the 14 indicators averaged over the 10 marine ecosystems under each of 3 fishing strategies (F_ltl for the low-trophic-level fishing strategy, F_htl for the high-trophic-level fishing strategy, and F_all for the all-trophic-level fishing strategy) under 2 types of primary productivity changes: (top) directional, and (bottom) random (the results for random were averaged over 3 levels of variability with σ = 0.1, 0.2, and 0.3). The full name and meaning of ecological indicators are provided in Table 1. Figure 9. Open in new tabDownload slide Proportion of GAM (generalized additive model) directions (increase and decrease) for each of the 14 indicators averaged over the 10 marine ecosystems under each of 3 fishing strategies (F_ltl for the low-trophic-level fishing strategy, F_htl for the high-trophic-level fishing strategy, and F_all for the all-trophic-level fishing strategy) under 2 types of primary productivity changes: (top) directional, and (bottom) random (the results for random were averaged over 3 levels of variability with σ = 0.1, 0.2, and 0.3). The full name and meaning of ecological indicators are provided in Table 1. Thresholds Overall, under all three fishing strategies, the occurrence of threshold for the non-linear-mixed type is not prevalent for most ecosystems (Figures 5–7; Supplementary Figures S1–S3). Exceptions include the Black Sea (under the F_ltl fishing strategy), Catalan Sea (under the F_htl fishing strategy), and Southeastern Australia and West Coast of Canada (under the F_all fishing strategy), where more than 50% of the non-linear-mixed GAM curves have thresholds (Table 2). The infrequency of threshold occurrences may imply that the root(s) of second derivatives either fall outside the fishing mortality range (0.25–1.5*FMSY) for most non-linear-mixed type (aside from the possibility of imaginary roots), indicating that future simulations may desire a broader range of fishing mortality rate in order to avoid missing true ecological thresholds. Table 2. Total number of non-linear-mixed GAM (generalized additive model) curves and identified thresholds for each of the ten ecosystems and under each of the three fishing strategies (F_ltl, F_htl, F_all), as well as the percentage of non-linear-mixed GAM curves that have thresholds. Ecosystem model . Number of non- linear-mixed GAM . Number of thresholds . Number of fake thresholds . Percentage of threshold occurrences . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . BlackSea-EwE 21 8 38 19 5 18 1 0 4 90 63 47 GulfGabes-OSMOSE 16 22 9 8 2 2 0 0 0 50 9 22 NorthSea-SizeS 24 7 11 11 4 7 0 0 0 46 57 64 CatalanSea.S-EwE 22 23 26 9 15 5 0 0 0 41 65 19 Australian.SE-Atlantis 26 12 46 8 0 24 0 0 16 31 0 52 Benguela.S-EwE 9 4 16 6 1 4 0 0 1 67 25 25 Canada.W-OSMOSE 6 6 20 2 4 19 0 0 0 33 67 95 Scotland.W-EwE 12 22 7 1 3 1 1 0 0 8 14 14 FloridaShelf.W-OSMOSE 7 10 12 1 5 4 0 0 0 14 50 33 ScotianShelf.W-EwE 21 21 4 0 2 3 0 0 1 0 10 75 Ecosystem model . Number of non- linear-mixed GAM . Number of thresholds . Number of fake thresholds . Percentage of threshold occurrences . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . BlackSea-EwE 21 8 38 19 5 18 1 0 4 90 63 47 GulfGabes-OSMOSE 16 22 9 8 2 2 0 0 0 50 9 22 NorthSea-SizeS 24 7 11 11 4 7 0 0 0 46 57 64 CatalanSea.S-EwE 22 23 26 9 15 5 0 0 0 41 65 19 Australian.SE-Atlantis 26 12 46 8 0 24 0 0 16 31 0 52 Benguela.S-EwE 9 4 16 6 1 4 0 0 1 67 25 25 Canada.W-OSMOSE 6 6 20 2 4 19 0 0 0 33 67 95 Scotland.W-EwE 12 22 7 1 3 1 1 0 0 8 14 14 FloridaShelf.W-OSMOSE 7 10 12 1 5 4 0 0 0 14 50 33 ScotianShelf.W-EwE 21 21 4 0 2 3 0 0 1 0 10 75 The scenarios with more than 50% chance of having thresholds are bolded, excluding those that have <19 (one-third of 56) non-linear-mixed GAM curves within a particular ecosystem. Open in new tab Table 2. Total number of non-linear-mixed GAM (generalized additive model) curves and identified thresholds for each of the ten ecosystems and under each of the three fishing strategies (F_ltl, F_htl, F_all), as well as the percentage of non-linear-mixed GAM curves that have thresholds. Ecosystem model . Number of non- linear-mixed GAM . Number of thresholds . Number of fake thresholds . Percentage of threshold occurrences . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . BlackSea-EwE 21 8 38 19 5 18 1 0 4 90 63 47 GulfGabes-OSMOSE 16 22 9 8 2 2 0 0 0 50 9 22 NorthSea-SizeS 24 7 11 11 4 7 0 0 0 46 57 64 CatalanSea.S-EwE 22 23 26 9 15 5 0 0 0 41 65 19 Australian.SE-Atlantis 26 12 46 8 0 24 0 0 16 31 0 52 Benguela.S-EwE 9 4 16 6 1 4 0 0 1 67 25 25 Canada.W-OSMOSE 6 6 20 2 4 19 0 0 0 33 67 95 Scotland.W-EwE 12 22 7 1 3 1 1 0 0 8 14 14 FloridaShelf.W-OSMOSE 7 10 12 1 5 4 0 0 0 14 50 33 ScotianShelf.W-EwE 21 21 4 0 2 3 0 0 1 0 10 75 Ecosystem model . Number of non- linear-mixed GAM . Number of thresholds . Number of fake thresholds . Percentage of threshold occurrences . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . F_ltl . F_htl . F_all . BlackSea-EwE 21 8 38 19 5 18 1 0 4 90 63 47 GulfGabes-OSMOSE 16 22 9 8 2 2 0 0 0 50 9 22 NorthSea-SizeS 24 7 11 11 4 7 0 0 0 46 57 64 CatalanSea.S-EwE 22 23 26 9 15 5 0 0 0 41 65 19 Australian.SE-Atlantis 26 12 46 8 0 24 0 0 16 31 0 52 Benguela.S-EwE 9 4 16 6 1 4 0 0 1 67 25 25 Canada.W-OSMOSE 6 6 20 2 4 19 0 0 0 33 67 95 Scotland.W-EwE 12 22 7 1 3 1 1 0 0 8 14 14 FloridaShelf.W-OSMOSE 7 10 12 1 5 4 0 0 0 14 50 33 ScotianShelf.W-EwE 21 21 4 0 2 3 0 0 1 0 10 75 The scenarios with more than 50% chance of having thresholds are bolded, excluding those that have <19 (one-third of 56) non-linear-mixed GAM curves within a particular ecosystem. Open in new tab The occurrence of inflection points resulting from local maxima or minima of the fifth order polynomial fits (not signalling thresholds) is rare except in the Southeastern Australia ecosystem, where there are 16 such type of inflection points (consisting 35% of the non-linear-mixed type contrasting with 52% of thresholds; Table 2). This implies that thresholds identified within the fishing mortality range are generally appropriate. Discussion The multi-ecosystem, multi-model simulation experiment conducted within the IndiSeas programme (Shin et al., 2018; Fu et al., 2019) has generated a comprehensive simulation database of ecological indicators under two ecosystem stressors, fishing pressure, and primary productivity change. The GAMs employed in the present study have allowed us to produce response curves for 14 major ecological indicators evaluated for 10 ecosystems in response to fishing mortality, revealing linearity or non-linearity in the responses, with decreasing or increasing trends and thresholds associated with non-linear responses. This comprehensive and comparative research thus provides a rich knowledge base of commonalities and differences among ecosystems in indicators’ dynamics and performance at different levels of fishing mortality under different primary productivity scenarios, as well as thresholds associated with non-linear responses of indicators, thereby providing essential knowledge to facilitate the move towards ecosystem-based fisheries management. Comparison with previous work Using the gradient forest method, Fu et al. (2019) investigated the performance of the same set of ecological indicators about their sensitivity, specificity, and threshold responses to both fishing pressure and primary productivity change. However, the gradient forest method does not allow quantitative descriptions of the directions (“decrease” or “increase”) and shapes of the indicators’ responses to pressures (e.g. fishing mortality). Samhouri et al. (2017) employed a multi-model inference framework, a combination of the gradient forest method and GAMs, for exploring directions, shapes, and thresholds in the responses of ecosystem state indicators to pressures of environmental change and human activities. The present study also uses GAMs to expand the previous work with the gradient forest method (Fu et al., 2019). Additionally, the present study has expanded the investigation from one fishing strategy (F_all) in Fu et al. (2019) to all three fishing strategies (F_ltl, F_htl, F_all). The expansion has been proven to be valuable as it has led to the important discovery that for all ten ecosystems, all aspects of the GAM curves, including model performances in the form of variance explained (Figures 1–3), shapes (Figures 5–7), and the response directions (Figure 9) are dependent on the different fishing strategies considered. Despite the importance of considering fisheries exploitation history, as pointed out by Shannon et al. (2014) and Fu et al. (2018), ecological indicators have hitherto been investigated largely in the absence of such knowledge (e.g. Blanchard et al., 2010; Coll et al., 2016). Future data collection and analysis related to fisheries exploitation should be more conscious about fishing patterns, i.e. what trophic levels fisheries target. Responses of indicators to fishing mortality: variance explained Our results demonstrate that responses of indicators to fishing mortality differ across fishing strategies in all ten ecosystems (Figures 1–3, 5–7, and 9). In particular, the RF2 values for all indicators tend to be highest under the F_all fishing strategy whereas lowest under the F_ltl fishing strategy, which may suggest that the impacts of fisheries exploitation on LTL taxa could be obscured because of the relatively closer coupling between the dynamics of LTL taxa and primary productivity. In order to further understand the commonalities and differences among the ten ecosystems, we explored how they would respond to primary productivity change. We conducted an additional 1680 GAMs (10 ecosystems × 14 ecological indicators × 3 fishing strategies × 4 scenarios of primary productivity change) with primary productivity as the predictor. Unlike the GAMs with fishing mortality as the predictor, the GAMs with primary productivity as the predictor generally resulted in low RP2 under all fishing strategies across all ecosystems with only a few exceptions (Supplementary Figures S10–S12), which add insight to the results for fishing mortality. One interesting exception is the southern Benguela ecosystem under the F_ltl fishing strategy (Supplementary Figure S10). The southern Benguela is an upwelling ecosystem with abundant lower trophic level taxa. Upwelling indices of this ecosystem have shown increased variability in the 1990s and 2000s, presumably reflecting climate change (Blamey et al., 2012), which have caused ecosystem changes mediated through the lower trophic level fish community of this wasp-waist foodweb (Cury et al., 2000). The exceptionally high RP2 values (but low RF2 values) for most indicators under the F_ltl fishing strategy suggests that this ecosystem is extremely sensitive to changes in primary productivity (see Ortega-Cisneros et al., 2018), and any potential ecosystem change evoked through lower trophic level fisheries is primarily driven by environmental changes. This conclusion is consistent with previous indicator-based studies that have concluded that interpretation of ecological indicator trends in response to fishing in the southern Benguela requires careful unpacking of changes in environmental conditions in this ecosystem (Shannon et al., 2010; Lockerbie et al., 2016). The importance of environmental drivers in upwelling systems has been highlighted in the northern Benguela ecosystem, where environmental conditions and overfishing together have caused apparent regime shifts (Heymans and Tomczak, 2016). The Black Sea ecosystem also displays unique features, particularly the large discrepancies in the GAM curves between the directional and random primary productivity change scenarios indicating that this ecosystem is very sensitive (or responsive) to changes in primary productivity, more so than to changes in fishing mortality, when the primary productivity changes are directional. Similar to the southern Benguela ecosystem, the changes in the Black Sea ecosystem are predominantly mediated by the wasp-waist control exerted on the foodweb by small pelagic fish with a significant interplay of bottom-up control (Akoglu et al., 2014), which determines the endurance of small pelagic fish under the influence of fisheries overexploitation. However, the large discrepancies disappear when comparing the results of the three different scenarios of random primary productivity changes (where σ = 0.1, 0.2, or 0.3), which suggests that the indicators considered in the present study are robust to environmental variability for the Black Sea ecosystem. The West Florida Shelf ecosystem contrasts with the southern Benguela ecosystem in that the West Florida Shelf has very high RF2 values for most indicators and very low RP2 values for all indicators. These results concur with the findings of Grüss et al. (2016a, 2016b), which show that the West Florida Shelf ecosystem is under such a strong top-down control that elevated natural mortality because of large environmental changes (e.g. harmful algal blooms) has negligible impacts on ecosystem structure and functioning compared with elevated mortality resulting directly or indirectly (via trophic interactions) from changes in fishing patterns. Other ecosystems such as the western Scotian Shelf and West Coast Scotland have an intermediate response. Although the community-level biomass indicators (B_ltl, B_htl, B_all) tend to have higher RP2 values compared with RF2 , other indicators, such as B/C and IVI, indicate strong responsiveness to fishing mortality. These ecosystems have a long exploitation history; however, their temperate location may be providing some resilience to exploitation (Frank et al., 2007). In contrast to other ecosystems in this study, the southern Catalan Sea stands out as one that has quite different indicator dynamics in response to fishing mortality and primary productivity change. This ecosystem has a long fishing history, and current fishing mortality rates in the southern Catalan Sea (northwestern Mediterranean Sea) are very high (Fernandes et al., 2017). All indicators except for Lifesp, TLc, TLcVar, and MTI have very low RF2 values under both directional and random primary productivity changes. Essentially, the southern Catalan Sea ecosystem has been fished down and the foodweb of this ecosystem is currently under bottom-up control. Consequently, the variations in many of the ecological indicators are no longer responsive to changes in fishing pressure, and are instead primarily driven by directional changes in primary productivity with high RP2 for most indicators when there is directional primary productivity change, as already highlighted in previous studies (Shannon et al., 2014; Coll et al., 2016; Lockerbie et al., 2017). As shown above, the five ecosystems (Black Sea, southern Benguela, southern Catalan Sea, western Scotian Shelf, and western Scotland) demonstrate diverse responses of indicators to fishing mortality and primary productivity change, despite the same ecosystem modelling approach EwE applied to these ecosystems. This suggests that the types of ecosystem modelling approach may not contribute to the differing results among the ten ecosystems. Because of the commonly low RP2 of the GAMs with primary productivity as the predictor, there is no need to explore other properties of the GAMs, such as shapes and directions. However, the generally low RP2 values could be caused by the rather narrow range of primary productivity (multiplier = 0.85, 0.9, 0.95, 1, 1.05, 1.1). Future research with a broader range of primary productivity scenarios are warranted for investigating the responses of ecological indicators to extreme climate change (primary productivity change in this case). Responses of indicators to fishing mortality: direction The majority of the GAM curves demonstrates decreasing trends with increasing fishing mortality (Figures 5–7 and 9; Supplementary Figures S4–S9). These decreasing trends reflect the expected behaviour of indicators in response to increasing fishing mortality (Rice and Rochet, 2005; Bundy et al., 2010; Shin et al., 2010). However, the trends (“decrease” or “increase”) of the GAM curves of a specific indicator can depend on the type of fishing strategy considered. For instance, the indicators Pred, Lifesp, and B_H2A increase with increasing fishing mortality under the F_ltl fishing strategy, but decrease under the F_htl and F_all fishing strategies. This implies that caution is required when interpreting the dynamics of these indicators when fishing primarily targets lower trophic level taxa, and even more so in ecosystems such as upwelling systems (e.g. the southern Benguela) where environmental variability strongly influences ecological indicators. On the other hand, other indicators, including IVI, TLc, MTI, B_ltl, B_L2A, and B_L2H, increase with increasing fishing mortality under the F_htl fishing strategy, but decrease under F_ltl and F_all in most ecosystems. The decrease of indicator TLc in response to fishing has been observed in many ecosystems (Pauly and Watson, 2005; Fu et al., 2012), resulting from “fishing down the foodweb.” It is likely that the ecosystems for which these studies were carried out are ecosystems where fishing activities primarily target HTL taxa. However, the present study suggests that the indicator TLc (along with IVI, TLc, MTI, B_ltl, B_L2A, and B_L2H) may increase in response to increasing fishing mortality if fishing activities in the ecosystem of interest primarily target LTL taxa, providing potentially false assurance of ecosystem health. Therefore, our multi-ecosystem, multi-model simulation experiment supports the assertion of Branch (2015) and Shannon et al. (2014), i.e. multiple working hypotheses are needed to explore how fishing affects marine foodwebs rather than assuming “fishing down the foodweb” applies in all cases. In contrast to the situations under the F_ltl and F_htl fishing strategies, all indicators except B_L2H decrease with increasing fishing mortality for most ecosystems under the F_all fishing strategy. Therefore, we conclude that when fishing expands over time on all-trophic-level taxa, the direction of change in ecological indicators is more predictable, with a generally decreasing trend in response to increasing fishing levels emerging. Responses of indicators to fishing mortality: linearity/non-linearity and thresholds Identifying non-linearity is a premise for applying early warning signals to ecosystem-based fisheries management that can be used in a similar way to limit reference points in the precautionary approach (Gabriel and Mace, 1999). However, formal tests of non-linearity are often hampered by typically short and autocorrelated time-series in empirical studies (Litzow and Hunsicker, 2016). Most studies on non-linearity have been based on the identification of threshold responses on GAMs curves fitted to empirical relationships between time-series of response (e.g. indicators) and drivers (e.g. environmental variables and fishing pressure; Large et al., 2013; Samhouri et al., 2017; Tam et al., 2017). The present study explores non-linearity in a literal sense because it deals with independent data points from ecosystem simulations for which autocorrelation is not a concern. Our results suggest that most indicators tend to respond to fishing mortality in a linear way, particularly for the community and biomass-based indicators (e.g. TLco, B_all, B_htl, B_ltl; Figure 8). Although the indicator B/C predominantly responds to fishing mortality in a non-linear way, the majority of trends in the B/C indicator is non-linear-monotonic under the F_all fishing strategy with distinct response direction and straightforward interpretation of response to fishing mortality. This is in line with the suggestion by Litzow and Hunsicker (2016) that most observed ecosystem changes may be parsimoniously explained by linear and reversible tracking of perturbation. For potentially non-linear relationships between response and drivers, GAMs and/or the gradient forest method have often been employed to quantify thresholds (e.g. Large et al., 2013, 2015; Samhouri et al., 2017; Tam et al., 2017; Fu et al., 2019). The previous studies of quantifying thresholds of ecological indicators in their responses to fishing pressure and environmental change have contributed important knowledge to identify operational reference points for moving towards ecosystem-based fisheries management by avoiding regime shifts caused by climate change and/or overfishing (Tam et al., 2017). Although the previous studies identified all inflection points (second derivative of the GAMs curve = 0) or significant shift in the gradient forest analysis as thresholds, the present study deliberately focuses only on the cases where fitted GAMs curves are non-linear-mixed and ignores all other cases of linear tracking and monotonic non-linearity. Fu et al. (2019), using the 14 indicators explored here, concluded that under the F_all fishing strategy, all 14 indicators had threshold responses and over 50% of the threshold responses were around 0.6*FMSY for most ecosystems. Compared with Fu et al. (2019), thresholds found in the present study under the F_all fishing strategy are generally similar, indicating general consistency between these two different approaches. However, by ignoring cases of linear tracking and monotonic non-linearity, the present study indicates that the occurrence of non-linearity is less frequent. In addition, thresholds identified using the fifth order polynomial curves are also infrequent for the non-linear-mixed type in most ecosystems under all fishing strategies (Figures 5–7; Supplementary Figures S1–S3). The elimination of threshold responses resulting from potentially linear and reversible tracking of perturbation is a step further towards ecosystem-based fisheries management by helping managers to identify sudden ecological transitions of real ecological concerns, because over-application of non-linearity and thresholds has caused confusion among managers (Litzow and Hunsicker, 2016). Overall, it is becoming clear that tailoring decision-support systems to the unique biological and environmental characteristics of an ecosystem is critical, and, once again, we stress the importance of taking cognizance of the fisheries exploitation history of an ecosystem (e.g. Shannon et al., 2014; Coll et al., 2016; Lockerbie et al., 2017,, 2018; Briton et al., 2019). Concluding remarks Using GAMs combined with linear regression and polynomial models to analyse the outputs from the multi-ecosystem, multi-model simulation experiment conducted under the IndiSeas programme is useful for uncovering indicators’ responses to fishing mortality in direction, non-linearity, and threshold. The present study has led to the following conclusions: (i) responses of indicators to fishing mortality in shape, direction, and threshold depend on the fishing strategies considered; (ii) the majority of the indicator response curves demonstrates negative relationships with fishing mortality, with a few exceptions which depend on fishing strategies, suggesting the importance of considering the history of fisheries exploitation when interpreting indicators for a particular marine ecosystem; (iii) most indicators tend to respond to fishing mortality in a linear way, particularly for community and biomass-based indicators, indicating that responses of these indicators to fishing mortality are more predictable; (iv) the infrequent occurrence of threshold for the non-linear-mixed response in most ecosystems may suggest that our approach has helped eliminate threshold responses of linear and reversible tracking of perturbation and focus attention only on thresholds of real ecological concerns. Although future simulations should be carried out over a broader range of fishing mortality rate to further explore ecological thresholds, the various responses (linearity, non-linearity, and threshold) demonstrated here should be used to start to define management targets (levels and/or directions) and recommendations for sustainable levels of fishing pressure. Supplementary data Supplementary material is available at the ICESJMS online version of the manuscript. Acknowledgements This is a contribution to the IndiSeas Working Group, co-funded by IOC-UNESCO (www.ioc-unesco.org) and EuroMarine (http://www.euromarinenetwork.eu), to the project EMIBIOS (FRB, contract no. APP-SCEN-2010-II) and to the IOC-UNESCO GOOS Programme and GOOS Biology and Ecosystems Panel. The work on Canada West Coast ecosystem was sponsored by Fisheries & Oceans Canada under the Aquatic Climate Change and Adaptation Services Programme. YJS and LV were supported by the project EMIBIOS (FRB, contract no. APP-SCEN-2010-II). AG was supported by NOAA’s Integrated Ecosystem Assessment (IEA) programme (http://www.noaa.gov/iea/). MC was supported by the Marie Curie Career Integration Grant Fellowships—PCIG10-GA-2011-303534—to the BIOWEB project and the Spanish National Project PELWEB (CTM2017-88939-R). All other authors were supported by their respective affiliations. We thank the participants of the IndiSeas Working Group for the useful discussions and insights since the early stages of the study. We are grateful to Huizhu Liu and Philippe Verley for their technical support on OSMOSE, Penny Johnson for her modelling work on Atlantis, Jeroen Steenbeek for his help with EwE, and Jennifer Houle for the size spectrum model. We are also thankful to Yu Bai for reviewing an earlier version of this manuscript and the R codes, as well as providing suggestions on threshold detection. At last but not the least, we would like to thank Editor Jason Link and three anonymous reviewers whose valuable comments have greatly improved the manuscript. References Akoglu E. 2013 . Nonlinear dynamics of the Black Sea ecosystem and its response to anthropogenic and climate variations. PhD thesis, Middle East Technical University, Ankara, Turkey. Akoglu E. , Salihoglu B., Libralato S., Oguz T., Solidoro C. 2014 . An indicator-based evaluation of Black Sea food web dynamics during 1960–2000 . Journal of Marine Systems , 134 : 113 – 125 . Google Scholar Crossref Search ADS WorldCat Alexander K. A. , Heymans J. J., MaGill S., Tomczak M., Holmes S., Wilding T. A. 2015 . Investigating the recent decline in gadoid stocks in the west of Scotland shelf ecosystem using a food-web model . ICES Journal of Marine Science , 72 : 436 – 449 . Google Scholar Crossref Search ADS WorldCat Andersen T. , Carstensen J., Hernandez G. E., Duarte C. M. 2009 . Ecological regime shifts: approaches to identification . Trends in Ecology and Evolution , 24 : 49 – 57 . Google Scholar Crossref Search ADS PubMed WorldCat Araújo J. N. , Bundy A. 2011 . Description of the ecosystem models of the Bay of Fundy, Western Scotian Shelf and NAFO Division 4X . Canadian Technical Report of Fisheries and Aquatic Sciences , 2952: xii + 189 p. Google Scholar OpenURL Placeholder Text WorldCat Araújo J. N. , Bundy A. 2012 . The relative importance of climate change, exploitation and trophodynamic control in determining ecosystem dynamics on the western Scotian Shelf . Canada. Marine Ecology Progress Series , 464 : 51 – 67 . Google Scholar Crossref Search ADS WorldCat Baudron A. R. , Serpetti N., Fallon N. G., Heymans J. J., Fernandes P. G. 2019 . Can the common fisheries policy achieve good environmental status in exploited ecosystems: the west of Scotland demersal fisheries example . Fisheries Research , 211 : 217 – 230 . Google Scholar Crossref Search ADS WorldCat Blamey L. K. , Howard J. A. E., Agenbag J., Jarre A. 2012 . Regime-shifts in the southern Benguela shelf and inshore region . Progress in Oceanography , 106 : 80 – 95 . Google Scholar Crossref Search ADS WorldCat Blanchard J. L. , Andersen K. H., Scott F., Hintzen N. T., Piet G., Jennings S. 2014 . Evaluating targets and trade-offs among fisheries and conservation objectives using a multispecies size spectrum model . Journal of Applied Ecology , 51 : 612 – 622 . Google Scholar Crossref Search ADS WorldCat Blanchard J. L. , Coll M., Trenkel V. M., Vergnon R., Yemane D., Jouffre D., Link J. S., Shin Y. J. 2010 . Trend analysis of indicators: a comparison of recent changes in the status of marine ecosystems around the world . ICES Journal of Marine Science , 67 : 732 – 744 . Google Scholar Crossref Search ADS WorldCat Boldt J. L. , Martone R., Samhouri J., Perry R. I., Itoh S., Chung I. K., Takahashi M. et al. 2014 . Developing ecosystem indicators for responses to multiple stressors . Oceanography , 27 : 116 – 133 . Google Scholar Crossref Search ADS WorldCat Boyce D. G. , Dowd M., Lewis M. R., Worm B. 2014 . Estimating global chlorophyll changes over the past century . Progress in Oceanography , 122 : 163 – 173 . Google Scholar Crossref Search ADS WorldCat Branch T. A. 2015 . Fishing impacts on food webs: multiple working hypotheses . Fisheries , 40 : 373 – 375 . Google Scholar Crossref Search ADS WorldCat Briton F. , Shannon L., Barrier N., Verley P., Shin Y. J. 2019 . Reference levels of ecosystem indicators at multispecies maximum sustainable yield . ICES Journal of Marine Science , doi:10.1093/icesjms/fsz104. Google Scholar OpenURL Placeholder Text WorldCat Brown C. J. , Schoeman D. S., Sydeman W. J., Brander K., Buckley L. B., Burrows M., Duarte C. M. et al. 2011 . Quantitative approaches in climate change ecology . Global Change Biology , 17 : 3697 – 3713 . Google Scholar Crossref Search ADS WorldCat Bundy A. , Shannon L. J., Rochet M. J., Neira S., Shin Y. J., Hill L., Aydin K. 2010 . The good (ish), the bad, and the ugly: a tripartite classification of ecosystem trends . ICES Journal of Marine Science , 67 : 745 – 768 . Google Scholar Crossref Search ADS WorldCat Burthe S. J. , Henrys P. A., Mackay E. B., Spears B. M., Campbell R., Carvalho L., Dudley B. et al. 2016 . Do early warning indicators consistently predict nonlinear change in long-term ecological data? Journal of Applied Ecology , 53 : 666 – 676 . Google Scholar Crossref Search ADS WorldCat Capuzzo E. , Lynam C. P., Barry J., Stephens D., Forster R. M., Greenwood N., McQuatters-Gollop A. et al. 2018 . A decline in primary production in the North Sea over twenty-five years, associated with reductions in zooplankton abundance and fish stock recruitment . Global Change Biology , 24 : e352 – e364 . Google Scholar Crossref Search ADS PubMed WorldCat Christensen V. , Walters C. 2004 . Ecopath with Ecosim: methods, capabilities and limitations . Ecological Modelling , 72 : 109 – 139 . Google Scholar Crossref Search ADS WorldCat Coll M. , Navarro J., Palomera I. 2013 . Ecological role of the endemic Starry ray Raja asterias in the NW Mediterranean Sea and management options for its conservation . Biological Conservation , 157 : 108 – 120 . Google Scholar Crossref Search ADS WorldCat Coll M. , Palomera I., Tudela S., Sardà F. 2006 . Trophic flows, ecosystem structure and fishing impacts in the south Catalan Sea, northwestern Mediterranean . Journal of Marine Systems , 59 : 63 – 96 . Google Scholar Crossref Search ADS WorldCat Coll M. , Shannon L. J., Kleisner K. M., Juan Jordà M. J., Bundy A., Akoglu G., Banaru D. et al. 2016 . Ecological indicators to capture the effects of fishing on biodiversity and conservation status of marine ecosystems . Ecological Indicators , 60 : 947 – 962 . Google Scholar Crossref Search ADS WorldCat Cury P. , Bakun A., Crawford R. J., Jarre A., Quinones R. A., Shannon L. J., Verheye H. M. 2000 . Small pelagics in upwelling systems: patterns of interaction and structural changes in wasp-waist ecosystems . ICES Journal of Marine Science , 57 : 603 – 618 . Google Scholar Crossref Search ADS WorldCat Fernandes P. G. , Ralph G. M., Nieto A., Criado M. G., Vasilakopoulos P., Maravelias C. D., Cook R. M. et al. 2017 . Coherent assessments of Europe’s marine fishes show regional divergence and megafauna loss . Nature Ecology & Evolution , 1 : 0170. Google Scholar Crossref Search ADS WorldCat Foley M. M. , Martone R. G., Fox M. D., Kappel C. V., Mease L. A., Erickson A. L., Halpern B. S. et al. 2015 . Using ecological thresholds to inform resource management: current options and future possibilities . Frontiers in Marine Science , 2 : 95. Google Scholar Crossref Search ADS WorldCat Frank K. T. , Petrie B., Shackell N. L. 2007 . The ups and downs of trophic control in continental shelf ecosystems . Trends in Ecology & Evolution , 22 : 236 – 242 . Google Scholar Crossref Search ADS PubMed WorldCat Fu C. , Gaichas S., Link J. S., Bundy A., Boldt J. L., Cook A. M., Gamble R. et al. 2012 . Relative importance of fishing, trophodynamic and environmental drivers in a series of marine ecosystems . Marine Ecology Progress Series , 459 : 169 – 184 . Google Scholar Crossref Search ADS WorldCat Fu C. , Large S., Knight B., Richardson A., Bundy A., Reygondeau G., Boldt J. 2015 . Relationships among fisheries exploitation, environmental conditions, and ecological indicators across a series of marine ecosystems . Journal of Marine Systems , 148 : 101 – 111 . Google Scholar Crossref Search ADS WorldCat Fu C. , Perry R. I., Shin Y-J., Schweigert J., Liu H. 2013 . An ecosystem modelling framework for incorporating climate regime shifts into fisheries management . Progress in Oceanography , 115 : 53 – 64 . Google Scholar Crossref Search ADS WorldCat Fu C. , Travers-Trolet M., Velez L., Grüss A., Bundy A., Shannon L. J., Fulton E. A. et al. 2018 . Risky business: the combined effects of fishing and changes in primary productivity on fish communities . Ecological Modelling , 368 : 265 – 276 . Google Scholar Crossref Search ADS WorldCat Fu C. , Xu Y., Bundy A., Grüss A., Coll M., Heymans J. J., Fulton E. A. et al. 2019 . Making ecological indicators management ready: assessing their ability to detect impacts of fishing and environmental change . Ecological Indicators , 105 : 16 – 28 . Google Scholar Crossref Search ADS WorldCat Fulton E. A. , Parslow J. S., Smith A. D., Johnson C. R. 2004 . Biogeochemical marine ecosystem models II: the effect of physiological detail on model performance . Ecological Modelling , 173 : 371 – 406 . Google Scholar Crossref Search ADS WorldCat Fulton E. A. , Smith A. D. M., Smith D. C., Johnson P. 2014 . An integrated approach is needed for ecosystem based fisheries management: insights from ecosystem-level management strategy evaluation . PLoS One , 9 : e84242. Google Scholar Crossref Search ADS PubMed WorldCat Gabriel W. L. , Mace P. M. 1999 . A review of biological reference points in the context of the precautionary approach. NOAA Technical Memorandum NMFS‐F/SPO‐40. pp. 34–45. Groffman P. M. , Baron J. S., Blett T., Gold A. J., Goodman I., Gunderson L. H., Levinson B. M. et al. 2006 . Ecological thresholds: the key to successful environmental management or an important concept with no practical application? Ecosystems , 9 : 1 – 13 . Google Scholar Crossref Search ADS WorldCat Grüss A. , Harford W. J., Schirripa M. J., Velez L., Sagarese S. R., Shin Y-J., Verley P. 2016b . Management strategy evaluation using the individual-based, multispecies modeling approach OSMOSE . Ecological Modelling , 340 : 86 – 105 . Google Scholar Crossref Search ADS WorldCat Grüss A. , Schirripa M. J., Chagaris D., Velez L., Shin Y-J., Verley P., Oliveros-Ramos R. 2016a . Evaluating natural mortality rates and simulating fishing scenarios for Gulf of Mexico red grouper (Epinephelus morio) using the ecosystem model OSMOSE-WFS . Journal of Marine Systems , 154 : 264 – 279 . Google Scholar Crossref Search ADS WorldCat Halouani G. , Lasram F. B. R., Shin Y-J., Velez L., Verley P., Hattab T., Oliveros-Ramos R. et al. 2016 . Modelling food web structure using an end-to-end approach in the coastal ecosystem of the Gulf of Gabes (Tunisia) . Ecological Modelling , 339 : 45 – 57 . Google Scholar Crossref Search ADS WorldCat Hastie T. , Tibshirani R. 1990 . Generalized Additive Models . Chapman and Hall/CRC Press , London, UK . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Heymans J. J. , Tomczak M. T. 2016 . Regime shifts in the northern Benguela ecosystem: challenges for management, in: Ecopath 30 years—modelling ecosystem dynamics: beyond boundaries with EwE . Ecological Modelling , 331 : 151 – 159 . Google Scholar Crossref Search ADS WorldCat Hunt G. , Coyle K. O., Eisner L. B., Farley E. V., Heintz R. A., Mueter F., Napp J. M. et al. 2011 . Climate impacts on eastern Bering Sea foodwebs: a synthesis of new data and an assessment of the oscillating control hypothesis . ICES Journal of Marine Science , 68 : 1230 – 1243 . Google Scholar Crossref Search ADS WorldCat Jackson J. B. C. , Kirby M. X., Berger W. H., Bjorndal K. A., Botsford L. W., Bourque B. J., Bradbury R. H. et al. 2001 . Historical overfishing and the recent collapse of coastal ecosystems . Science , 293 : 629 – 638 . Google Scholar Crossref Search ADS PubMed WorldCat Kelly R. P. , Erickson A. L., Mease L. A., Battista W., Kittinger J. N., Fujita R. 2015 . Embracing thresholds for better environmental management . Philosophical Transactions of the Royal Society B: Biological Sciences , 370 : 20130276. Google Scholar OpenURL Placeholder Text WorldCat Kirby R. , Beaugrand G., Lindley J. 2009 . Synergistic effects of climate and fishing in a marine ecosystem . Ecosystems , 12 : 548 – 561 . Google Scholar Crossref Search ADS WorldCat Large S. I. , Fay G., Friedland K. D., Link J. S. 2013 . Defining trends and thresholds in responses of ecological indicators to fishing and environmental pressures . ICES Journal of Marine Science , 70 : 755 – 767 . Google Scholar Crossref Search ADS WorldCat Large S. I. , Fay G., Friedland K. D., Link J. S. 2015 . Critical points in ecosystem responses to fishing and environmental pressures . Marine Ecology Progress Series , 521 : 1 – 17 . Google Scholar Crossref Search ADS WorldCat Link J. S. , Yemane D., Shannon L. J., Coll M., Shin Y-J., Hill L., Borges M. F. 2010 . Relating marine ecosystem indicators to fishing and environmental drivers: an elucidation of contrasting responses . ICES Journal of Marine Science , 67 : 787 – 795 . Google Scholar Crossref Search ADS WorldCat Lintz H. E. , McCune B., Gray A. N., McCulloh K. A. 2011 . Quantifying ecological thresholds from response surfaces . Ecological Modelling , 222 : 427 – 436 . Google Scholar Crossref Search ADS WorldCat Litzow M. A. , Hunsicker M. E. 2016 . Early warning signals, nonlinearity, and signs of hysteresis in real ecosystems . Ecosphere , 7 : e01614-10 . Google Scholar Crossref Search ADS WorldCat Lockerbie E. , Shannon L. J., Jarre A. 2016 . The use of ecological, fishing and environmental indicators in support of decision making in southern Benguela fisheries . Ecological Modelling , 69 : 473 – 487 . Google Scholar OpenURL Placeholder Text WorldCat Lockerbie E. M. , Coll M., Shannon L. J., Jarre A. 2017 . The use of indicators for decision support in northwestern Mediterranean Sea fisheries . Journal of Marine Systems , 174 : 64 – 77 . Google Scholar Crossref Search ADS WorldCat Lockerbie E. M. , Lynam C. P., Shannon L. J., Jarre A. 2018 . Applying a decision tree framework in support of an ecosystem approach to fisheries: IndiSeas indicators in the North Sea . ICES Journal of Marine Science , 75 : 1009 – 1020 . Google Scholar Crossref Search ADS WorldCat Möllmann C. , Diekmann R., Müller-Karulis B., Kornilovs G., Plikshs M., Axe P. 2009 . Reorganization of a large marine ecosystem due to atmospheric and anthropogenic pressure: a discontinuous regime shift in the central Baltic Sea . Global Change Biology , 15 : 1377 – 1393 . Google Scholar Crossref Search ADS WorldCat Moore J. C. 2018 . Predicting tipping points in complex environmental systems . Proceedings of the National Academy of Sciences of the United States of America , 115 : 635 – 636 . Google Scholar Crossref Search ADS PubMed WorldCat Nicolaou T. , Constandinou G. 2016 . A nonlinear causality estimator based on non-parametric multiplicative regression . Frontiers in Neuroinformatics , 10 : 19. Google Scholar Crossref Search ADS PubMed WorldCat Ortega-Cisneros K. , Cochrane K. L., Fulton E. A., Gorton R., Popova E. 2018 . Evaluating the effects of climate change in the southern Benguela upwelling system using the Atlantis modelling framework . Fisheries Oceanography , 27 : 489 – 503 . Google Scholar Crossref Search ADS WorldCat Pauly D. , Watson R. 2005 . Background and interpretation of the ‘Marine Trophic Index’ as a measure of biodiversity . Philosophical Transactions of the Royal Society B , 360 : 415 – 423 . Google Scholar Crossref Search ADS WorldCat Poloczanska E. S. , Burrows M. T., Brown C. J., Garcia Molinos J., Halpern B. S., Hoegh-Guldberg O., Kappel C. V. et al. 2016 . Responses of marine organisms to climate change across oceans . Frontiers in Marine Science , 3 : 62. Google Scholar Crossref Search ADS WorldCat R Core Team, 2019 . R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing, Vienna, Austria . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Reed J. , Shannon L., Velez L., Akoglu E., Bundy A., Coll M., Fu C. et al. 2017 . Ecosystem indicators—accounting for variability in species' trophic level . ICES Journal of Marine Science , 74 : 158 – 169 . Google Scholar Crossref Search ADS WorldCat Rice J. C. , Rochet M-J. 2005 . A framework for selecting a suite of indicators for fisheries management . ICES Journal of Marine Science , 62 : 516 – 527 . Google Scholar Crossref Search ADS WorldCat Samhouri J. F. , Andrews K. S., Fay G., Harvey C. J., Hazen E. L., Hennessey S. M., Holsman K. et al. 2017 . Defining ecosystem thresholds for human activities and environmental pressures in the California Current . Ecosphere , 8 : e01860. Google Scholar Crossref Search ADS WorldCat Scheffer M. , Carpenter S. R. 2003 . Catastrophic regime shifts in ecosystems: linking theory to observation . Trends in Ecology and Evolution , 18 : 648 – 656 . Google Scholar Crossref Search ADS WorldCat Scheffer M. , Carpenter S., Foley J. A., Folke C., Walker B. 2001 . Catastrophic shifts in ecosystems . Nature , 413 : 591 – 596 . Google Scholar Crossref Search ADS PubMed WorldCat Serpetti N. , Baudron A. R., Burrows M. T., Payne B. L., Helaouët P., Fernandes P. G., Heymans J. J. 2017 . Impact of ocean warming on sustainable fisheries management informs the ecosystem approach to fisheries . Scientific Reports , 7 : 13438. Google Scholar Crossref Search ADS PubMed WorldCat Shannon L. , Coll M., Bundy A., Gascuel D., Heymans J. J., Kleisner K., Lynam C. P. et al. 2014 . Trophic level-based indicators to track fishing impacts across marine ecosystems . Marine Ecology Progress Series , 512 : 115 – 140 . Google Scholar Crossref Search ADS WorldCat Shannon L. J. , Christensen V., Walters C. 2004 . Modelling stock dynamics in the southern Benguela ecosystem for the period 1978–2002 . African Journal of Marine Science , 26 : 179 – 196 . Google Scholar Crossref Search ADS WorldCat Shannon L. J. , Coll M., Yemane D., Jouffre D., Neira S., Bertrand A., Diaz E. et al. 2010 . Comparing data-based indicators across upwelling and comparable systems for communicating ecosystem states and trends . ICES Journal of Marine Science , 67 : 807 – 832 . Google Scholar Crossref Search ADS WorldCat Shannon L. J. , Neira S., Taylor M. 2008 . Comparing internal and external drivers in the southern Benguela and the southern and northern Humboldt upwelling ecosystems . African Journal of Marine Science , 30 : 63 – 84 . Google Scholar Crossref Search ADS WorldCat Shin Y-J. , Bundy A., Shannon L. J., Blanchard J. L., Chuenpagdee R., Coll M., Knight B. et al. ; the IndiSeas Working Group 2012 . Global in scope and regionally rich: an IndiSeas workshop helps shape the future of marine ecosystem indicators . Reviews in Fish Biology and Fisheries , 22 : 835 – 845 . Google Scholar Crossref Search ADS WorldCat Shin Y-J. , Cury P. 2004 . Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing . Canadian Journal of Fisheries and Aquatic Sciences , 61 : 414 – 431 . Google Scholar Crossref Search ADS WorldCat Shin Y-J. , Houle J. E., Akoglu E., Blanchard J., Bundy A., Coll M., Demarcq H. et al. 2018 . The specificity of marine ecological indicators to fishing in the face of environmental change: a multi-model evaluation . Ecological Indicators , 89 : 317 – 326 . Google Scholar Crossref Search ADS WorldCat Shin Y-J. , Shannon L. J., Bundy A., Coll M., Aydin K., Bez N., Blanchard J. L. et al. 2010 . Using indicators for evaluating, comparing, and communicating the ecological status of exploited marine ecosystems. 2. Setting the scene . ICES Journal of Marine Science , 67 : 692 – 716 . Google Scholar Crossref Search ADS WorldCat Tam J. C. , Link J. S., Large S. I., Andrews K., Friedland K. D., Gove J., Hazen E. et al. 2017 . Comparing apples to oranges: common trends and thresholds in anthropogenic and environmental pressures across multiple marine ecosystems . Frontier in Marine Science , 4 : 282 . Google Scholar Crossref Search ADS WorldCat Wood S. N. 2017 . Generalized Additive Models: An Introduction with R , 2nd edn. Chapman and Hall/CRC Press , London, UK . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Wood S. N. , Pya N., Säfken B. 2016 . Smoothing parameter and model selection for general smooth models (with discussion) . Journal of the American Statistical Association , 111 : 1548 – 1575 . Google Scholar Crossref Search ADS WorldCat © International Council for the Exploration of the Sea 2019. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © International Council for the Exploration of the Sea 2019.