Spatiotemporal changes in the Atlantic cod (Gadus morhua) stock at Flemish Cap (1993–2019) and their relationships with demersal communities

Krerkkrai Songin^*1, Graham Pierce² and Fran Saborido-Rey³

Instituto de Investigaciones Marinas (IIM), CSIC, Vigo, Spain
²Email: g.j.pierce@iim.csic.es
³Email: fran@iim.csic.es

^*1Corresponding author: Krerkkrai Songin, tel: +34603370660; e-mail: ksongin@iim.csic.es

 

Songin, K., Pierce, G., and Saborido-Rey, F. 2024. Spatiotemporal changes in the Atlantic cod (Gadus morhua) stock at Flemish Cap (1993–2019) and their relationships with demersal communities. J. Northw. Atl. Fish. Sci., 55: 59–78. https://doi.org/10.2960/J.v55.m748

Abstract

Changes in Atlantic cod (Gadus morhua) abundance at Flemish Cap, likely due to exploitation and perhaps also to changing environmental conditions, have been well documented since 1980s. While the ecological implications of cod fluctuations have been explored in relation to dominant and commercially important species including redfish (Sebastes spp.), northern shrimp (Pandalus borealis) and Greenland halibut (Reinhardtius hippoglossoides), the broader ecological impacts, e.g. on less abundant species, remain less well explored. This study aimed to analyse spatiotemporal variation in the distribution and abundance of cod, and identify associated changes in distribution and abundance of other species with various trophic relationships to cod. This analysis used a delta Generalized Additive Model (GAM) approach, incorporating binomial and quasi-Poisson GAMs fitted to EU bottom trawl survey data from 1993 to 2019. Trophic species and guilds were defined based on the sizes and feeding habits of each species, as established in previous studies. Atlantic cod is considered to comprise of two trophic species: cod under 46 cm and larger cod. Model predictions were used to construct distribution maps and estimate distribution range and annual total abundance. Bottom temperature was a more important predictor in abundance (quasi-Poisson) models than in presence (binomial) models. The observed decline in cod abundance was associated with contraction in the distribution range. Significant negative correlations were identified between cod trophic species and all but one of the other trophic species in the same trophic guilds, for both distribution range and abundance. Species in other trophic guilds that rely on northern shrimp as prey also exhibited negative correlations with cod. The abundances of the main prey of cod, namely juvenile redfish and northern shrimp, showed negative correlations with cod abundance but no relationship was seen for distribution range. The abundance of large Acadian redfish (S. fasciatus) and large beaked redfish (S. mentella), which are major prey species of cod, was positively correlated with that of large cod, suggesting that the abundance of these prey species depends more on external variables, such as intense exploitation, than on their predator-prey relationships. These findings highlight the importance, for fishery management, of considering both the direct effects of fishing mortality and the indirect effects via trophic relationships.

 

Keywords: demersal fish, deep-sea, ecology, trophic relationship, Northwest Atlantic

PDF, Supplementary Materials

 

---

Introduction

The collapse of Atlantic cod (Gadus morhua) abundance in the Northwest Atlantic in the late 20th Century has been well documented (Boudreau et al., 2017) but less attention has been given to associated changes in the spatial dimension (Pérez-Rodríguez et al., 2012, 2017; Garrido et al., 2023a) or to wider ecosystem effects (Dawe et al., 2012).

As the abundance of a population declines, ultimately its distribution range will also decline (Rose et al., 2000; Thorson et al., 2016). Identifying a fished population's spatial response to external stress can provide valuable information for spatial management measures, such as protected areas, which could ensure preservation of the stock when depleted. Fishery closures, even if temporary, can be employed as precautionary measures to maintain a healthy stock condition (Gell and Roberts, 2003). The Northwest Atlantic Fisheries Organization (NAFO) has implemented seasonal area closure for cod fisheries in Division 3M during the spawning period, i.e. in the first quarter of the year (NAFO, 2021). Permanent closure has also been employed in parts of NAFO Division 3M to conserve the biodiversity of vulnerable ecosystems (NAFO, 2022). Better understanding of the relationships between fish abundance and distribution and their changes under various circumstances could help improve spatial management.

In a complex ecosystem, fluctuations of key species such as cod, whether caused by human activities or varying environmental conditions can, in turn, influence their competitors and prey (Pérez-Rodríguez et al., 2017). Areas that have experienced overfishing of high trophic level species often show regime shifts and trophic cascades that are challenging to reverse and may even be irreversible (Fisher et al., 2015). The nature of such changes is, however, not necessarily obvious. Thus, the decline of cod, a higher trophic level species, would lead to both reduced pressure on its prey and reduced competition for other carnivores, effects that would ultimately be antagonistic.

Flemish Cap is an underwater mountain with a summit at a depth of 125 m, located near the Canadian Exclusive Economic Zone (EEZ) to the east of Grand Bank. It is one of the historically most significant fishing grounds for cod. The area is regulated by NAFO as part of management area Division 3M. It is semi-isolated from the continental shelf by a 1 100 m deep channel known as Flemish Pass, which limits migration between the seamount and the continental shelf (Stein, 2005; Vázquez et al., 2014). Records of cod exploitation and biomass in the area clearly showed signs of overfishing and stock collapse in the late 1990s (González-Troncosoet al., 2022). A fishery moratorium was implemented in 1999 but the overall abundance of cod continued to decline and did not show clear signs of recovery until the late 2000s, suggesting the involvement of unfavourable environmental conditions in delaying the recovery (Ruiz-Díaz et al., 2022). The population subsequently recovered sufficiently to reopen the fishery in 2010 and has since been closely monitored to maintain spawning stock biomass (SSB) above the reference limit set by NAFO Scientific Council (NAFO, 2019).

There is evidence that the wide fluctuations in cod abundance over the years have had significant impacts on the ecosystem. Flemish Cap hosts a diverse range of over one hundred identified fish species, most being demersal (Vázquez et al., 2013). Since 1988, the EU has conducted fisheries-independent surveys in the area, providing data that have been used to investigate the impact of fishing and the environment on fish stocks and community structure (Pérez-Rodríguez et al., 2012). The complex trophic interactions in the area have been studied for common fish species, revealing shared feeding habits such that demersal fish can be grouped into four trophic guilds (Pérez-Rodríguezet al., 2011). Guild I members predominantly prey on small invertebrates, mainly northern shrimp (Pandalus borealis) and hyperiids. Guild II exhibits the most diverse invertebrate diet and some of the members rely largely on northern shrimp but still have greater dietary variation than Guild I members. Guild III consists exclusively of redfish with diets strongly based on pelagic invertebrates and northern shrimp. Guild IV is predominantly piscivorous, with one of the main prey items being redfish. Recent studies modelling the dominant system of cod-redfish-shrimp, considering their biomass and energy flows in the area, have shown an intensely intertwined relationship among them, and the change in one species due to fishing or natural mortality could profoundly affect the others (Pérez-Rodríguez et al., 2022). In the period of cod depletion, a major assemblage shift was detected, with other normally less abundant species in the area such as Greenland halibut (Reinhardtius hippoglossoides) and wolffishes (Anarhichus spp.) becoming more dominant (Pérez-Rodríguez et al., 2017).

While changes in the demersal fish community have been observed in recent decades, the spatial aspects of these changes have been less well explored (Pérez-Rodríguez et al., 2012; Pham et al., 2019). A previous study exploring changes in fish distribution through spatial modelling approaches identified contraction in the distributions of cod and another commonly co-occurring species in the area, American plaice (Hippoglossoides platessoides), when their biomass is low, and a significant increase in the dominance of Greenland halibut (Hendrickson and Vázquez, 2005). It is expected that the spatial structure of the demersal community as a whole, including species that are less abundant and/or of lesser commercial interest, will also have been influenced by the fluctuations in cod abundance, and these wider ecological impacts of cod fluctuations still need to be investigated. Understanding the effects on the demersal community of fluctuations in key species can help inform ecosystem-based fisheries management (Hilborn, 2011).

This study aims to explore the ecological impacts of cod stock fluctuations by investigating the correlations in distribution and abundance between cod and other trophic species. Specifically, we (i) investigate changes in the spatial distribution and abundance of cod in Flemish Cap over periods of differing fishing intensity during 1993 to 2019 and compare abundance estimates based on the modelling approach with those from traditional calculations, (ii) explore the influence of environmental variables on changes in abundance and distribution, and (iii) examine changes in distribution and abundance of other species with varying trophic associations to cod.

