Sea temperature is the primary driver of recent and predicted fish community structure across Northeast Atlantic shelf seas

Climate change has strongly influenced the distribution and abundance of marine fish species, leading to concern about effects of future climate on commercially harvested stocks. Understanding the key drivers of large‐scale spatial variation across present‐day marine assemblages enables predictions of future change. Here we present a unique analysis of standardised abundance data for 198 marine fish species from across the Northeast Atlantic collected by 23 surveys and 31,502 sampling events between 2005 and 2018. Our analyses of the spatially comprehensive standardised data identified temperature as the key driver of fish community structure across the region, followed by salinity and depth. We employed these key environmental variables to model how climate change will affect both the distributions of individual species and local community structure for the years 2050 and 2100 under multiple emissions scenarios. Our results consistently indicate that projected climate change will lead to shifts in species communities across the entire region. Overall, the greatest community‐level changes are predicted at locations with greater warming, with the most pronounced effects at higher latitudes. Based on these results, we suggest that future climate‐driven warming will lead to widespread changes in opportunities for commercial fisheries across the region.

Sectors of the Northeast Atlantic have warmed by over 1°C over the last 50 years, and further substantive warming is projected by 2069-2098 (Tinker et al., 2016). Studies have already shown changes to species distributions and communities in these waters (Beaugrand et al., 2009;Bedford et al., 2020;Chivers et al., 2017;Edwards et al., 2021;Simpson et al., 2011). Recent work indicates that fish species in the region will undergo poleward movements to track preferred thermal habitats under further warming (Cheung et al., 2010(Cheung et al., , 2011Fernandes et al., 2020;Jones et al., 2013), consistent with predictions across the globe (Cheung et al., 2009;Okunishi et al., 2012). However, due to highly specific depth and habitat requirements, the future distributions and abundance of species may be prevented from tracking their optimal thermal niche, including in the North Sea (Rutterford et al., 2015).
Thus, many species may need to acclimate or adapt to changing thermal environments to avoid suffering population abundance declines (Rutterford et al., 2015). Such climate-driven changes to the abundance and distributions of commercial importance have the potential to drive socio-economic impacts (Montero-Serra et al., 2015;Payne et al., 2021;Simpson et al., 2011;Teixeira et al., 2014;Wabnitz et al., 2018) and fisheries may need to undertake adaptation measures to minimise detrimental impacts (Garza-Gil et al., 2011;Jones et al., 2014).
Driven by the socio-economic importance of Northeast Atlantic shelf-sea fisheries, spatially explicit trawl surveys of regional marine fish communities have been undertaken since the 1960s. However, each trawl survey only covers a discrete section of the Northeast Atlantic continental shelf, and surveys also vary in sampling methods and seasonal coverage (Cadigan & Dowden, 2010;Harley & Myers, 2001;ICES, 2020). Analyses of wider spatial patterns therefore require data standardisation. A straightforward approach has been to use simplified metrics such as presence-absence (Simpson et al., 2011), or frequency-of-occurrence within catches (Montero-Serra et al., 2015). Although these methods can provide valuable insights into species distributions, ideally analyses of abundance data would also account for variation in differences in catch efficiency of each species between gears (Moriarty et al., 2020;Walker et al., 2017). Here we study relative abundance using a method that explicitly accounts for the effects of contrasting survey methods in single species statistical models using data from geographically overlapping surveys (Maltby et al., 2020). We then use these data to determine the key environmental predictors of marine fish community structure across the region, and to predict the effects of projected climate scenarios from 2050 to 2100.

| Fish community survey data
We collated data from 23 fisheries trawl surveys from across the shelf seas of the Northeast Atlantic, from southern Spain to northern Norway. The study area includes the Baltic Sea, but excludes the Barents Sea and Mediterranean Sea. To include surveys in the standardisation analysis, each had to overlap with at least one other survey ( Figure S1 and Table S1). To ensure adequate survey overlap, two surveys with similar gear and sampling seasons were grouped (from Western Scotland and Norway respectively; Table S1). We selected the components of surveys that took place from 2005 to 2018, inclusive, to align the fish abundance data with available environmental data for the study region. Only hauls with a duration greater than 5 min were included. In total, 31,502 survey hauls were included in the analysis, using data recorded at the coordinates of the haul.
Survey hauls were assigned to 1° latitude × 1° longitude grid cells, to align with environmental data that were available at the same spatial scale. Only those grid cells with at least three hauls from each survey undertaken in the area were included ( Figure 1). This enabled us to maintain high spatial coverage, while avoiding cells with the lowest sampling effort (Rutterford et al. 2015). One grid cell (50.5°N, 4.5°W) was treated as two cells in the analysis, as it was fully split by land with different surveys on the northern and southern side.
Average catch per unit effort (CPUE, individuals per hour) was calculated for each species for individual hauls in each grid cell.
These abundance data were then fourth-root transformed to reduce the influence of outliers. Where a grid cell was sampled, but a species not recorded throughout the study period zero abundance was assumed. Only species-level data from individuals in the following classes were included: ray-finned fishes, elasmobranchs, chimaeras, hagfishes and lampreys. Where species names were not consistent, the accepted species names were assigned using the World Register of Marine Species database (WoRMS, 2019). The analysed data included survey data for 487 species for 340 grid cells spanning 36° N to 71° N and 15° W to 31° E ( Figure S2).

| Standardised fish abundance
We standardised fish species abundance data to control for variation between survey gears and methods using a general linear model (CPUE ~grid cell + survey). Using the lsmeans package version 2.30.0 (Lenth, 2016) in R Core Team (2018), we generated least-square mean CPUE estimates per grid cell and survey (following Maltby et al., 2020). We then compared the modelled estimates to reported catches using Spearman's rank correlations for all 487 species reported across the study area. Only species where the average r value between the modelled data and the survey catch-across all surveys-was >0.4, and total catch exceeded 100 individuals, were considered to have a model fit that reliably estimated catch levels. Species where CPUE was only recorded in one survey or less than three grid cells were removed from the analysis at this stage. We then used the model to generate standardised least-square mean species CPUE estimates for each grid cell ( Figure S2), limited to zero. After this process, a total of 198 of a total 487 species were standardised using a general linear model and a post hoc calculation of least-square means. These 198 species comprised over 99% of total fish (individuals) caught and reported in all surveys.

| Environmental data (2005-2017)
World Oceans Atlas environmental data from the National Oceanic and Atmospheric Administration (NOAA) were available for 193 grid cells overlapping with fisheries survey data including: dissolved oxygen (Garcia et al., 2018b), nitrate (μmol/kg), phosphate (μmol/kg) (Garcia et al., 2018a), salinity (PSU)  and temperature (°C) at depth data  ( Figure 1). To reduce the influence of outliers, several variables were log 10 transformed, including depth, phosphate at depth and surface dissolved oxygen data. Mean annual temperature and salinity data were available for the specific time period 2005-2017, other variables were based on mean annual values from available observations. Data from 0 m (surface) and the maximum depth were used. Maximum depth of cells was extracted from Bio-ORACLE (Tyberghein et al., 2012) and married with NOAA records (as above) for that specific depth band (data sources detailed in Table S1).

| Principal coordinate analysis of species assemblage
Principal coordinate analysis (PCoA) based on a Bray-Curtis similarity matrix, was used to identify the primary axes of spatial variation across the fish species assemblage using the R packages vegan version 2.5.7 (Oksanen et al., 2020) and ape version 5.4.1 (Paradis & Schliep, 2019). The eig_perc function (Mikryukov & Mahé, 2021) was used to estimate the proportion of variance in species composition explained by the PCoA. PCoA was used as a repeatable ordination method that was reliably able to capture variation in only the species abundance matrix. Multiple linear regression models were generated to identify potential environmental drivers of variation in community structure for PCoA axes 1, 2 and 3. Variables were selected using stepwise AIC forward regression, and using one step forward in the R package olsrr version 0.5.3 (Hebbali, 2020).
The standardised species CPUE data were split into functional groupings, and again the primary axes of spatial variation were identified using PCoA, with potential environmental drivers of variation PCoA axes 1, 2 and 3 tested using multiple regression linear models. These functional groupings included: pelagic (70 species), demersal (125 species), small size at maturity <30 cm (60 species) and large size at maturity >30 cm (45 species), as extracted from Fishbase (2019) using the R package rfishbase (Boettiger et al., 2012).
To explore variation between exploited and non-exploited species, average catches for the period 2006-2018 from official nominal catch data (ICES, 2006(ICES, -2018 were used to categorise species as non-exploited <100 t (23 species) and commercially exploited (76) landings >10,000 t. Species landed in quantities between these values (>100 and < 10,000 t) were not assessed to ensure that the analysis only captured species that can be clearly defined as commercially exploited or non-exploited (see associated data files for species groupings).

| Projecting future abundances using generalised additive models
Generalised additive models (GAMs) were used to predict species CPUE for the 193 cells of the study area using sea surface temperature (SST), sea surface salinity (SSS) from 2002 to 2009, available at the 1° latitude x 1° longitude grid cell scale, and mean cell depth from Bio-ORACLE (Tyberghein et al., 2012), using the following model structure and R package mgcv version 1.8.31 (Wood, 2011).
For each species the ability of the model to predict the training data was tested, and species were included where r values were >0.3 alongside significant p-values (<.05) for at least one of the predictor variables (Table S2). These models were able to reliably selfpredict abundances for 152 of the species (88%) ( Table S2).
For those 152 species we used the GAM model to predict future CPUE for each species and cell of the study area abundance.
Bio-ORACLE data were the only dataset available for the spatial extent of the study at the time of this analysis. RCPs represent possible future greenhouse gas concentrations and are used to describe different possible climate futures, ranging from a stringent mitigation scenario (RCP2.6), two intermediate scenarios (RCP4.5 and RCP6.0) and one scenario with very high GHG emissions (RCP8.5) (IPCC, 2014). These were used to explore the range of potential conditions, and associated climate impacts, the marine environment may face. The GAM prediction data for all time periods for the 152 species were analysed using PCoA, as described above.

| Community structure (2005-2018)
Our reconstructions of community structure using PCoA of standardised data showed clear gradients along the primary axes of variation (PCoA 1-3; Figure 2). Sea surface temperature was the environmental variable that most closely associated with the first axis of variation (PCoA axis 1), which explained 29% of variation in the data (Table 1). Surface salinity was most closely associated with the second axis of variation (PCoA axis 2), which explained 21% of variation in the data ( Table 1). The steepest salinity gradient across the study region was between the high-salinity Kattegat and the lowsalinity Baltic Sea (Berg et al., 2015). Depth was most closely associated with the third axis of variation (PCoA axis 3), which explained 12% of variation in the data ( Table 1). The studied depth variation was substantial, ranging from 10 m coastal waters to 4900 m off the shelf edge. When data were restricted to species groupings based on their ecological and life-history traits (pelagic, demersal, small size at maturity, large size at maturity, high commercial catch, low commercial catch), again we consistently found strong associations between the primary axis of variation PCoA1 and sea surface temperature (Table S3).

| Community structure (2050 and 2100)
The dominant roles of sea surface temperature, salinity and depth on fish assemblage structure justifies their use as predictors in abundance-based single species generalised additive models (GAMs). We predicted changes in local community composition towards those associated with warmer conditions under all RCP scenarios to 2100, except for one grid cell (55.5 N, −8.5 E) in the scenario with the lowest forcing (Figure 3a, RCP 2.6). Overall, the model predicted greater shifts in community structure, as measured using PCoA axis scores, under scenarios with greater climate forcing. In the region of study, RCP 6.0 caused lower climate change forcing than RCP 2.6 and 4.5 up to 2050, although greater forcing by 2100 ( Figures S3 and Figure S4, Figure 3), which explained why this higher emissions scenario predicted lower impacts on species assemblages than RCP 2.6 and 4.5 prior to 2100. Higher latitude areas with the lowest temperatures were projected to experience the most warming and greatest change in community assemblages (Figure 3).

| DISCUSS ION
Collectively, the use of large-scale fisheries-independent data has enabled us to map the major gradients in the abundance of marine fish communities across the Northeast Atlantic. Within the sampled region, we found temperature to be the primary driver of community structure, consistent with previous studies across subsets of the study area including the North Sea (Dulvy et al., 2008;Rutterford et al., 2015), southwestern UK waters (Genner et al., 2004;Maltby et al., 2020) and across broader regions of the European shelf (Montero-Serra et al., 2015;Simpson et al., 2011). We also found salinity to have an important role in structuring the fish assemblage of the region, largely by limiting the distributions of many fully marine species (Hiddink & Coleby, 2012;Mackenzie et al., 2007). The brackish environment of the Baltic tended to favour low-salinity species such as European smelt (Osmerus eperlanus) (Sendek & Bogdanov, 2019), as well as having relatively large abundance of marine species that have been able to locally adapt to the low-salinity environment, such as Atlantic cod (Gadus morhua) (Berg et al., 2015) F I G U R E 2 Spatial variation in the structure of the Northeast Atlantic fish assemblage (198 species, 193 grid cells) as captured by principal coordinates analysis (PCoA) scores for (a) PCoA axis 1 that primarily associates with sea surface temperature, (b) PCoA axis 2 that primarily associates with sea surface salinity and (c) PCoA axis 3 that primarily associates with depth (Table 1). and European flounder (Platichthys flesus) (Florin & Hoglund, 2008).
Although salinity gradients are most extreme in the Baltic, salinity was also associated, to a lesser extent, with variation across the study area when the Baltic was excluded (Table S4). We also found depth to have an important role in structuring community composition. The key role of depth in limiting the abundance of species in the region is well understood (Campbell et al., 2011;Magnussen, 2002) and linked to physiological adaptations to specific environmental conditions (Brown & Thatje, 2014). For example, deeper water species have visual systems adapted to low light levels (de Busserolles et al., 2020), in addition to physiological adaptations enabling tolerance to high hydrostatic pressure (Somero, 1992) and cooler temperatures (Baalsrud et al., 2017;Brown & Thatje, 2014).
Nitrate levels were also found to have an association with community structure (Table 1). Marine nitrates come from multiple sources and in our study region, high nitrate levels are present in the southern North Sea, which has been linked to terrestrial runoff (Desmit et al., 2020;Menesguen et al., 2018). Additionally, high nitrates are present in deep sea environments off the continental shelf, linked to natural nitrification that takes place below the euphotic zone (Dore & Karl, 1996). Although nitrate levels are unlikely to affect fish species directly it may be that factors covarying with nitrate influence community assemblages. For example, nutrient inputs alter phytoplankton productivity and oxygen levels in the southern North Sea (Capuzzo et al., 2018) and upwelling areas (Alonso-Perez and Castro, 2014), and riverine inputs reduce salinity while increasing nutrient concentrations (Desmit et al., 2015).
The resolution of large-scale assemblage-level patterns in species abundance in our study region have hitherto been limited to planktonic species, sampled using the Continuous Plankton Recorder (Beaugrand et al., 2009;Bedford et al., 2020;Edwards et al., 2021). Those studies have identified clear environmentally driven gradients in the distributions of species, and enabled detailed comparisons of the shifts in assemblage composition in the region in response to climate change (Beaugrand et al., 2009;Bedford et al., 2020). In this study, we were able to identify the major drivers of spatial patterns of variation in the fish assem- comparisons. To ensure adequate replication and overlap of survey data, and to align with climate models for the area, we were limited to one time period of fish survey data (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). As survey data time series extend and become increasingly accessible (Moriarty et al., 2017), and further datasets are made available (Maureaud et al., 2021), there will be enhanced opportunities to explore temporal changes in the distribution of fish species, and to further refine climate change prediction models. are likely to depend on many more variables. Habitat availability (Rutterford et al., 2015), including at different life-history stages (Attrill & Power, 2002;Hermant et al., 2010), as well as dispersal capacity (Sunday et al., 2015), fisheries impacts and ecological interactions (Genner et al., 2004), are likely to combine to determine the capacities of species for distributional shifts. Our data-driven analytical approach assumes that models capture the outcomes of biotic and abiotic interactions, without relying on detailed understanding of the causative processes. Certainly the use of further predictor variable data, potentially at the sub-regional scales, and finer resolutions, are likely to improve the variance explained by models; for example, incorporating substrate characteristics  and past and present fishing activity (Jennings & Blanchard, 2004).
However, incorporating variables such as habitat dependencies is complex as it can vary across key life stages for each species, is not established for many species and is also spatially variable. These issues make habitat difficult to capture on a grid cell scale. Moreover, a previous study across a section of our study area showed that available habitat metrics did not add substantial explanatory power to species distribution models of demersal fish species (Maltby et al., 2020).
There are additional limitations when predicting future distributions of species communities using available climate projections and fisheries survey data, due to data availability. Climate model projections were not reliably available for many of the coastal cells in the analysis, which meant that shallower and often more variable environments were not captured in the analysis. Many fish species depend on coastal areas for at least part of their life history and so improvements to climate models should enhance predictions for these ecologically and commercially important areas. Additionally, the trawl survey sampling frequency has not been consistent across the study area, particularly in the more offshore and northerly sectors of our study area. This likely increases the variance around species abundance estimates in these areas, and more consistent sampling would enable more robust comparisons to be made across the whole sampling area.
Uncertainty in climate and fish distribution data needs to be considered when generating model predictions. Species distribution projections rely on the assumption that climate variables are accurately captured in climate models, and previous projections of fish distributions under climate change have shown considerable variation when using data from different climate ensemble models, particularly as models project further into the future (Maltby et al., 2020). Nevertheless, the use of ensembles of climate models that use different approaches could help to identify locations where results are most consistent and likely reliable. There may also be variation in predicted abundance depending on the species distribution modelling approach. Here we used GAM models, but the use of more than one modelling approach to predict species distributions could help to ensure that uncertainty in future projections is better understood. As further analysis and standardisation approaches for fisheries survey data are developed, then uncertainty in both data and modelling methods should reduce, enabling projections to inform ecological and fisheries policy and management with increasing confidence.
Overall, this study suggests species will undergo distributional changes under increasing warming, with implications for community composition, although the capacity to track preferred temperatures may be impaired by wider environmental constraints not captured in our models (Edwards et al., 2021;Rutterford et al., 2015). Persistence of species with limited capacity to move will depend partly on their ability to acclimate to new thermal regimes. However, in the long term it may also require evolutionary adaptation to changing conditions, through selection on standing genetic variation, or via novel mutations Riebesell & Gattuso, 2015). For example, demersal species in the Northeast Atlantic may need to adapt to warmer or deeper conditions where there is no available habitat of suitable depth and temperature (Rutterford et al., 2015).
Under different climate change scenarios, the capacity for, and rate of, distributional change between species is likely to result in unpredictable changes to the fish assemblage. Efforts to limit and mitigate anthropogenic pressures exerted on marine systems will also influence the likelihood of sustaining populations (Duarte et al., 2020).
The Northeast Atlantic survey data used in this study could be further explored with detailed species-level analyses. Exploring the future distribution and abundance of individual species would provide tangible evidence to guide management of commercially important (Baudron et al., 2020) and threatened (Sguotti et al., 2016) species in a changing climate. Although this study reports on a subset of the full surveyed community, the relative responses of species could provide insights into the stability of interconnected assemblages under a changing climate. The inter-survey standardisation approach we use, however, precludes absolute abundance estimates, as needed for stock assessments, where greater detail about species catchability and gear variations would need to be incorporated, as well as information on size and age structures of populations (Fraser et al., 2007;Walker et al., 2017).
In conclusion, our study clearly identifies sea temperature as the primary driver of fish community composition across the Northeast Atlantic continental shelf. This is one of most productive regions of the world for marine fisheries, with annual landings in excess of 5.8 million tonnes and valued at USD 8.5 billion (Pauly et al., 2020). Our analysis implies that fisheries at higher latitudes, within the study area will be most impacted by further climate change. Some of the changes will likely result in opportunities for fisheries to exploit alternative species as warmer water species replace cooler water species (Simpson et al., 2011;Montero-Serra et al., 2015). We have not, however, projected diversity or functional group changes at specific locations. Nonetheless, there is potential to use our standardisation approaches to explore changes in species compositions across an extensive latitudinal range, and build on previous analyses of diversity change over smaller latitudinal ranges (e.g. Magurran et al., 2015).
The key finding from this study is that temperature is undoubtedly one of the most prominent drivers of community structure in marine demersal fish species of Northeast Atlantic shelf seas, and likely more broadly across latitudinal gradients. Warming seas over recent decades are already implicated in changes to distributions and abundance of species across this study area (Barange et al., 2018;Montero-Serra et al., 2015;Teixeira et al., 2014). With further predicted changes to fish communities, evaluating the effects of temperature on stocks may help to support appropriate management (e.g. Watson et al., 2022), although this may be challenging given that precise mechanisms by which temperatures affect stocks are typically not clear. Where species distributions are predicted to change, and potentially move between existing management areas, internationally coherent management approaches will be needed to enable sustainable harvesting and conservation under a changing climate.

ACK N O WLE D G E M ENTS
We thank the national survey teams and contributors to the International Council for the Exploration of the Sea (ICES) fisheries survey datasets. We also thank the Centre for Environment

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare that there is no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available at https://doi.org/10.5281/zenodo.7605008.