Disentangling the effects of fishing and temperature to explain increasing fish species richness in the North Sea

Climate change and fishing have impacted marine species richness (R) at global and local scales. It has previously been shown that R of the fish community in the North Sea has increased since the early 1980’s. Over the same period, ocean temperature has increased, and fishing mortality has decreased in the North Sea. Because these are confounded over time, either trend could plausibly explain the increase in R. Therefore, a logic-based approach was used to disentangle the effects of temperature and fishing on R, using spatio-temporal models fitted to survey data. To investigate the effect of temperature on R, fish species were subset by thermal affinity, as either Lusitanian (warm) or Boreal (cold) species. To investigate the effect of fishing mortality on R, species were subset by management category as either quota (assumed to be targeted) or non-quota species. Trends in these subsets were plotted separately to investigate which subsets of the fish community have contributed to the overall R increase. Over three decades, fish R increased by an average of 2.5 species per haul. These increases were predominantly of Lusitanian non-quota species (1.9). A small increase was observed in quota species (0.6); however, this increase was driven by quota-Lusitanian species (0.4). Our results suggest that temperature rather than fishing mortality was the driver of R increase in the North Sea since 1991 and highlight the importance of long-term monitoring in detecting ecological responses to climate change at the community level.


Introduction
Biodiversity is an important indicator of ecosystem health, impacting food webs and ecosystem productivity, as well as providing cultural, social and economic value to communities (Hooper et al. 2005;Cardinale et al. 2012). Global biodiversity is in decline and negatively correlated with anthropogenic pressures, such as climate change and habitat loss (Butchart et al. 2010;Tittensor et al. 2014). Projections of terrestrial, freshwater and marine scenarios of environmental change suggest that the decline in global biodiversity will continue over the 21st Century (Sala and Knowlton 2006;Pereira et al. 2010). Rates of species extinction have increased significantly since 1900, and this trend is unlikely to be reversed in the short term (Pereira et al. 2010;Batt et al. 2017;Chase et al. 2019).
Species richness (R), i.e., the number of species present at a given point in time for a defined area, is a simple and widely used metric for measuring biodiversity. Despite the well-documented decline in global R, at smaller spatial scales in the marine environment, R has increased over recent years (Dencker et al. 2017;Batt et al. 2017). Demographic processes such as extirpation and colonisation affect R at local scales but do not necessarily lead to changes in R globally (Sax et al. 2002;Chase et al. 2019). For example, a species newly colonising an area does not increase global R if it has previously existed elsewhere. Similarly, though a species may disappear from a region, it may survive elsewhere, therefore, global R would not change. If the number of new species expanding their distribution into an area is greater than the number leaving, then R will increase locally (Hiddink and Coleby 2012;Chase et al. 2019). This could be particularly true for marine systems, where connectivity may be higher than in terrestrial systems due to reduced barriers to marine dispersal (Carr et al. 2003;Osland et al. 2021).
Climate change and warming ocean temperatures are both likely to impact marine R (Hillebrand et al. 2018). Poleward shifts in species distributions have already occurred due to temperature increases (Dulvy et al. 2008;Hillebrand et al. 2018). Range expansions into waters with relatively lower temperature were more common than range contractions away from relatively warmer waters, suggesting that increases in regional R are more likely in the short term (Hillebrand et al. 2018). This is explained by the difference in speed at which colonisation and extirpation occur. Colonisation can occur and be observed quickly, particularly in marine environments where fish are good at tracking thermal conditions. In comparison, extirpation happens relatively slowly, as temperature increases are often not of a magnitude to cause sudden die-offs and competitive displacement by colonising species also takes time (Hillebrand et al. 2018). In temperate regions, where average sea surface temperature (SST) is below 20 °C, R is likely to increase due to poleward range expansion of warm-water species (Worm and Lotze 2016;Chaudhary et al. 2021). In eight out of nine North American sea regions, covering both the Atlantic and Pacific, trends in R between 1984 and 2014 were positive (Batt et al. 2017). Increasing temperatures could explain this increase as several of the species monitored within the study were sensitive to temperature change and expanded their ranges (Batt et al. 2017). A recent study on marine fish species in the North Atlantic between 1977 and 2013 found that temperature accounted for greater than 50% of the deviance in R explained in their models (Gislason et al. 2020). In the central Baltic and Kattegat, both the ecosystem-scale R and local R increased over a 10-year period (Hiddink and Coleby 2012). In the Kattegat this increase was linked to rising ocean temperatures (Hiddink and Coleby 2012).
Fishing is widely assumed to negatively impact biodiversity but can be complicated to interpret. Fishing can reduce biodiversity, including R, by altering size and community 1 3 structure of marine ecosystems, depleting commercial fish stocks and causing localized extinctions (Myers, Hutchings and Barrowman, 1997;Greenstreet and Rogers 2006;Thurstan and Roberts 2010;Niklitschek et al. 2010;Hiddink and Coleby 2012;Crépin et al. 2012;Eero et al. 2012;Wing and Jack 2013). Stock depletion may lead to reduced R at local levels due to sparser distribution of individuals resulting in a reduced likelihood of occurrence in sampling. Declines in large predatory species due to high fishing pressure are well documented (Myers and Worm 2003;Ferretti et al. 2008Ferretti et al. , 2010. For example, barndoor skate in the North-West Atlantic have seen sharp declines and they appear to have been extirpated from the majority of their range at < 1000 m depth where fishing occurs (Casey and Myers 1998). In the Baltic Sea, cod stocks have been overexploited and extirpated from parts of their range and, even after some levels of recovery, have not re-colonised their full historic range (Eero et al. 2012). Roundfish declines in some heavily fished areas have led to community composition changes and replacement by shellfish (Frank et al. 2005;Jackson 2008). Fishing can also cause declines in non-target species that are taken as bycatch, leading to reductions in abundance and R (Niklitschek et al. 2010). Conversely, through removal of large predatory fish or dominant species, fishing could potentially lead to increases in biodiversity as smaller species increase in number, known as meso-predator release (Hiddink and Coleby 2012;Ellingsen et al. 2015). It has also been argued that fishing does not cause direct extinctions of target species, since fishing becomes economically unviable before extinction (Dulvy et al. 2003;Burgess et al. 2013). However, this may not apply to multi-species fisheries given that several species contribute to maintaining the economic viability of the fishery (Dulvy et al. 2003;Burgess et al. 2013).
Surveys of fish abundance in the North Sea have been conducted annually since 1965, coordinated by the International Council for the Exploration of the Sea (ICES) and known as the North Sea International Bottom Trawl Survey (NS-IBTS) (ICES 2020). The NS-IBTS has been used to quantify spatial and temporal trends in abundance of bottom-dwelling (demersal) fish in the North Sea (Dencker et al. 2017). As all the fish species caught are identified to species level where possible, NS-IBTS data is ideal for quantifying indices of demersal fish biodiversity, including R. R is highest in the northern North Sea and along the coast of the United Kingdom (UK) and lowest in the German Bight (Dencker et al. 2017). Generally, R is expected to be greater in areas with greater benthic complexity and greater productivity which are typically greater in coastal areas (Tews et al. 2004;Koivisto and Westerbom 2010;Verdiell-Cubedo et al. 2013;Henseler et al. 2019). Therefore, it is not surprising that higher species richness was observed in coastal areas of the North Sea. Over the last 30 years, R in the North Sea has increased by an average of around 3 species per haul (Dencker et al. 2017). These increases are not spatially uniform and are concentrated in the Western coastal regions and entry points to the North Sea (Dencker et al. 2017).
There are several possible explanations for increasing R in the North Sea. The North Sea is one of the fastest warming areas with SST increasing by 1.31 °C between 1982 and 2006, the second highest change of any large marine ecosystem (Belkin 2009;Capuzzo et al. 2018). The winter bottom temperature for the North Sea increased by 1. 6 °C between 19806 °C between and 20046 °C between (Dulvy et al. 2008. Over the last 30 years, because both SST and sea bottom temperature (SBT) have increased further, both pelagic and demersal fish species will have experienced warming (Fig. 1a, b). Species associated with more southerly ranges have been increasingly found in the North Sea (Dulvy et al. 2008). Over the same period, fishing mortality has also declined (Fig. 1c). The average fishing mortality (Fbar) of the fully exploited age classes of several commercially important species (fish of ages preferentially 1 3 targeted by fisheries), has declined over the last 30 years due to the full implementation of sustainable management measures based on the principle of maximum sustainable yield (Mesnil 2012;ICES 2022). High fishing pressure in the North Sea coincided with a decline in multiple indices of Scottish groundfish biodiversity, including R, between 1925 and 1996 (Greenstreet and Rogers 2006). These decreases were largest in areas with high fishing effort suggesting that fishing was the main driving factor (Greenstreet and Rogers 2006). Improvements in fisheries management may have also led to increases in large fish community R between 2000 and 2010 (Greenstreet et al. 2011).
Though R is the simplest measure of biodiversity it can provide important detail on community changes. R is more sensitive to the early stages of species colonisation in comparison to abundance-based measures and is easily interpretable. Similarly, increases in R per haul from survey data could indicate a recovery of species impacted by fishing. Increasing abundance and expanding species range back into previously occupied areas also has the potential to increase R. Understanding the factors driving increase in R, previously reported by Dencker et al. (2017), is important for ecosystem-based management. Since both an increase in temperature and a decrease in fishing mortality are confounded in the North Sea ( Fig. 1), disentangling their separate effects on R using a purely statistical approach is not straightforward (McLean et al. 2019). The aim of this study is to disentangle the separate effects of temperature and fishing on R using a simple logic-based approach, comparing temporal trends in subsets of the demersal fish community in the North Sea corresponding to functional groupings of species. If temperature is the main driver of R increase in the North Sea (Fig. 1a, b), then it could be expected that an increase in Lusitanian (warm) species will be detected over the time period that warming has occurred. If declining fishing mortality is the main driver of R increases in the North Sea (Fig. 1c), then it is expected that increases will be detected principally in quota fish species which are now subject to reduced fishing mortality relative to the beginning of the study period, though this does not preclude synchronous increases in non-quota fish species due to reduced fishing effort. If increases are seen across all subsets of the fish community, then it is likely that both factors Fig. 1 Showing a sea surface temperature (SST), b, sea bottom temperature (SBT) and c fishing mortality (Fbar). SST data was taken from the Hadley ISST1 database (Rayner et al. 2003) and averaged for ICES Division 4a-c fitted with a simple GAM using a spline smoother. SBT data was taken from the ICES DATRAS data portal for the NS-IBTS and a random sample of 250 points was taken from each year then fitted with a GAM using a spline smoother. Fbar data was taken from the ICES WGNSSK report, the latest version of which was published in 2022 (ICES 2022) calculated as species-specific F bar for multiple species then fitted with a GAM using a spline smoother. Species included were Gadus morhua, Melanogrammus aeglefinus, Microstomus kitt, Pleuronectes platessa, Pollachius virens, Solea solea, Scophthalmus maximus, Merlangius merlangus, Glyptocephalus cynoglossus and Scophthalmus rhombus contribute to increases in R. Lastly, the spatio-temporal trends of R across these different subsets of the fish community were investigated. If the fish community in the North Sea followed poleward trends previously reported in other regions, it is expected that both the Boreal and the Lusitanian community will have shifted northwards. The temporal, spatial and spatio-temporal trends in the different subsets of R were then interpreted with regards to conservation of fish biodiversity in a changing North Sea ecosystem and the potential implications for fisheries management.

Data
Data were obtained from the NS-IBTS for quarter 1 from 1991 to 2019. The survey study area covers ICES divisions 3a, 4a-c and part of 7d, also referred to as ICES roundfish areas 1-10 ( Fig. 2). ICES rectangles within division 7d (roundfish area 10) were removed from the dataset as they were not sampled pre-2007. The quarter 1 survey is conducted annually between January and March by several different participating countries using standardised survey methods(ICES 2020). Though a similar survey in quarter 3 has been conducted annually since 1991, it was not included here due to the sparse and variable spatial coverage in the early years of the study period. Compiled fish abundance data were downloaded in February 2020 from the ICES data portal (https:// www. ices. dk/ data/ data-porta ls/ Pages/ DATRAS. aspx) in the format catch per unit effort per length per haul per hour for areas 1-10. Since 1991, only the Grand Ouverture Verticale (GOV) trawl were used by the IBTS to standardise for gear type. The NS-IBTS survey area is divided into ICES statistical rectangles of 1° longitude and 0.5° latitude. Countries involved in conducting surveys are allocated several rectangles, aiming to survey each rectangle at least twice per year in quarter 1 by two different countries, though this is not always achieved. Only ICES rectangles which were surveyed in at least 23 out of 29 years were included (80% of the study period). This left a total of 11,012 hauls from the study period. From this dataset the variables year, latitude, longitude, number of fish at length per species were used. All non-fish species were removed from the dataset reducing the dataset to only fish species including elasmobranchs. In 2018 the maximum sampling depth in the North Sea was changed from 200 to 250 m to reflect that a number of hauls before 2018 already exceeded this maximum depth. However, this is not thought to have impacted the results as there were never more than five hauls deeper than 200 m in a given year with a median of three.
Though the NS-IBTS aims to identify all individuals to species level where possible, in some hauls individuals are only identified to genus or family level. Values of R have the potential to be higher in later years if taxonomic resolution improved over the study period i.e., individuals only identified to genus or family level early in the study period were identified to species-level in more recent years. Therefore, observations of individuals to genus/family-level, and to species-level within each family were compared over the study period. Only Ammodytes sp. and Syngnathus sp. showed potential improvements in taxonomic resolution. To be consistent, these were aggregated to genus level throughout and treated as a single species. Species from three other genera: Alosa sp., Callionymus sp. and Pomatoschistus sp. were aggregated due to uncertainty around species level identification. Four species were removed from the dataset due to likely misidentification, for example, native ranges confined to the Americas or due to not being accepted distinct species 1 3 (supplementary material). Other minor adjustments to the dataset as downloaded were also made and are detailed within the supplementary material. This left a total of 162 species identified in the dataset. Total R per haul of the fish community (R t ) was then calculated using this dataset.

Biodiversity subsets and indices
To apply our logic-based approach to disentangle the effects of ocean warming and fishing on fish R, species were subset by thermal affinity and management category. To investigate the effect of ocean warming, species richness was calculated for either Lusitanian, R l , Boreal, R b , or Atlantic species as per Yang (1982). Species not included in Yang (1982) were categorised using their native ranges as reported via FishBase, a global database on finfish (https:// www. fishb ase. se/). This resulted in 60 species categorised as Boreal, 88 1 3 categorised as Lusitanian and 14 categorised as Atlantic within the dataset. Atlantic species richness made up a very low proportion of R t per haul in the North Sea, only present in 4129 of 11,012 hauls, 37.5%, and showed no significant change over the study period. Therefore, the results for the Atlantic species are not included for the thermal affinity comparison, though Atlantic species were included in the management subsets when present.
To investigate the effect of fishing mortality, species richness was calculated for each management subset, i.e., whether they were quota-regulated species in the North Sea (R q ) or not (R n ). For the purpose of this study, species were classified as quota-regulated if they were assigned species-specific TAC, or TAC shared across two species for ICES divisions 3.a or 4.a-c in the EU council regulation on fishing opportunities for 2020 (EU 2020). This included TAC set for the Lophiidae family of which there are two species within the North Sea. Though there is TAC set for the order Rajiformes, these were not included within the quota subset due to the breadth of species this TAC is spread across. Though some species found during the North Sea IBTS are allocated TACs for areas outside ICES divisions of 3.a and 4.a-c, these were not considered as 'quota' species in this study. This subset was used for identifying the impacts of fishing as species which are assigned species specific quota would show the best signal for reductions in fishing mortality in the North Sea. Twenty-five species were categorised as quota species for the North Sea leaving 137 nonquota species. Quota species richness was further subset according to thermal affinity as quota-Lusitanian (R ql ) or quota-Boreal (R qb ) to analyse whether potential changes in the quota community due to climate change (Lusitanian, n = 11, Boreal, n = 12, and Atlantic n = 2, Table 1). Similarly, non-quota species were subset by Lusitanian (R nl n = 77) and Boreal (R nb n = 48).

Analyses
Species richness for each haul was computed from the NS-IBTS abundance data to give R i,h , where i is the subset of data modelled (t, l, b, q, n, ql, qb, nl, nb) and h is a given haul. Generalised additive models (GAMs) were used to represent the general underlying spatial and temporal trends in R i,h over the time series. GAMs were used as they apply smoothers to the data, allowing for a non-linear relationship between time, space and R i,h . Four simple GAMs were compared: with either time (Year; Eq. 1) or space (denoted by latitude, Lat, and longitude, Long; Eq. 2) as the sole predictor of R i,h a spatio-temporal model allowing for change in both space and time but with changes in space being constant in time (Eq. 3) or a spatio-temporal model allowing for simultaneous changes in both time and space (Eq. 4). A thin plate regression (s) spline was used as the smoother for the temporal model and a tensor (te) product was used to create the models with a spatial component. All models were fitted using the restricted maximum likelihood method assuming a Poisson distribution. These three models were used to calculate a predicted value of R ( R i,h ) are given by: where i,h is an independent and identically distributed (i.i.d) random error. AIC and % deviance explained of each model for each subset was then compared (Akaike 1987;Burnham and Anderson 2004). Predicted R per haul for the midpoint of each ICES rectangle included in the dataset and for all years in the time series (1991-2019), R i,h,rec,y , and 95% confidence intervals (hereafter 95% CI) were computed using the predict.gam function where rec denotes a given ICES statistical rectangle and y denotes year. To represent the spatial distribution of R i,h,rec in the North Sea, we calculated the mean R i,h,rec,y by subset across the whole time series ( R i,rec ) as: To illustrate the overall temporal trend and compare between subsets, the mean species richness per subset for the whole of the North Sea for each year ( R i,y ), was calculated as: as well as the 95% CI. 95% CI were used to infer significant change in R i,y , i.e., where 95% CI shows no overlap between the start and end of the study period, the change is considered to be significant (p < 0.05).
The difference in R i,y between 1991 and 2019 was calculated to illustrate the magnitude of change in R i,y over the entire study period was calculated as: h,rec,1991 was the R per haul for a given ICES rectangle in 1991 and n h,rec,1991 was the number of rectangles in that year. R i,h,rec,2019 was the R per haul for a given ICES rectangle in 2019 and n h,rec,2019 the number of rectangles in that year. R i,1991 was then subtracted from R i,2019 to estimate the difference in R i between the start and the end of the study period for the North Sea as a whole, R i , using 95% CI to specify whether changes were significant.
Spatio-temporal change in R i,h,rec between 1991 and 2019 was calculated using the raster calculator function in R. R i,h,rec,1991 and R i,h,rec,2019 , the predicted R for each rectangle in the years 1991 and 2019, were used to create rasters. The 1991 raster was then subtracted from the 2019 raster to illustrate change in predicted R per rectangle, R i,h,rec , between the start and end of the study period. This method was chosen to illustrate the total change between the start of the study period to the end rather than rate of change as this varied between years.

Model selection
For all subsets, spatio-temporal models fit best (Table 2), with time explaining around 5% of the deviance, and space explaining around 40% of the deviance. For the majority of subsets Eq. 4 had the lowest AIC value and explained a greater percentage of deviance, though in some cases Eq. 3 had the lower AIC (Table 2). However, model Eq. 4 was selected for all subsets for consistency, as it had lower AIC in a greater number of models (Akaike 1987;Burnham and Anderson 2004).

Spatial variation of R
Mean fish community richness per rectangle, R t,rec , was highest in the North West around the Shetland islands, in the Skagerrak and Kattegat, and along the UK coastline   (Fig. 3a). R t,rec was lowest in the central North Sea and along the coast of continental Europe. Mean fish community richness per rectangle of Lusitanian species R l,rec was highest in ICES rectangles close to the coast, particularly those nearest the English Channel and Orkney (Fig. 3b). For boreal species, R b,rec was highest in the Skagerrak and Kattegat and showed a generally northerly distribution (Fig. 3c). For quota species, R q,rec was greater in the Northern North Sea and moderately high in the Skagerrak and Kattegat, whereas R n,rec was higher in the South West, UK coastal areas and the Skagerrak and Kattegat (Fig. 3d, e). R ql,rec was highest in the Skagerrak and Kattegat and distributed similarly to R l,rec , with the lowest values in the central and East (Fig. 3f). R qb,rec was greatest in the Northern North Sea and particularly low in the Southern North Sea (Fig. 3g). R nl,rec was distributed similarly to R l,rec , greatest at the Western entry points to the North Sea, particularly near the English Channel, and the coastal areas (Fig. 3h). R nb,rec was greatest in the Skagerrak and Kattegat and in coastal regions and was low in the central North Sea away from the coasts, and particularly low in the Northernmost areas of the North Sea (Fig. 3i). Patterns in R t,rec were similar to those observed in R l,rec , R n,rec , R ql,rec and R nl,rec . Distributions of R b,rec , R q,rec and R qb,rec were all similar to each other but distinct from R t,rec .

Fig. 3
Spatial distributions of average R using model predictions for the different subsets of the fish community, a R t , the whole dataset, b R l , Lusitanian, c R b , Boreal, d R q , quota, e R n , non-quota, f R ql , quota-Lusitanian, g R qb , quota-Boreal h R nl , non-quota-Lusitanian i R nb , non-quota-Boreal, with minimum and maximum values on the colour scale differing between subsets Temporal and spatio-temporal trendsR t,y increased significantly by 2.5 (18.1%) species between 1991 and 2019 (Fig. 4a). The values of R t,rec describing changes over the study period were not uniform in the North Sea (Fig. 4b). The central North Sea showed the least change, indicated by low R t,rec , with a slight decrease at the very centre. Large values of R t,rec particularly in the Northern North Sea to the east of Shetland and off the South East coast of the UK near the English Channel indicate an increase in total species richness (Fig. 4b). Increases were also observed across all coastal areas of the North Sea.

Temporal and spatio-temporal trends in the thermal affinity subset
Subsetting by thermal affinity found a significant increase in R l,y of 2.4 (43.6%) species per haul over the study period, whereas R b,y declined by − 0.4 species per haul. There was clear overlap in the 95% CI for R b,y shown between the 1991 and 2019, therefore it cannot be concluded that R b,y has decreased significantly (Fig. 4c). These trends mean that Lusitanian species have become an equal component of R t,y as Boreal species in the North Sea since the mid 2000s. Spatio-temporally, values of R l,rec showed increases throughout the North Sea over the study period, though increases were highest off the South East coast of the UK and in the Skagerrak and Kattegat (Fig. 4d). R b,rec showed slight increases in the North-Eastern North Sea and in the Southern North Sea (Fig. 4e). Decreases were found in the central North Sea of greater than 2 species per haul and of a similar magnitude in the Skagerrak and Kattegat.

Temporal and spatio-temporal trends in the management subset
Subsetting by management found a slight increase in R q,y over the study period of 0.6 species per haul (Fig. 4f), but because there was overlap in the 95% CI of 1991 and 2019 these changes cannot be considered significant. R n,y accounted for a greater proportion of the change in R t,y in the North Sea, significantly increasing by 1.8 (28.5%) species per haul (Fig. 4f). Spatio-temporal trends showed a clear spatial difference in R q,rec within the North Sea (Fig. 4g). R q,rec is high in the northern North Sea, with much lower changes observed in the Southern North Sea. R n,rec showed similar patterns to those seen in R t,rec , though, the greatest increases were found off the Southeast coast of the UK with smaller increases in the Northern North Sea and the coastal regions (Fig. 4h).

Temporal and spatio-temporal trends in the management-thermal affinity subset
Subsetting by quota species by thermal affinity showed an increase in R ql,y of 0.4 (16.0%) while R qb,y decreased by 0.1 (Fig. 4i). The increase in R ql,y was significant, whereas there was overlap in the 95% CIs of R qb,y . Splitting non-quota species results in R nl,y increasing significantly by 1.9 (63.3%) and R nb,y decreased by 0.2 (Fig. 4l), though not significantly. Spatio-temporally, R ql,rec showed slight increases in the northern North Sea and the Skagerrak and Kattegat, similar to those seen in R q,rec (Fig. 4j). R qb,rec showed slight increases in the Northern North Sea and slight decreases in the central North Sea, with little change occurring in the South (Fig. 4k). R nl,rec showed very similar patterns to R n,rec with increases seen throughout the North Sea, though the greatest increases occurred in the South East (Fig. 4m). R nb,rec showed slight increases in the South East, slight decreases in the central and North-Western North Sea, and little change in the north-eastern North Sea (Fig. 4n).

Temporal trends in the thermal affinity subset
Overall, R of the fish community increased between 1991 and 2019, confirming an earlier study (Dencker et al. 2017). Increases of R in the North Sea could be interpreted as positive sign of recovery, however, increases in R were driven by increases in Lusitanian species with an affinity for warmer waters. The increase in Lusitanian species in the North Sea since 1991 means that they make up an equal proportion of R in the North Sea fish community as Boreal species. Therefore, we conclude that recent increases in R in the North Sea were driven by increasing temperatures, rather than decreases in fishing mortality. Past research has shown the ratio of Lusitanian to Boreal species R increased during periods of higher sea temperature in a small section of the Southern North Sea (Ter Hofstede and Rijnsdorp 2011). Our results show that these increases have been observed throughout the North Sea and at the individual-haul level which illustrate the large magnitude of the changes occurring in the fish community.
Recent increases in R suggest a fish community response to climate change through an increase in Lusitanian species. These results are consistent with the prediction that temperate sea regions with an average SST below 20 °C will see increases in R due to ocean warming (Chaudhary et al. 2021). Changes in the composition of the fish community towards warmer water species have previously been referred to as tropicalization (Vergés et al. 2014). Originally this term was used to mean shifts from temperate fish communities to species typically associated with tropical environments. However, in some cases it has been used more generally to refer to shifts towards species typically associated with warmer waters, such as observed in the North Sea (McLean et al. 2021). These shifts have also previously been referred to as sub-tropicalization in the North Sea and Baltic Sea (Montero-Serra,Edwards and Genner, 2015). Both terms may misrepresent the situation in the North Sea, where many Lusitanian species may already occur in the southern North Sea and have expanded their range, or occurred rarely but are becoming increasingly common. Fish community changes towards species typically associated with warmer waters have already been observed in Australia, where this has led to food web changes with increased herbivory and damage to kelp forests (Wernberg et al. 2016;Vergés et al. 2016). There is comparatively less research on how shifts towards warmer thermal affinities will affect ecosystem functioning in cold-temperate regions such as the North Sea, though it could be expected that food web and energy pathways may be altered, as has been predicted for the Arctic (Frainer et al. 2017;Friedland et al. 2020). Increases have also been found in the community temperature index (CTI), which represents the abundance-weighted mean thermal affinity of Fig. 4 Temporal trends in R of a the fish community, c thermal affinity subsets, f management subsets, i quota-thermal affinity subsets and j non-quota-thermal affinity subsets. Changes in R between 1991 and 2019 at the spatio-temporal level of b the fish community, d Lusitanian, e Boreal, g quota, h non-quota, j quota-Lusitanian, k quota-Boreal, m non-quota-Lusitanian and n non-quota-Boreal 1 3 the community. Fish communities in the North Sea, US coastal regions and the Barents Sea have all seen increases in their CTI (McLean et al. 2021). Increases in the CTI have generally been assumed to be driven primarily by increases in abundance of species with warmer thermal affinities. However, in the Barents and Bering Seas increases in the CTI were driven by deborealization, i.e., the loss of Boreal species (McLean et al. 2021). Our results are consistent with the interpretation that increases in Lusitanian species have had a greater impact on species richness in the North Sea than deborealization, since increases were predominantly seen in Lusitanian species with little decline observed in Boreal species. This may be because current warming has not driven current temperatures above the thermal tolerance of Boreal species yet. Although deborealization was not detected in the North Sea, Boreal species have been shown to move to colder and deeper waters in previous studies (Dulvy et al. 2008). Eventually with further temperature increase, Boreal species will be constrained by lack of continuously deeper, cooler waters in the North Sea leading to a risk of future declines in Boreal species (Rutterford et al. 2015;McLean et al. 2021).

Temporal trends in the management subset
While there was an increase in quota R, increases in non-quota R made up a greater proportion of the increase in total R. Furthermore, increases in quota R resulted from increases in the Lusitanian component of the quota community, giving further evidence that increases in temperature are the principal drivers of the observed increases in R. Conversely, the Boreal quota community showed slight declines, which is consistent with expected trends under increasing temperatures. Though increases in Lusitanian quota richness were small in magnitude they may provide the first indication of future trends. Trends in non-quota R, when subset further by thermal affinity, were also driven by increases in the Lusitanian component, though in this case the increases were significant. Since both the non-quota and quota community showed similar trends, where increases were observed in Lusitanian species but not in Boreal, it is unlikely that these trends represent a change in fisheries selectivity of non-quota species. Therefore, the observed change of R in the North Sea are primarily due to increasing R of Lusitanian non-quota species which are adapted for warmer temperatures in the North Sea brought about by climate change, rather than improved fisheries management and reduced fishing mortality.
Increases in Lusitanian species, particularly Lusitanian-quota species, have implications for both fisheries management and industry. The current system of quota allocation in the North Sea uses the concept of relative stability, assigning quotas for commercial species to countries based on historical catch levels between 1973 and 1978 (Morin 2000). Our results show that the fish community has substantially changed from that in 1991, and by extrapolation from the fish community in 1973-1978, due to increases of Lusitanian species richness. As such, fishers are likely to encounter Lusitanian species for which they have limited quota more regularly. Lusitanian species for which quota are already allocated using relative stability, such as hake (Merluccius merluccius), may increase in abundance creating issues for fisheries if catches also increase (Baudron and Fernandes 2015). Zonal attachment, the proportion of fish stock biomass present in a country's Exclusive Economic Zone, is already out of sync with the UK's quota allocation in the North Sea (Fernandes and Fallon 2020). Further shifts in community composition could exacerbate this issue and have economic consequences for the fishing industry. Equally, future deborealization could lead to more restrictive TAC for boreal species which will also have implications for the fishing industry. These changes in the North Sea fish community may necessitate changes to the policy of relative stability as fish distributions change. This would mean changes to the setting and allocation of quotas to reflect current fish distributions, potentially along the lines of zonal attachment (Fernandes and Fallon 2020). Any changes to quota allocation however, are likely to be politically complicated.

Spatial and spatio-temporal trends in R subsets
Increases in Lusitanian species were found throughout the North Sea, with only a small area in the central North Sea showing little increase. This differs from a more simplistic prediction of a poleward gradient in increases in Lusitanian R, with greater increases in the South compared to the North. However, water currents into the North Sea from the Atlantic bring more saline water in between Orkney and Shetland, which may also facilitate increases in Lusitanian R in the northern North Sea alongside recent temperature increases (Mathis et al. 2015;Tian et al. 2016;Quante et al. 2016). The connectivity to adjacent ICES areas 6 and 7, which exhibit greater taxon richness than the North Sea, may also contribute to increases in richness through colonisation from these areas (Heessen et al. 2015). Spatial distribution of Lusitanian R showed the highest values at the Western entry points to the North Sea and in coastal regions. Areas where Lusitanian species were typically lower in number may also be expected to show larger increases in the future as Lusitanian species already present in the North Sea begin to colonise further. Therefore, it is not surprising that extent of spatio-temporal changes in Lusitanian R can differ from the distribution of Lusitanian R. Boreal species showed a slight increase over the study period in the Southern region of the North Sea, contrary to the expected poleward shift. This increase was driven by non-quota Boreal species and may be due to other environmental or anthropogenic factors such as an increase in offshore windfarms, providing habitat for fish species, or reduced competition by quota-Boreal species which shifted their distributions.
Slight increases in both components of the quota community were found in the Northern North Sea. This may be due to poleward movements of Boreal quota species, as has been observed in cod (Engelhard et al. 2014), and Lusitanian species entering through the Northern North Sea. Interestingly, increases in quota R in the Northern North Sea coincide with areas with the highest demersal landings (fish caught which are subsequently taken to port to be sold) by the UK fishing fleet (UK Sea Fisheries Statistics, Marine Management Organisation 2020). When total fishing effort for all nations fishing in the North Sea are considered, the picture is more complex. Otter trawl effort is fairly spread out while beam trawl effort is concentrated in the Southern North Sea suggesting total fishing effort is higher in the Southern North Sea (Couce et al. 2020). Spatio-temporal changes in quota R are not consistent with this pattern in fishing, particularly that areas where the largest decreases in quota-Boreal R have been seen are areas where fishing disturbance has been lowest (OSPAR 2017). This is consistent with our findings that reductions in fishing mortality are unlikely to have been the cause of recent increases in R.

Usefulness and future research
This study highlights the utility in simplistic, functional subsets to disentangle different factors likely to drive increases. It compliments previous modelling approaches which have aimed to quantify the impacts warming and fisheries management as well as predicting under future climate scenarios (Serpetti et al. 2017;Beaugrand et al. 2022). Our results are consistent 1 3 with research in the Celtic Sea which found that the environment likely had greater impact on demersal community structure than changes in fishing (Mérillet et al. 2020). Though it could be argued that R is only one simple measure of biodiversity, it has the advantage of being particularly sensitive to colonisation, a key component of climate change effects on fish communities. It is also easily interpretable as a measure of community change where other metrics, such as evenness, combine signals which are more difficult to interpret in the context of climate change. Our results showing the shift in species richness towards increasing levels of Lusitanian species can be considered an early indicator of the effects of climate change on the fish community.
However, other metrics are required to further quantify the impacts of climate change on the fish community within the North Sea. For example, a previous study showed the biomass of Boreal species in the North Sea was greater than that of Lusitanian species (Yang 1982). With recent increases in Lusitanian R in the North Sea, it is important to understand whether similar trends have occurred in biomass. Since this research is now outdated, future research to quantify whether biomass of Lusitanian species is increasing or even exceeds Boreal species is necessary to determine whether Lusitanian species are the dominant thermal affinity group within the fish community in terms of abundance. This has also been highlighted for the Celtic Sea where cold-water species may be susceptible to biomass loss and northward shifts in range due to increasing sea bottom temperature (Mérillet et al. 2020). Conversely, warm-water species have become more numerous within the Celtic Sea and could become the dominant species group if shifts continue (Mérillet et al. 2020). In the Eastern Mediterranean, increases in the abundance of warm affinity species and decreases in the abundance of cold affinity species have been observed, showing a strong reaction to climate change (Givan et al. 2018). Abundance measures are likely to be more sensitive to declines in Boreal species, as species will become rarer before they are extirpated due to temperature increases. Therefore, quantifying changes in biomass is important for measuring further implications of climate change on the North Sea fish community.
The large fish indicator (LFI), the proportion of fish over a certain size in the community, is another measure of fish biodiversity and an Ecological Quality Objective in the North Sea (Greenstreet et al. 2011). The impact of increasing Lusitanian species richness on the LFI is not yet understood, though an increase in the proportion of Lusitanian species in the North Sea, which are typically smaller, could lead to an overall decline in LFI . Other measures such as Typical Length and Mean Maximum Length which have also been suggested for monitoring Good Environmental Status may also differ between the two thermal affinities (Lynam and Rossberg 2017). Increasing temperatures may also reduce the overall size of the Boreal community which may have implications for reproductive output and population stability (Baudron et al. 2014;Tu et al. 2018). Further research could look to use the logical community subset approach for investigating how various size-based indicators have over a period characterised by changes in temperature and fishing mortality, and whether previous targets will be achievable in a changing North Sea. This could help to illustrate potential future impacts of increasing Lusitanian species richness on the LFI in the North Sea.

Conclusion
By using a simple, logic-based approach to investigate trends in R of different subsets of the community, it is possible to identify key drivers of community change. These simple subset methods are useful where explanatory variables are confounded and can provide 1 3 further detail on changes within the fish community. While R has increased significantly in the North Sea over the last 30 years, these increases are largely driven by increases in the number of Lusitanian species. Increases found in quota R were mostly due to increases in Lusitanian species. Both results suggest that climate change is the main factor behind recent increases in R in the North Sea. Spatio-temporal patterns of R are more complicated but show a widespread increase in Lusitanian species richness and some suggestion of a poleward shift in quota species, particularly Boreal quota species.