Material and methods

Survey sample collection

Flemish Cap bottom trawl surveys have been conducted by the EU since 1988, following the protocol described by NAFO and led by four collaborating institutes from Spain and Portugal: the Institute of Marine Research (IIM-CSIC), the Spanish Institute of Oceanography (IEO), AZTI-Tecnalia Foundation and the Institute of Fisheries and the Sea (IPMA) (Vázquez et al., 2014). Sampling was conducted using two vessels, RV Cornide de Saavedra from 1988 to 2002, reaching a maximum depth of 730 m, and RV Vizconde de Eza since 2003. The depth limit was extended to 1460 m in 2004 to cover a larger area where species preferring deeper water resided. Due to the vessel replacement, the catch data of the focal species including Atlantic cod, redfish, American plaice, Greenland halibut and roughhead grenadier (Macrourus berglax) prior to 2004 were transformed using calibration factors calculated by González-Troncoso and Casas (2005) to standardize the fishing effort to that of the current vessel. The study area was divided into 34 strata for stratified random sampling (Fig. 1). The initial 19 strata are limited to a depth of 730 m and are consistently present each year. The additional strata were sampled from 2004 onward, except strata 26 and 27, which were irregularly sampled due to a dense sponge ground and were only present from 2004 to 2007. The sampling effort was approximately one haul per 100 square miles per year for every stratum except for stratum 33, the smallest stratum, which received a sampling effort close to two hauls per 100 square miles per year.

A Lofoten trawl was used for sampling due to its durability and similarity to the gear used by commercial fisheries in the area. The cod-end mesh size was 35 mm, which is considered to be appropriate for retaining juveniles of commercial species. The trawl opening was 0.0075 mile wide, and the boat maintained a speed of 3.5 knots for 30 minutes after the net made contact with the bottom. Trawling was conducted between 6:00 and 22:00 h daily.

The captured fish were initially identified to species before individual measurements were taken, with the exception of juvenile redfish (<15 cm fork length) due to practical challenges associated with species identification for this genus. Fish size was generally recorded as total length (TL) in centimetres (to the nearest 1 cm), except for redfish, for which fork length (FL) was recorded, and Macrourid fishes (roughhead grenadier and marlin-spike grenadier (Nezumia bairdii)), for which anal length (AL) to the nearest 0.5 cm was recorded. The size of northern shrimp was measured by carapace length in millimetre. In cases where a particular species was highly abundant in the catch and it was impractical to measure all individuals, subsamples were selected for measurement based on a stratified random procedure. The target subsample size was 50 individuals (25 fish for each identified sex) for each species in each haul, except for juvenile redfish for which the target was 20 individuals. The overarching goal was to collect 10 samples for each class by the conclusion of each survey.

Trophic species

To investigate the ecological impact of fluctuations in Atlantic cod abundance, 14 other fish species and northern shrimp were chosen for the study. All fish species and juveniles of redfish were classified into trophic species based on their size and feeding habits, following a previous study on the main fish species in the area (Pérez-Rodríguez et al., 2011) (Table 1). These authors classified the trophic species into four trophic guilds (I to IV) based on observed dietary similarity between 1993 to 2008. In the present study, we include an additional guild (V) for northern shrimp as it is one of the key prey species in Flemish Cap ecosystem that should be modelled for its distribution and abundance.

Bottom temperature

July monthly average bottom temperatures (BT) spanning from 1993 to 2019 were obtained from Copernicus Marine Service (CMEMS). The monthly average temperature was used instead of the in-situ data to minimise the contribution of daily BT variation within each survey. The consistency of bottom temperature mapping by CMEMS is also advantageous for predictive spatial models. The BT data were based on reanalyses from the CMEMS Global Ocean Reanalysis Products, which incorporate real-time global forecasting (Drévillon et al., 2021). The standard grids used had a spatial resolution of 1/12°. Temperature values were attributed to trawls based on the specific grids corresponding to the sampling locations.

Data analysis

Density calculation

Density, expressed as number of individuals per square mile (n/mile²) was calculated for each sampling location by dividing the total number of individuals caught in the haul by the trawled area (i.e., trawling distance × trawl opening width).

Spatiotemporal GAM

The delta GAM approach, as described by Grüss et al. (2014), was utilized to analyse the spatial structure of abundance for each trophic species over time. This approach combines a binomial GAM that describes the probability of presence and a quasi-Poisson GAM that describes the density when a trophic species is present. The models were fitted using the “mgcv” package (Wood, 2023) in R version 4.1.1 following the equations:

g(η) = s(latitude) +s(longitude) + s(depth) + s(year) + s(BT) (1)

g(η) = s(latitude, longitude) +s(depth) + s(year) + s(BT) (2)

g(η) = s(latitude, longitude, year) + s(depth, year) + s(BT) (3)

where η is the binary response of presence/absence data for the binomial GAM, and the density response when non-zero data are present for quasi-Poisson GAM, g represents the model link function, and s is a regression spline. The spline used for the modelling is the “thin-plate regression spline”, chosen for its flexibility and data driven suitable for when there is no prior knowledge of the form of relationships among the variables. The Restricted Maximum Likelihood (REML) method was used to estimate the smoothness parameters of the smooth functions in GAMs. To determine the best-fitting model for each trophic species, interactions among certain variables were explored by comparing models (1), (2) and (3). A logit link function was used for binomial GAMs, and a log link was used for quasi-Poisson GAMs. Note that coordinates were transformed from the recorded decimal degrees to Universal Transverse Mercator (UTM) coordinate system before being fitted to the models.

The best-fit binomial GAMs were selected for each trophic species using the Akaike Information Criterion (AIC) (Akaike, 2011). Deviance explained was used to identify best-fit quasi-Poisson GAMs since the quasi-Poisson family does not provide the necessary likelihood information to compute the AIC. The adequacy of the basis dimension (k) was checked for all best-fit GAMs using the gam.check function. Delta GAM predictions were generated by multiplying the probability of presence estimated by the best-fit binomial models with the density estimated by the best-fit quasi-Poisson models at the same coordinates. Thus, the delta predicted density (y, n/mile²) was calculated using the equation:

y = p*a (4)

wherep represents probability of presence, and a represents density if present.

Model prediction and validation

Best-fit delta GAMs were validated by performing non-parametric Spearman’s correlation analysis to compare observed and predicted densities. A test with 1 000 bootstraps sampling with replacement was used to estimate the Spearman´s correlation coefficient (ρ) and conclude whether ρ was significantly different from zero. In addition, the magnitude of differences between the observed density and that predicted by delta GAM is observed through value ranges and distributions.

The best-fit delta GAMs were used to predict density by year, coordinates, depth and temperature across the 1/20° gridded maps. The resolution was chosen based on the efficiency in computational effort while providing a fine spatial scale for observing spatial variation in presence probability and density. The depth inputs for the prediction were extracted from the seamless gridded topographic and bathymetric bare-earth evaluation data set ETOPO 2022 (NOAA, 2023). The temperature inputs for the model prediction are from the previously extracted CMEMS data. Following the sampling depth fitted to the model, the predictions from 1993 to 2003 were made for depths up to 730 m while the prediction since 2004 are for depths up to 1 460 m.

The annual total abundances of the focal species, Atlantic cod, as predicted by best-fit delta GAMs, were compared to the conventional calculations based on the same survey data as published by NAFO (González-Troncoso et al., 2022). The predicted total abundance was calculated as followed:

Total abundance = ∑i=0nyi*GAi (5)

where y_iis density at gridiand GA_iis the area of grid square _i.

As the published abundance values are provided by age class, the predicted total abundance of GM1was compared to summed abundances of fish aged 1 to 3 years old, while the abundance of GM2 was compared to the summed abundance of older age classes. The comparison was performed using Wilcoxon signed rank test.

Mapping and distribution area calculation

Chronological (annual) maps were constructed for all trophic species for 1993–2019, following the available temperature map data. Mapping was carried out using QGIS version 3.20. The distribution area of each trophic species was defined as the area with a density of at least 10 individuals per square mile. The centroids of polygons weighted by density were then calculated.

