a, Bathymetry with an overlay showing locations of the 84 1° × 1° latitude–longitude cells in which fish abundance, distribution and sea temperature were reported and predicted. b, Mean sea surface temperature (SST, red) and near-bottom temperature (NBT, black) in the study cells from 1980–2060 in summer (July–September, solid line) and winter (January–March, dashed line) from the QUMP_ens_00 northwest European shelf seas climate model. Mean decadal values (as used in the model) are overlaid in the corresponding colours for SST and NBT for each season.
Fish surveys.
We used two long-term monitoring surveys that give detailed descriptions of the distribution and abundance of demersal (bottom-dwelling) fish in the North Sea. The Centre for Environment, Fisheries and Aquaculture Science UK (Cefas) time series is a summer survey (August–September) conducted since 1980. The survey encompasses 69 1° × 1° latitude–longitude cells with at least three hauls conducted in each decade. The International Council for the Exploration of the Sea (ICES) International Bottom Trawl Survey (IBTS) time series is a winter survey (January–March) conducted since 1980. The survey encompasses 84 1° × 1° cells with at least three hauls conducted in each decade. Both surveys are conducted using otter trawling gear (Granton trawl for pre-1992 Cefas surveys, otherwise Grande Ouverture Verticale (GOV) trawls). Raw catch data were fourth-root transformed to reduce skewness that is inherent in ecological abundance data.
Our study focused on the ten most abundant demersal species targeted by commercial fisheries or taken as bycatch (Fig. 2c), which together accounted for 68% of commercial landings (by weight) in the North Sea fishery from 1980–2010 (www.ices.dk/marine-data/dataset-collections/Pages/Fish-catch-and-stock-assessment.aspx). For both surveys, we grouped data into three ten-year time slices and one three-year time slice for the analyses: 1980–1989, 1990–1999, 2000–2009 and 2010–2012. The limited 2010–2012 time slice was used only for testing predictions from the GAMs. To ensure a balanced design, mean values for each decadal time period were used. This method controls for the variable numbers of survey hauls taken in each cell and ensures that longer-term responses to climate change are identified rather than year-on-year variability. All data were fourth root transformed before being subject to GAM modelling, and individual cell predictions were back-transformed before calculation of correlation coefficients.
Depth.
We used mean 1° × 1° cell in situ measures of depth taken during the hauls for each survey (Supplementary Fig. 3), which closely matched data from the 1° × 1° resolution GEBCO Digital Atlas (summer survey, r = 0.91; winter survey, r = 0.90; www.gebco.net/data_and_products/gebco_digital_atlas)3.
Temperature and salinity.
We calculated sea surface temperature (SST), near-bottom temperature (NBT) and salinity (Supplementary Fig. 3) for the period 1980–2012 using the UK Meteorological Office Hadley Centre QUMP_ens_00 standard model for the northwest European shelf seas. Modelled temperatures closely matched data from the Hadley Centre global ocean surface temperature database (HadISST1.1; 92 cells, Pearson’s r = 0.84; www.metoffice.gov.uk/hadobs/hadisst). Data from the QUMP_ens_00 model were provided as monthly means for 1° × 1° cells, enabling mean winter (January–March), summer (July–September) and mean annual values to be calculated (Fig. 1).
Fishing pressure.
We calculated a spatially explicit metric of fishing pressure for each ten-year time slice by combining annual multispecies fishing mortality (F) estimates for North Sea demersal species (mean estimates of regional F for cod, dab, haddock, hake, lemon sole, ling, long rough dab, plaice, saithe and whiting, weighted by spawning-stock biomass, from ICES stock assessments; www.ices.dk/datacentre/StdGraphDB.asp)3 with mean otter and beam trawling effort for each 1° × 1° cell based on hours of fishing27 (Supplementary Fig. 3). This integrated metric, combining temporal trends in fishing mortality and spatial distribution of fishing effort, enabled us to test the importance of fishing pressure as a predictor of abundance.
Identifying key predictors.
We used GAM models, coded using the mgcv package in R (www.r-project.org), to test the performance of GAMs for predicting changes in fish species’ distribution and identify the importance of different variables to these predictions. The s smooth was used with k = 7 for all variables to limit the degrees of freedom in line with the number of data points. The Gaussian model was used. Assessment of the plots for each variable using the gam.plot function showed that increasing the k value did not improve the model fit to each variable. The gam.check function was used to check the k index was above or close to 1, with non-significant p values. Analysis of the residuals showed no obvious deviations from normal distributions, and the response to fitted values relationship was close to linear.
Data from 2000–2009 were used to test sets of variables, as this period had the greatest survey intensity. To identify variables that most strongly influenced prediction we first developed a model with all variables (annual temperatures, seasonal temperatures, depth, salinity and fishing), and a subsequent five models each excluding one set of variables (Supplementary Table 1). Sea surface and near-bottom temperatures from both the summer and winter were grouped together to characterize seasonal fluctuations. This suite of potentially correlated variables captured the extremes of temperatures that all species may experience at different life stages, and ensured that thermal conditions with and without the seasonal thermocline, annually varying ocean currents and land mass effects are all included. We compared the performance of models based on: the strength of correlation r between observed and predicted data, weighted AIC (ref. 28) using data from the AIC function in R, and using generalized cross validation (GCV, through summary.gam in R). Inclusion of interaction terms between depth and seasonal temperature extremes either reduced or had little influence on model performance (Supplementary Table 2 and summaries based on Akaike weights in Supplementary Fig. 4).
Model development.
We developed predictive GAMs with a set of variables that were effective across all species. The correlation coefficient r, AIC values and GCV values of modelled and observed data were compared. Across-species inclusion of depth, seasonal temperature, annual temperature, salinity and fishing effort all improved the predictions (Fig. 2a). A key finding from this model development stage is that variables that are readily measured and projected in climate models effectively predict species distributions. On average, models that excluded fishing effort were most similar to the all-variable models (Supplementary Table 1 and Fig. 2a). As this metric had little predictive value, and we have no robust models of future fishing effort, we excluded it when making future predictions.
Training period and predictive performance.
To assess the influence of the duration of training data on predictive power, GAMs trained on sets of one, two and three decades of data for each species were used to predict ten years into the future (Supplementary Fig. 1), and the associations between predicted and known data compared. We also assessed the performance of the model to predict further into the future within the historic records available (Supplementary Fig. 1). We compared predictions with known abundance data for each species for each forecasting period (0 to 30 years).
Forecasting future distributions.
We used surface and near-bottom annual and seasonal temperature projections from the QUMP_ens_00 model, surface and near-bottom salinity, and average depths from surveys between 1980–2012 as the environmental variables for our predictions. We predicted fish abundances for sequential decades from 2000–2009 to 2050–2059 (Supplementary Figs 5 and 6) using environmental variables (Supplementary Figs 3 and 7), and observed fish abundances from 2000–2009. Throughout the projection period many cells do not experience temperatures outside of the range used to train the model (Supplementary Fig. 2). For the widespread species in this study it is therefore likely that at least parts of the population have experienced future conditions. However, we recognize that in future projected conditions the climate in some areas of the North Sea will depart from existing variability in the model training period. As it is not possible to test the model beyond present thermal conditions using known data, some caution should be taken in interpreting projections for cells as they begin to experience temperatures beyond those at present in the region (Supplementary Fig. 2).