Correlation of cod to other trophic species

Although the delta GAMs are constructed using all available data and the predictions for mapping since 2004 were made at all depths, only the annual distribution area and total abundance at depths of less than 730 m (i.e., the depth ranged sampled prior to 2004) were used for correlation analyses between cod and other trophic species. This is due to the fact that cod does not generally occur at greater depths and its correlations to other species can be expected to be much more discernible if the analysis is confined to depths 730 m or less than if we were to include deeper areas.

The correlations between annual distribution area of cod (≤730 m depth) and those of other species throughout the study period were estimated using Spearman’s correlation test with 1000 bootstraps. The correlations between annual total abundances for ≤730 m depth were also tested using the same method. Due to the potential lagged effects of cod population abundance on other species, the correlations paired cod total abundance or distribution area of a given year to the abundance or distribution area of other species from the following year. The one-year lag was chosen as it is the most immediate possible time to assess the direct impact of cod on other trophic species with the least interference from recruitment.

Results

Best-fit Delta GAMs

A total of 3 941 hauls was included in the binomial GAM fitting. For most trophic species, the best-fit binomial GAMs revealed significant interactions between coordinates and year, as well as between depth and year (p-value <0.05) (Table 2). The only exceptions were MB1, MB2 and HP, as their best-fit binomial GAMs lacked interaction terms. There are 12 trophic species for which the temperature term has a significant effect in the best-fit binomial GAMs.

The number of hauls with non-zero density varied among species, ranging from 876 to 2 649 hauls (Table 2). The best-fit quasi-Poisson model for all trophic species except AD featured significant interaction terms between spatial variables and year (p<0.05). In contrast to binomial models, most of the best-fit quasi-Poisson GAMs included a significant temperature term.

Both binomial and quasi-Poisson models explained over 50% of the deviance for most species, irrespective of their trophic guild. However, trophic species such as AR, AD, AL2 and AM2 were less effectively described by both GAMs. Despite belonging to different guilds, these species share the characteristic of lower abundance (<120 n/mile2)compared to the others while still commonly occurring, as can be seen through the high number of non-zero haul (>1 300 presences) (Table 2). GC is another species that is less well described by the binomial GAM but its density is well-fitted by quasi-Poisson GAM.

Spearman’s correlation tests indicate significant (p<0.01) and positive correlations between all observed and predicted values (Table 3). The majority of the Spearman’s p exceed 0.7 for both presence-absence and density tests. Strong positive correlations were found between observed non-zero density and predicted values even for trophic species for which the best-fit quasi-Poisson GAMs had low values of deviance explained, which suggests that the models were reasonably reliable for prediction. Only observed density and quasi-Poisson GAM predictions of SM1, SM2 and AD show p below 0.5. The correlation test also revealed strong and significant (p<0.01) positive correlations (p>0.7) between the observed density (0 included) and delta GAM prediction for most species.

Due to the large number of absences (0 density) in sampling, which could significantly skew the value distribution to the median of 0, the predicted density was compared to observed density only for hauls with presence (non-zero density). The values of observed non-zero density and the predictions by delta GAM are similar in range and distribution in all trophic species. The detailed comparisons between the observations and predictions can be found in the supplementary materials (Fig. S1). The median of observed values falls within the predicted interquartile range for all trophic species except AD, which also show the lowest p between observations and predictions. The observed and predicted density distribution and range are more similar in species that show stronger correlations between them.

Bottom temperature effects

The bottom temperature at sampling locations during the study period typically fell within a narrow interquartile range of 1.5ºC>, between 3 and 4.5ºC (Fig. 2). However, the minimum-maximum range in within the entire period is approximately 3.3ºC, ranging from 1.9–5.2ºC. The range within one year can be as large as 2.2ºC, as seen in 2005. The spatial variation in bottom temperature displayed a distinct pattern of warmer southwest waters and colder northeast waters in most years (Fig. 3). Nevertheless, there are exceptional years showing a clear departure from this pattern such as seen in 1994, 2000 and 2018.

Trophic species for which there were significant relationships between occurrence and temperature share a common pattern whereby occurrence increases with temperature upw to 3ºC (Fig. 4a). Trajectories at warmer temperatures diverge and can be categorized into three types. The first type shows a continuous increase with low deceleration, exemplified by PC, SF2 and PB. The second type reveals a trajectory that plateaus, as demonstrated by NB1, NB2, SM1 and RH2. The third type features a trajectory that re-accelerates from a plateau when the temperature surpasses a certain threshold (around 5ºC), as seen in AR, SJ, SF1, SN1 and GM2. It is important to note that the models were fitted to data with higher numbers in moderate temperatures and fewer instances for extreme highs and lows. This may impact the fitting accuracy as can be seen in wide 95% confidence intervals.

Fig. 4. (a) Smoothers for the partial effect of bottom temperature on occurrence (binomial GAMs), and (b) smoothers for the partial effect of bottom temperature on density (quasi-Poisson GAMs). Shaded areas represent 95% confidence intervals, with each mark along the inside of the x-axis corresponding to a single observation (the “rugplot”). Note that the scale of the x and y-axes varies among plots for display purposes.

Significant effects of BT on density revealed by quasi-Poisson GAMs exhibited more diverse trends than binomial GAM fits. Only AL1 and GM2 show a continuous increase in density with temperatures warmer than 3ºC (Fig. 4b). In contrast to binomial GAM fits, multiple trophic species show a clear negative effect of warming temperature on density. These include MB1, NB2, HP, AD and AM2. In some trophic species, a distinct peak in density can be identified within commonly observed temperature ranges. These distinct peaks appear in SJ between 4–5ºC, and SN1 along with SN2 between 3–4ºC.

Distribution and abundance prediction

Densities for all trophic species were estimated on 1/20° grids using the corresponding best-fit delta GAMs. Density maps of every second year were generated for GM1 (Fig. 5a) and GM2 (Fig. 5b), illustrating spatiotemporal variations in density. The mapping results for all trophic species in all years can be accessed in the supplementary materials (Fig. S2–S27).

The mapping results for GM1 reveal a declining trend from the initial period (Fig. 5a). The density and distribution area reached their lowest levels in 1999 when the fishing moratorium was implemented. The stock remained stagnant until the mid-2000s when signs of recovery started to appear. Both density and distribution area increased until the moratorium was lifted in 2010. However, the recovered high abundance experienced a subsequent decline after 2012. Notably, during periods of relatively low abundance, there was a contraction in the distribution area with small cod remaining mainly at the top part of the cap. The weight centroids of polygons slightly shifted away from near the top of Flemish Cap towards the deeper parts in the northeast direction when the density was high.

The density and distribution of GM2 generally followed a similar trend to GM1, experiencing the decline in the early years (Fig. 5b). GM2 abundance bottomed in 2001, two years after GM1 reached its lowest level (Fig. 6). The distribution area continued to contract until it reached its smallest coverage in 2004. The stock then recovered and reached its highest abundance in 2014, three years after GM1 achieved its highest level. Signs of decline can again be seen in the later years. Although the trend is similar to GM1, the area coverage is more stable for GM2 with less contraction when overall abundance decline. The changes in the weight centroids of polygons are much less perceivable than in GM1

he annual total abundances of Atlantic cod in Flemish Cap estimated by delta GAMs were validated by comparing them to results from NAFO Division 3M cod abundance estimation. The GAM predictions closely resemble the reported abundances (Fig. 6). Wilcoxon signed-rank tests showed no significant differences in total abundance between NAFO estimated Atlantic cod aged 1–3 years and delta GAM prediction of GM1 (p-value= 0.09) or between NAFO estimated Atlantic cod older than 3 years and delta GAM prediction of GM2 (p-value = 0.79). Although the sign-rank tests indicate that there is no significant difference in the medians between the two estimation methods, the magnitude of differences is notable in GM2. Compared to the NAFO estimation, the result from delta GAM of GM2 is discernibly higher than the NAFO swept-area estimation when the reported abundance is low (late-2000s) while being lower when the reported number is high (mid-2010s).

Cod correlations to other trophic species

Total area of Flemish Cap (≤730 m depth) was estimated to be 10 862.10 miles^2. The occupied portion of this area by each trophic species was estimated for each year (Fig. 7a). The area occupied by cod fluctuated greatly. For GM1, its smallest occupied area was 69% smaller than the maximum coverage. For GM2, the smallest occupied was 42% smaller than the maximum coverage. The occupied areas of most members of Guild III (SJ, SF1, SF2, SM1 and SM2) and Guild V were more stable than others over time.

Spearman’s correlation analysis revealed significant strong positive correlations between the GM1 and GM2 distribution area. The correlation between GM2 and a later year’s GM1 is stronger than the correlation between GM1 and a later year’s GM2 (Table 4). Both size classes displayed significant negative correlations with most of the other trophic species in a following year, regarding both abundance and distribution area (Table 4). Both GM1 and GM2 were significantly negatively correlated to most trophic species within their respective guilds. The only exception is NB2 not being significantly correlated to GM1. Most members of Guild II are also significantly negatively correlated to GM1 and GM2, except for shallow water flatfishes GC and HP, which showed strong positive correlations. Most of Guild III members were not significantly correlated to cod, except for positive correlations between SF2 to both GM1 and GM2; negative correlations of SM1 to both GM1 and GM2; and negative correlation between SN1 and GM2. Guild V was not significantly correlated with either GM1 nor GM2.

The annual total abundance (≤730 m depth) for all species was estimated, and fluctuations over time are illustrated in Fig. 7b. The changes in total abundance for most species mirrored the expansions and contractions of their distribution area, with the exception of redfish in Guild III and northern shrimp in Guild V, for which total abundance fluctuations were much more pronounced than distribution area changes. Spearman’s correlation analysis results for abundance were similar to those for distribution area in most cases (Table 4). A strong, significant positive correlation was found between GM1 and GM2, both of which were consistently negatively correlated with members of Guilds I and IV, except for AD. The cod abundance correlations to Guild II were similar to distribution area correlation but weaker. In contrast to the results for distribution area, both GM1 and GM2 showed a strong, significant negative correlation with PB abundance. A significant negative correlation between SJ and both cod classes was found. Other notable differences in correlation analysis results are the non-significance of SM1 and SN1, while a significant positive correlation was shown between SM2 and GM2. Notably, the correlations for abundance were generally weaker than those for distribution area. Overall, the correlations were all stronger for GM2 than GM1, except the correlation to PB.

Discussion

The present study employed a variety of statistical tools to investigate the spatiotemporal dynamics of Atlantic cod and other demersal species in Flemish Cap. Our findings shed light on the changes in the cod abundance and distribution over time due to fishing pressure and environmental variation, and reveal how the fluctuation of this top predator impacted the entire demersal community. Moreover, the study highlights the heterogeneity of these impacts among different trophic groups.

Delta GAM for spatiotemporal modelling

The estimation of cod abundance using the delta GAM approach in this study yielded results that aligned with observations and traditional methods (González-Troncoso et al., 2022). Nevertheless, it is important not to disregard the noticeable differences in estimations between the two methods, especially evident in large cod. The swept-area method estimates abundance by multiplying the density by the stratum area for each stratum, before combining the abundance from all strata (Vázquez et al., 2014). On the other hand, the GAMs estimate the number of individuals in fine grids. Then the numbers of individuals in all fine grids are combined to derive the total abundance. The use of fine grids by the modelling approach could pick up finer scale changes that the swept-area extrapolation did not. Moreover, the modelling approach derived the estimates from the use of covariates, unlike the relatively simple multiplication by the swept-area method which did not consider spatiotemporal factors. Although the delta GAM approach can offer valuable spatial information on a finer scale, further exploration is required to assess its accuracy and precision when determining abundance, compared to the traditional estimation method, before considering its use for stock assessment. We consider sampling effort bias among strata to be negligible, but spatial treatment or post-stratification schemes that could help reduce model uncertainty in the future are worth exploring. Integrating fishery data in addition to the survey data for more input and validation could also improve the modelling, as estimations by NAFO also use commercial fishery data (Garrido et al., 2023a).

To gain further insights into the drivers of overall abundance and distribution changes, it is essential to integrate environmental factors into the modelling and mapping approach. Previous studies have explored depth-related demersal assemblages in Flemish Cap (Nogueira et al., 2017). However, our best-fit models indicate that other environmental variables could also affect the distribution. Sea bottom temperature was clearly shown to be a significant factor affecting density and distribution. All models also included significant terms for spatial coordinates which suggested the presence of unexplained spatial variation due to spatial variation in other environmental variables that were not included in modelling. Environmental variables such as sediment types and benthic faunas, including deep-sea corals and sponges, which vary spatially, could play a role in habitat preference (Murillo et al., 2011, 2016, 2020b, 2020a). Future analyses that incorporate these factors could further explain the assemblage variation found in this study. In this study we used year as a fixed factor to arbitrarily represent possible temporally dynamic factors such as fishing mortality, currents, salinity, primary productivity and zooplankton. In the future investigation when these factors are explicitly included in the model, a modification of GAMs to generalized additive mixed models (GAMMs) to treat year as a random effect might be needed (Lin and Zhang, 1999).

While the spatiotemporal variation in bottom temperature is typically low, as expected in deep-sea environments (Yasuhara and Danovaro, 2016), incorporating this variable into the models has proven important for predicting occurrence and abundance of several species. This supports the hypothesis that changes in the distribution patterns of Flemish Cap are associated with oceanographic conditions(Cerviño et al., 2005). Other investigation also found abiotic factors such as the North Atlantic Oscillation, which is linked to local oceanographic variables including temperature, salinity and currents, to be associated with changes and trends in Flemish Cap demersal community (Pérez-Rodríguez et al., 2012). Binomial GAMs with a significant temperature term unanimously show bottom temperature below 3ºC to be less favourable for fish occurrence but varying effects of temperature can be seen on density. In most of the trophic species for which there appeared to be a negative relationship between temperature and density (when the species was present), including small roughhead grenadier, American plaice, northern wolffish (A. denticulatus) and large spotted wolffish (A. minor, temperature apparently did not affect occurrence. The lack of significant temperature term in binomial GAMs of several species could be because the temperature effects are relatively weak and/or because such effects are partially confounded with depth or spatial variables.

Cod abundance and distribution fluctuation

The findings of this study reinforced the notion that a significant decline in the overall population of a species could come with a contraction of its distribution range, which had previously been shown in multiple taxa across multiple areas (Thorson et al., 2016; Orio et al., 2019). In the case of Flemish Cap cod, abundance and distribution can be observed across three distinct periods: the pre-moratorium, the moratorium and the post-moratorium. In the pre-moratorium period, the abundance and distribution area of cod markedly declined, evidently due to high fishing mortality pressure on low SSB until 1999 when the moratorium was enforced (González-Troncoso et al., 2022).

Throughout the moratorium, both abundance and distribution area further declined until the mid-2000s, marking the species' lowest recorded abundance despite the cessation of fishing pressure. Unfavourable environmental conditions and limited food availability likely hampered recruitment and delayed recovery (Ruiz-Díaz et al., 2022). Notably, bottom temperatures also declined from 1999 to their lowest point in 2006. This potentially affected recovery as the models indicated that cod generally thrives in warmer temperatures. During this low abundance period, the distribution area was predominantly confined to the top of the cap at depths shallower than 270 m. This aligns with a previous spatial modelling study describing cod distribution changes from 1988 to 2002 (Hendrickson and Vázquez, 2005). A noticeable distribution expansion occurred only in the late-2000s when species’ density showed signs of recovery, coinciding with the increase of temperature since 2006.

The recovery from the mid to late-2000s prompted the lifting of the moratorium in 2010 and fishing mortality was reintroduced (Garridoet al., 2023a). In the post-moratorium period, both the abundance and distribution area began to decline once more but remained within the “safe” zone, maintaining a biomass above the reference point and fishing mortality below the target(NAFO, 2023b). After the abundance peaked in 2011 for small cod and 2014 for large cod, the distribution contracted again, following the decline in stock status. This pattern suggests that the most preferred habitat for cod is in the shallowest part of the mount, and they only occupy deeper space when the area becomes sufficiently congested, likely reflecting intraspecific competition (Planque et al., 2011; Thorson et al., 2016).

Observing changes in spatial distribution by size classes is advantageous, as different size classes do not play the same ecological role and are subject to different fishing pressure (Kindsvater and Palkovacs, 2017). This study reveals a clear decline in larger cod density that started several years prior to the decline of small cod. This is likely due to the higher catchability of older and larger fish, which were intensively exploited and depleted before the impact escalated to cause recruitment overfishing (Walters and Maguire, 1996; González-Troncoso et al., 2022). Conversely, the signs of recovery appeared earlier in small cod, as expected, preceding the recovery of the larger size class in subsequent years.

Cod fluctuation effects on the community

The spatiotemporal changes observed in other species, many of which are less threatened by fisheries, could be caused by changes in the trophic structure rather than direct fishing impacts (Nogueira et al., 2016). The severe depletion of cod, particularly the larger size class that predominantly preys on fish, notably juvenile redfish, between the mid-1990s and mid-2000s, might have opened up more access to resources for other piscivores of Guild IV (Pérez-Rodríguez et al., 2017). A prominent shift in distribution was evident in large Greenland halibut, which typically resides in deeper areas(Nogueira et al., 2017). Its presence and density drastically increased in shallower regions when large cod populations collapsed, causing significant changes in the assemblage structure. This finding aligns with a previous spatial modelling study, which clearly illustrated the contraction of the distributions of cod and American plaice from 1988 to 2002, and their replacement by Greenland halibut (Hendrickson and Vázquez, 2005). A multispecies model also describes Flemish Cap ecological shifts from a system in which the dominant species by biomass were cod and redfish to a system characterised by the predominance of other fish species by the late-1990s, before the subsequent recovery (Pérez-Rodríguez et al., 2017). This study's wider scope sheds light on other species that share trophic niches with cod, large and small, which showed increases in abundance and expansion in their distribution range in the absence of Atlantic cod. These species include thorny skate (A. radiata), roughhead grenadier, marlin-spike grenadier, longfin hake (P. chesteri), Arctic eelpout (L. reticulatus), Greenland halibut and wolffishes.

While small cod primarily feed on northern shrimp and hyperiids, their increased presence can still have adverse effects on other fish species due to competition. The spatial distribution of northern shrimp remained relatively stable over shallow areas, yet its overall abundance fluctuated greatly, showing a strong negative correlation to small cod. As the cod population recovered to the point of lifting the moratorium, the shrimp fishery had to enter its own moratorium in 2011 due to the species being unable to sustain both the fishery and increased predation mortality (Casas Sánchez, 2023). Other trophic species in Guild I, which largely consume northern shrimp, exhibited strong negative correlations with small cod in terms of both distribution area and abundance (Pérez-Rodríguez et al., 2011). Multiple species in Guild II, known for having the most diverse diets, also demonstrated negative correlations with small cod due to northern shrimp being a significant part of their diet. In contrast, benthic fishes in Guild II that rely less on northern shrimp, namely American plaice which primarily preys on ophiuroids and witch flounder (G. cynoglossus), which mainly feeds on polychaetes (Link et al., 2002; Gonzálezet al., 2005), appeared to benefit from the increase in the cod population. These species could be less threatened by cod competition and predation but more susceptible to other predators, such as the large Greenland halibut, a benthic flatfish known to prey on other flatfish (Hovde et al., 2002). It should also be noted that some of the correlation results might arise from the species reacting to the changes in bottom temperature or other unaccounted variables. The densities of some species that show strong negative correlations to cod such as small roughhead grenadier, large marlin-spike grenadier and large spotted wolffish also show an opposite relationship with temperature compared to that seen in cod.

While the total abundance of northern shrimp and juvenile redfish is negatively correlated with cod, reflecting their established prey-predator relationship(Pérez-Rodríguez and Saborido-Rey, 2012; Pérez-Rodríguez et al., 2017), their distribution areas do not exhibit a clear correlation. The distribution area of northern shrimp and redfish, including juveniles, remained consistently high throughout the study period, irrespective of the abundance situation. This suggests a lesser density-dependent impact on their distribution, at least within the depth range of ≤730 m. A study investigating density-dependent effects across multiple taxa in areas including Northwest Atlantic also found the area occupied by scorpaeniform fishes to be less positively correlated to abundance than was seen in gadiform and pleuronectiform fishes (Thorson et al., 2016).

Apart from juvenile redfish and small Acadian redfish, most members of Guild III display no significant negative abundance correlations to cod. This could be primarily explained by the size of the fish, as larger redfish are less susceptible to predation. Redfish larger than 27 cm were not found in the stomach contents of cod, as indicated by a previous study in Flemish Cap (Iglesias et al., 2012). This size effect is evident in the significant negative correlation between large cod and small Acadian redfish, which has a smaller upper limit of 21.5 cm. Small beaked redfish and Golden redfish (S. norvegicus), having larger upper size limits of 27.5 cm, show a negative but non-significant correlation to large cod.

The weak but significant positive correlations found between large cod and large redfish might appear unexpected at first glance, as the higher abundance of cod would be expected result in high consumption of juvenile redfish and subsequently lower recruitment in the later years. However, these positive correlations could be influenced by external variables rather than reflecting their trophic relationship, especially the effect of overexploitation which occurred over the same period (Ávila de Melo et al., 2019). Since the surveys began in 1988, redfish SSB underwent a drastic decline up until the mid-1990s, largely due to a sharp increase in fishing mortality (Ávila de Melo et al., 1998). A significant portion of juveniles was also removed as bycatch by shrimp fisheries between 1993 and 1995, contributing to the low recruitment of large redfish in subsequent years (Pérez-Rodríguez et al., 2017). The SSB then rebounded swiftly due to strong recruitment from the early to late-2000s and remained at a high level while fishing mortality was suppressed (Gonçalves et al., 2023). The coincidental collapse of cod and redfish due to the overexploitation, followed by a similar period of recovery, could overshadow the negative impact of cod predation on young redfish. It should be noted that the lag used in the analysis was only one year for assessing the immediate direct impact of cod on other species. The results may not demonstrate the longer-term effect that could be associated to cod such as reduced recruitment that may take more than one year to realise. A further investigation with different lags, depending on life history, to discern cod effects on other species through recruitment should be conducted.

Although unfavourable abiotic conditions were suspected to be the main cause of impaired recruitment, the shift in assemblage and trophic structure after the collapse could also have contribute significantly to the recovery delay (Borovkov et al., 2005; Kuparinen et al., 2014; González-Troncosoet al., 2022). Cod abundance fluctuations may not have a one-way effect on other species, as recent studies showed that food availability can also impact cod recovery (Sguotti et al., 2018; Ruiz-Díaz et al., 2022). The decline of the spawning stock due to the lower food availability could lead to recruitment reduction and delayed recovery (Pérez-Rodríguez et al., 2010). High fishing mortality on their prey could lead to a decrease in cod biomass, making it practically impossible to have high yields and maintain healthy populations of cod, redfish and shrimp concurrently (Pérez-Rodríguez et al., 2022). The replacement of cod by other predators, such as Greenland halibut and wolffishes, which are strongly negatively correlated to cod, might hinder the remaining cod's ability to compete for food and habitat, impeding their immediate reoccupation of the area after the moratorium. Such phenomena are known as trophic cascades, where the trophic structure experiences a fundamental change and is difficult to restore the original system state due to the removal of key species (Fisher et al., 2015; Ripple et al., 2016). Only after the recovery of shrimp, then redfish, to a high abundance level, did the cod population begin to recover(Pérez-Rodríguez et al., 2017).

Management implications

The cod fisheries management in NAFO Division 3M predominantly relies on catch quotas derived from single-species stock assessment (NAFO, 2023a). Since the fishery reopened in 2010, the cod population has been closely monitored to ensure sustainably low fishing mortality (González-Troncoso et al., 2022). While precautionary approaches are considered in management, the reliance on single-species assessments may overlook the ecosystem's complexity, as illustrated by several studies in the area (Pérez-Rodríguez et al., 2011, 2017; Pérez-Rodríguez and Saborido-Rey, 2012; Nogueira et al., 2018). Tests on multiple harvest control rules have emphasized the necessity for target fishing mortalities for cod, redfish and shrimp to be jointly designed (Pérez-Rodríguez et al., 2022). Although the spatial modelling results of this study at this stage do not offer a mechanism to replace the current stock assessment, they do expand the evidence about ecological interactions of cod to include a broader range of species, beyond those commercially targeted. The pressure exerted by cod, whether through competition or direct predation, can significantly impact non-commercially targeted species, including thorny skate, roughhead grenadier, marlin-spike grenadier, longfin hake, wolffishes, Arctic eelpout and presumably some other less abundant species that are not included in this study due to the lack of samples. Although maintaining cod SSB at a safe level should remain a key management priority, the unavoidable impact of cod on other species highlights the necessity of ensuring low fishing pressure on redfish and shrimp to maintain their ecological status (Gonçalves et al., 2023). The ongoing moratorium on shrimp is likely to be crucial for maintaining the structure of the ecosystem - and future reopening should carefully consider the impacts on shrimp biomass as well as their predators (Casas Sánchez, 2023).

Variations in the spatial distribution and abundance of species can offer valuable insights for spatial management(Stelzenmüller et al.<, 2013). For example, the cod population typically resides at the top of the mount, and only spills over into deeper waters when density is high. Protecting the shallowest part could help ensure cod survival, particularly small individuals, and allow catches to only remove the excess fish that spill over from the protected area (Di Lorenzo et al., 2016). The NAFO Scientific Council reports until 2020 showed that the long line fishery consistently operates mainly around the top of Flemish Cap while bottom trawl fisheries were already operating less intensely in the same area, especially those not directly targeting cod (NAFO, 2015, 2020). Therefore, limiting fisheries access to the shallow area may benefit the cod population. Investigation should be carried out into the impacts of closure on the fisheries.

Recent management measures, aiming specifically to protect cod, were implemented in 2021 (NAFO, 2021). First, the closure of the cod fishery during the first quarter of the year, which aligns with the cod spawning period in the area. Second, the employment of sorting grids or grates with a maximum bar spacing of 55 mm for trawls directly targeting cods to reduce small cod bycatch. Since the new measures were put in place, the proportion of the cod catch that is smaller than 45 cm has decreased. However, longlines fishing patterns shifted to concentrate even more in the central part, although the spatial change might be attributed to the lower abundance and TAC rather than solely in response to the measures(Garrido et al., 2023c). This observation suggests that non-spatial measures could reduce juvenile catches but when the stock is in a poor state and concentrated at the top of the cap, the fishing activity would still put more pressure on their refuge area. An explicit special management measure aiming to protect cod, especially when vulnerable, should be considered to ensure their survival and promote recovery.

Fisheries targeting other commercial species, such as the redfish trawl fishery, which catches the highest percentage of cod bycatch, should avoid operating over the shallow areas of less than 270 m depth to reduce small cod bycatch (Garrido et al., 2023c). The shallow area is already not typically densely populated by all redfish species as modelled in this study and observed in other investigations (Nogueira et al., 2017). Therefore, protecting the shallow area from redfish trawls, even partially or temporally, should be considered to help protecting cod. It could follow an example of temporal area closure for shrimp fishery within the 200 m depth contour from 1^st June to 31^st December, which has already been established independently from the overall shrimp moratorium (NAFO, 2023a). However, the socio-economic impact of closure, even temporary, compared to the ecological benefits should be studied for each fishery before decision making.

Protecting marine biodiversity is an important management objective to maintain ecosystem functioning and services, and spatial management is a tool that can help achieve this (Mangano et al., 2015). Some trophic species with inherently lower numbers, including wolffishes, thorny skate, longfin hake, Arctic eelpout and small Greenland halibut, were typically found at the bottom edge of their distribution range when their abundance drastically declined, forming a ring around 400 to 600 m. Accounting for these spatial patterns can be valuable in spatial management efforts to reduce bycatch risk and preserve biodiversity (Komoroske and Lewison, 2015; Hazen et al., 2018). In efforts to conserve biodiversity, NAFO has already instituted numerous closure areas to protect vulnerable marine ecosystems around Flemish Cap, which have been recently reassessed and extended (NAFO, 2022). Most of these closures are situated at depths exceeding 500 m and primarily focus on benthic invertebrate diversity. Integrating the insights from this study could add another layer to future closure assessments, specifically emphasising rare fish species. This inclusion, beyond protecting biodiversity, could also aid vessels in avoiding bycatch that exceeds quotas, potentially reducing the need for time-consuming and costly measures such as exiting the division or conducting tow trials (NAFO, 2023a).

Conclusion

This study utilized a statistical modelling approach to reveal spatiotemporal changes in cod and other demersal species, indicating that fishing and environmental variations affect not only cod abundance but also their distribution range. Cod distribution range greatly contracted when the abundance was low. The results also highlight the secondary effects of cod stock fluctuation on trophic structure and show potential ecosystem-wide impacts from exploiting key species. Nevertheless, not all trophic groups responded uniformly to changes in cod population. The distribution changes for juvenile redfish and northern shrimp ≤730 m depth were small, yet their abundance significantly benefitted from the reduction in cod populations. Cod competitors, including Greenland halibut, thorny skate, roughhead grenadier, marlin-spike grenadier, longfin hake, Atlantic wolffish (A. lupus), spotted wolffish and Arctic eelpout also benefitted from cod depletion. Large redfish do not show a sign of benefiting from cod depletion despite their young being a main cod prey. This could be due to external variables including fishing and environmental pressure having stronger effects. Other species with minimal trophic relationships with cod but residing at a similar depth, such as American plaice and witch flounder, could benefit from cod recovery because their predators could be supressed by cod. These ecological ties should be considered for future wider multispecies interaction modelling and management decision-making.

Acknowledgements

Our sincere thanks to all personnel across all institutes who collected the data in every survey. We thank NAFO for providing access to cod abundance estimates by ages. We also thank CMEMS for making public their data on bottom temperature and NOAA for providing high-resolution bathymetry data. KS thanks the Thai government for funding and the IIM CSIC for hosting his stay.

References

Akaike, H. 2011. Akaike’s Information Criterion. In M. Lovric (Ed.) International Encyclopedia of Statistical Science, 1^st edition, pp. 25–25. Ed. by M. Lovric. https://doi.org/10.1007/978-3-642-04898-2.

Ávila de Melo, A., Saborido-Rey, F., and Alpoim, R. 1998. An Assessment of Redfish in NAFO Div. 3M Including an Approach to Precautionary Management Based on Spawning Biomass and Growth. NAFO SCR Doc., No. 98, Serial No. N3044. 51 p.https://www.nafo.int/Portals/0/PDFs/sc/1998/scr-98-053.pdf.

Ávila de Melo, A., Saborido-Rey, F., Fabeiro, M., Rábade, S., González Troncoso, D., González-Costas, F., Pochtar, M., et al. 2019. An assessment of beaked redfish (S. mentella and S. fasciatus) in NAFO Division 3M (including an update for the most recent level of natural mortality). NAFO SCR Doc., No. 19, Serial No. N6932. 81 p. https://www.nafo.int/Portals/0/PDFs/sc/2019/scr19-016.pdf.

Borovkov, V. A., Vaskov, A. A., and Karsakov, A. L. 2005. The impact of water circulation on the year-class abundance dynamics of redfish and cod on the Flemish Cap. Journal of Northwest Atlantic Fishery Science, 37: 127–134. https://doi.org/10.2960/J.v37.m553.

Casas Sánchez, J. 2023. Division 3M Northern shrimp (Pandalus borealis) – Interim Monitoring Update. NAFO SCR Doc., No. 23, Serial No. N7453. 12 p. https://www.nafo.int/Portals/0/PDFs/sc/2023/scr23-053.pdf.

Cerviño, S., Gil, J., and Sánchez, R. 2005. Changes in Flemish Cap cod distribution and its relationship with environmental changes. NAFO SCR Doc., No. 16, Serial No. N5097. 12 p. https://www.nafo.int/Portals/0/PDFs/sc/2005/scr05-016.pdf

Dawe, E. G., Koen-Alonso, M., Chabot, D., Stansbury, D., and Mullowney, D. 2012. Trophic interactions between key predatory fishes and crustaceans: Comparison of two Northwest Atlantic systems during a period of ecosystem change. Marine Ecology Progress Series, 469: 233–248. https://doi.org/10.3354/meps10136.

Di Lorenzo, M., Claudet, J., and Guidetti, P. 2016. Spillover from marine protected areas to adjacent fisheries has an ecological and a fishery component. Journal for Nature Conservation, 32: 62–66. https://doi.org/10.1016/j.jnc.2016.04.004.

Fisher, J. A., Casini, M., Frank, K. T., Möllmann, C., Leggett, W. C., and Daskalov, G. 2015. The importance of within-system spatial variation in drivers of Marine Ecosystem regime shifts. Philosophical Transactions of the Royal Society B: Biological Sciences, 370: 1–8. https://doi.org/10.1098/rstb.2013.0271.

Garrido, I., González-Troncoso, D., and González-Costas, F. 2023a. Assessment of the cod stock in NAFO Division 3M. NAFO SCR Doc., No. 23, Serial No. N7396. 51 p. https://www.nafo.int/Portals/0/PDFs/sc/2023/scr23-009.pdf.

Garrido, I., González-Troncoso, D., and González-Costas, F. 2023b. Analysis of the Flemish Cap cod fishery: monitoring of the consequences of the management decisions. NAFO SCR Doc., No. 23, Serial No. N7398. 36 p. https://www.nafo.int/Portals/0/PDFs/sc/2023/scr23-011.pdf.

Gonçalves, P., Alpoim, R. and Ávila de Melo, A. 2023. Redfish Div. 3M Biological , Points and advice under the PA alternative framework. NAFO SCR Doc., No. 25, Serial No. N7413. 29 p. https://www.nafo.int/Portals/0/PDFs/sc/2023/scr23-025.pdf.

González, C., Román, E., and Paz, X. 2005. Condition and feeding of American plaice (Hippoglossoides platessoides) in the North Atlantic with emphasis on the Flemish Cap. Journal of Northwest Atlantic Fishery Science, 37: 87–102. https://doi.org/10.2960/J.v37.m556.

González-Troncoso, D., González-Costas, F., and Garrido, I. 2022. Assessment of the Cod Stock in NAFO Division 3M. NAFO SCR Doc., No. 25, Serial No. N7298. 58 p. https://www.nafo.int/Portals/0/PDFs/sc/2022/scr22-025.pdf.

Grüss, A., Drexler, M., and Ainsworth, C. H. 2014. Using delta generalized additive models to produce distribution maps for spatially explicit ecosystem models. Fisheries Research, 159: 11–24. https://doi.org/10.1016/j.fishres.2014.05.005.

Hazen, E. L., Scales, K. L., Maxwell, S. M., Briscoe, D. K., Welch, H., Bograd, S. J., Bailey, H., et al. 2018. A dynamic ocean management tool to reduce bycatch and support sustainable fisheries. Science Advances, 4: 1–7. https://doi.org/10.1126/sciadv.aar3001.

Hendrickson, L., and Vázquez, A. 2005. Density-dependent changes in the spatial distributions of Atlantic cod (Gadus morhua), American plaice (Hippoglossoides platessoides), and Greenland halibut (Reinhardtius hippoglossoides) on the Flemish Cap during 1988–2002. Journal of Northwest Atlantic Fishery Science, 37: 53–72. https://doi.org/10.2960/J.v37.m566.

Hilborn, R. 2011. Future directions in ecosystem based fisheries management: A personal perspective. Fisheries Research, 108: 235–239. https://doi.org/10.1016/j.fishres.2010.12.030.

Hovde, S. C., Albert, O. T., and Nilssen, E. M. 2002. Spatial, seasonal and ontogenetic variation in diet of Northeast Arctic Greenland halibut (Reinhardtius hippoglossoides). ICES Journal of Marine Science, 59: 421–437. https://doi.org/10.1006/jmsc.2002.1171.

Iglesias, G., González-Costas, F., and González-Troncoso, D. 2012. Atlantic cod predation on redfish in Flemish Cap. NAFO SCR Doc., No. 27, Serial No. N6053. 10 p. https://www.nafo.int/Portals/0/PDFs/sc/2012/scr12-027.pdf.

Kindsvater, H. K., and Palkovacs, E. P. 2017. Predicting Eco-evolutionary Impacts of Fishing on Body Size and Trophic Role of Atlantic Cod. Copeia, 105: 475–482. https://doi.org/10.1643/OT-16-533.

Komoroske, L. M., and Lewison, R. L. 2015. Addressing fisheries bycatch in a changing world. Frontiers in Marine Science, 2: 1–11. https://doi.org/10.3389/fmars.2015.00083.

Kuparinen, A., Keith, D. M., and Hutchings, J. A. 2014. Increased environmentally driven recruitment variability decreases resilience to fishing and increases uncertainty of recovery. ICES Journal of Marine Science, 71: 1507–1514. https://doi.org/10.1093/icesjms/fsu021.

Lin, X., and Zhang, D. 1999. Inference in generalized additive mixed models by using smoothing splines. Journal of the Royal Statistical Society. Series B: Statistical Methodology, 61: 381–400. https://doi.org/10.1111/1467-9868.00183.

Link, J. S., Bolles, K., and Milliken, C. G. 2002. The feeding ecology of flatfish in the Northwest Atlantic. Journal of Northwest Atlantic Fishery Science, 22: 1–17. https://doi.org/10.2960/J.v30.a1.

Mangano, M. C., O’Leary, B. C., Mirto, S., Mazzola, A., and Sarà, G. 2015. The comparative biological effects of spatial management measures in protecting marine biodiversity: A systematic review protocol. Environmental Evidence, 4: 1–8. https://environmentalevidencejournal.biomedcentral.com/articles/10.1186/s13750-015-0047-2

Murillo, F. J., Durán Muñoz, P., Altuna, A., and Serrano, A. 2011. Distribution of deep-water corals of the Flemish Cap, Flemish Pass, and the Grand Banks of Newfoundland (Northwest Atlantic Ocean): Interaction with fishing activities. ICES Journal of Marine Science, 68: 319–332. https://doi.org/10.1093/icesjms/fsq071.

Murillo, F. J., Serrano, A., Kenchington, E., and Mora, J. 2016. Epibenthic assemblages of the Tail of the Grand Bank and Flemish Cap (northwest Atlantic) in relation to environmental parameters and trawling intensity. Deep-Sea Research Part I: Oceanographic Research Papers, 109: 99–122. https://doi.org/10.1016/j.dsr.2015.08.006.

Murillo, F. J., Weigel, B., Bouchard Marmen, M., and Kenchington, E. 2020a. Marine epibenthic functional diversity on Flemish Cap (north-west Atlantic)—Identifying trait responses to the environment and mapping ecosystem functions. Diversity and Distributions, 26: 460–478. https://doi.org/10.1111/ddi.13026.

Murillo, F. J., Kenchington, E., Koen-Alonso, M., Guijarro, J., Kenchington, T. J., Sacau, M., Beazley, L., et al. 2020b. Mapping benthic ecological diversity and interactions with bottom-contact fishing on the Flemish Cap (northwest Atlantic). Ecological Indicators, 112: 1–13. https://doi.org/10.1016/j.ecolind.2020.106135.

NAFO. 2015. Report of the 8th Meeting of the NAFO Scientific Council (SC) Working Group on Ecosystem Science and Assessment (WG-ESA). NAFO SCS Doc., No. 19, Serial No. N6549. 176 p. https://www.nafo.int/Portals/0/PDFs/sc/2015/scs15-19.pdf.

NAFO. 2019. Report of the NAFO Scientific Council Flemish Cap (NAFO Div. 3M) Cod Stock Management Strategy Evaluation (MSE). NAFO SCS Doc. No. 4 Serial No. N6911 34 p. https://www.nafo.int/Portals/0/PDFs/sc/2019/scs19-04.pdf.

NAFO. 2020. Report of the 13th Meeting of the NAFO Scientific Council Working Group on Ecosystem Science and Assessment (WG-ESA). NAFO SCS Doc., No. 23, Serial No. 7148. 270 p. https://www.nafo.int/Portals/0/PDFs/sc/2020/scs20-23.pdf.

NAFO. 2021. Northwest Atlantic Fisheries Organization Conservation and Enforcement Measure 2021. NAFO/COM Doc., No. 1, Serial No. N7153. 194 p. https://www.nafo.int/Portals/0/PDFs/COM/2021/comdoc21-01.pdf.

NAFO. 2022. Report of the 15th Meeting of the NAFO Scientific Council Working Group on Ecosystem Science and Assessment (WG-ESA). NAFO SCS Doc., No. 25, Serial No. N7367. 90 p. https://www.nafo.int/Portals/0/PDFs/sc/2022/scs22-25.pdf.

NAFO. 2023a. Northwest Atlantic Fisheries Organization Conservation and Enforcement Measures 2023. NAFO/COM Doc., No. 1, Serial No. N7368. 204 p. https://www.nafo.int/Portals/0/PDFs/COM/2023/comdoc23-01.pdf.

NAFO. 2023b. NAFO Precautionary Approach Working Group (PA-WG) Meeting Report. NAFO SCS Doc., No. 7, Serial No. N7382. 14 p. https://www.nafo.int/Portals/0/PDFs/sc/2023/scs23-07.pdf.

NOAA. 2023. ETOPO 2022 15 Arc-Second Global Relief Model. NOAA National Centers for Environmental Information. https://doi.org/10.25921/fd45-gt74.

Nogueira, A., González-Troncoso, D., and Tolimieri, N. 2016. Changes and trends in the overexploited fish assemblages of two fishing grounds of the Northwest Atlantic. ICES Journal of Marine Science, 73: 345–358. https://doi.org/10.1093/icesjms/fsv172.

Nogueira, A., Paz, X., and González-Troncoso, D. 2017. Demersal groundfish assemblages and depth-related trends on Flemish Cap (NAFO division 3M): 2004–2013. Fisheries Research, 186: 192–204. https://doi.org/10.1016/j.fishres.2016.08.016.

Nogueira, A., Pérez-Rodríguez, A., González-Troncoso, D., and Saborido-Rey, F. 2018. Could population and community indicators contribute to identify the driver factors and describe the dynamic in the Flemish Cap demersal assemblages? Fisheries Research, 204: 33–40. https://doi.org/10.1016/j.fishres.2018.01.019.

Orio, A., Bergström, U., Florin, A. B., Lehmann, A., Šics, I., and Casini, M. 2019. Spatial contraction of demersal fish populations in a large marine ecosystem. Journal of Biogeography, 46: 633–645. https://doi.org/10.1111/jbi.13510.

Pérez-Rodríguez, A., Morgan, M., Rideout, R., Domínguez-Petit, R., and Saborido-Rey, F. 2010. Study of the relationship between total egg production, female spawning stock biomass, and recruitment of Flemish Cap cod (Gadus morhua). Ciencias Marinas, 37: 675–687. https://doi.org/10.7773/cm.v37i4b.1785.

Pérez-Rodríguez, A., Koen-Alonso, M., González-Iglesias, C., and Saborido-Rey, F. 2011. Analysis of common trends in the feeding habits of main demersal fish species on the Flemish Cap. NAFO SCR Doc., No. 77, Serial No. N6009. 30 p. https://www.nafo.int/Portals/0/PDFs/sc/2011/scr11-077.pdf.

Pérez-Rodríguez, A., and Saborido-Rey, F. 2012. Food consumption of Flemish Cap cod (Gadus morhua) and redfish (Sebastes spp.) using generic bioenergetic models. NAFO SCS Doc., No. 68, Serial No. N6136. 15 p. https://www.nafo.int/Portals/0/PDFs/sc/2012/scr12-068.pdf.

Pérez-Rodríguez, A., Koen-Alonso, M., and Saborido-Rey, F. 2012. Changes and trends in the demersal fish community of the Flemish Cap, Northwest Atlantic, in the period 1988–2008. ICES Journal of Marine Science, 69: 902–912. https://doi.org/10.1093/icesjms/fss019.

Pérez-Rodríguez, A., Howell, D., Casas, M., Saborido-Rey, F., and Ávila-De Melo, A. 2017. Dynamic of the Flemish cap commercial stocks: Use of a gadget multispecies model to determine the relevance and synergies among predation, recruitment, and fishing. Canadian Journal of Fisheries and Aquatic Sciences, 74: 582–597. https://doi.org/10.1139/cjfas-2016-0111.

Pérez-Rodríguez, A., Umar, I., Goto, D., Howell, D., Mosqueira, I., and González-Troncoso, D. 2022. Evaluation of harvest control rules for a group of interacting commercial stocks using a multispecies MSE framework. Canadian Journal of Fisheries and Aquatic Sciences, 79: 1302–1320. https://doi.org/10.1139/cjfas-2021-0069.

Pham, C. K., Murillo, F. J., Lirette, C., Maldonado, M., Colaço, A., Ottaviani, D., and Kenchington, E. 2019. Removal of deep-sea sponges by bottom trawling in the Flemish Cap area: conservation, ecology and economic assessment. Scientific Reports, 9: 1–13. https://doi.org/10.1038/s41598-019-52250-1.

Planque, B., Loots, C., Petitgas, P., Lindstrøm, U., and Vaz, S. 2011. Understanding what controls the spatial distribution of fish populations using a multi-model approach. Fisheries Oceanography, 20: 1–17. https://doi.org/10.1111/j.1365-2419.2010.00546.x.

Ripple, W. J., Estes, J. A., Schmitz, O. J., Constant, V., Kaylor, M. J., Lenz, A., Motley, J. L., et al. 2016. What is a trophic cascade? Trends in Ecology and Evolution, 31: 842–849. https://doi.org/10.1016/j.tree.2016.08.010.

Rose, G. A., DeYoung, B., Kulka, D. W., Goddard, S. V., and Fletcher, G. L. 2000. Distribution shifts and overfishing the northern cod (Gadus morhua): A view from the ocean. Canadian Journal of Fisheries and Aquatic Sciences, 57: 644–663. https://doi.org/10.1139/f00-004.

Ruiz-Díaz, R., Dominguez-Petit, R., and Saborido-Rey, F. 2022. Atlantic Cod Growth History in Flemish Cap Between 1981 and 2016: The Impact of Fishing and Climate on Growth Performance. Frontiers in Marine Science, 9: 1–15. https://doi.org/10.3389/fmars.2022.876488.

Sguotti, C., Otto, S. A., Frelat, R., Langbehn, T. J., Ryberg, M. P., Lindegren, M., Durant, J. M., et al. 2018. Catastrophic dynamics limit Atlantic cod recovery. Proceedings of the Royal Society B: Biological Sciences, 286: 1–10. https://doi.org/10.1098/rspb.2018.2877.

Stelzenmüller, V., Lee, J., South, A., Foden, J., and Rogers, S. I. 2013. Practical tools to support marine spatial planning: A review and some prototype tools. Marine Policy, 38: 214–227. https://doi.org/10.1016/j.marpol.2012.05.038.

Thorson, J. T., Rindorf, A., Gao, J., Hanselman, D. H., and Winker, H. 2016. Density-dependent changes in effective area occupied for sea-bottom-associated marine fishes. Proceedings of the Royal Society B: Biological Sciences, 283: 1–10. https://doi.org/10.1098/rspb.2016.1853.

Vázquez, A., Casas, J. M., Brodie, W. B., Murillo, F. J., Mandado, M., Gago, A., Alpoim, R., et al. 2013. List of species as recorded by Canadian and EU Bottom Trawl Surveys in Flemish Cap. No. 5., Serial No. N6154. 13 p. https://www.nafo.int/Portals/0/PDFs/sc/2013/scr13-005.pdf.

Vázquez, A., Miguel Casas, J., and Alpoim, R. 2014. Protocols of the EU bottom trawl survey of Flemish Cap. NAFO Scientific Council Studies, 46: 1–51. https://doi.org/10.2960/S.v46.m1

Walters, C., and Maguire, J. J. 1996. Lessons for stock assessment from the northern cod collapse. Reviews in Fish Biology and Fisheries, 6: 125–137. https://doi.org/10.1007/BF00182340.

Wood, S. N. 2023, December 21. mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation. https://doi.org/10.32614/CRAN.package.mgcv.

Yasuhara, M. and Danovaro, R. 2016. Temperature impacts on deep-sea biodiversity. Biological Reviews, 91: 275–287. https://doi.org/10.1111/brv.12169.

Citation: Songin, K., Pierce, G., and Saborido-Rey, F. 2024. Spatiotemporal changes in the Atlantic cod (Gadus morhua) stock at Flemish Cap (1993–2019) and their relationships with demersal communities.
J. Northw. Atl. Fish. Sci., 55: 59–78. https://doi.org/10.2960/J.v55.m748