Three decades of increasing fish biodiversity across the north-east Atlantic and the Arctic Ocean

Observed range shifts of numerous species support predictions of climate change models that species will shift their distribution northwards into the Arctic and sub-Arctic seas due to ocean warming. However, how this is affecting overall species richness is unclear. Here we analyse scientific research trawl surveys from the North Sea to the Arctic Ocean collected from 1994 to 2020, including 193 fish species. We found that demersal fish species richness at the local scale has doubled in some Arctic regions, including the Barents Sea, and increased at a lower rate at adjacent regions in the last three decades, followed by an increase in species richness and turnover at a regional scale. These changes in biodiversity paralleled an increase in sea bottom temperature. Within the study area, Arctic species’ probability of occurrence generally declined over time. However, the increase of species from southern latitudes, together with an increase of some Arctic species, ultimately led to an enrichment of the Arctic and sub-Arctic marine fauna due to increasing water temperature consistent with climate change. Significance Statement Global modelling studies suggest increased species arrivals from lower latitudes and local expirations at high latitudes due to global warming. Our analysis of 20,670 standardized scientific trawl surveys with 193 fish species from the north-east Atlantic and Arctic Oceans found an increase in species richness in the region parallel to an increase in sea bottom temperature. Some Arctic species declined in probability of occurrence over time, but some increased. This, together with the increase of southern-latitude species led to an enrichment of the Arctic and sub-Arctic marine fauna attributed to climate change.


Figures captions
Climate warming constitutes one of the main faces of climate change, and is having a direct impact on 71 species, communities and ecosystems. In the oceans, the average increase in temperature in the last 72 140 years has been of 1°C (1, 2). Marine ectotherms occupy most of their potential latitudinal range with 73 regard to thermal tolerance and therefore move to higher latitudes following the displacement of their 74 thermal niche (3,4). However, these changes can occur at different rates across species, depending on 75 traits such as their dispersal potential, thermal niche and capacity to exploit new resources, which may 76 lead to changing community composition (5, 6). Understanding how these changes occur is crucial for 77 effective conservation and management strategies. Yet, to date, consistent empirical evidence of a 78 generalization of these shifts in the Arctic fish community is lacking.

80
Arctic and sub-Arctic ecosystems are among the most rapidly warming regions in the world, with some 81 areas warming four times faster than, and seas warming twice, the global average (2,7,8), and their 82 species composition may be changing accordingly (9). Until now, a doubling in species richness have 83 been reported in some areas of the North Sea, though not in others, and increases of smaller magnitude 84 have also been reported around North America (10-12). However, fish community analyses are mostly 85 restricted to non-polar latitudes, and they rarely exceed the 62 °N of the north Bering's Sea. Studies 86 northern than this only exist in the Barents Sea, where some boreal species arrived recently (9,13,14) 87 (Table S1). Of these, one examined distributional shifts in a fish community of 49 species from 2004 to 88 2017, and it reported an increase of less than 50% in species richness (14). This represents the only 89 work reporting empirical evidence of a regional increase in fish species richness with climate change in 90 Arctic latitudes. Thus, although an increase in species richness is predicted into the Arctic Ocean, and 91 several model projections exist in the literature, including species extirpations (15-17), the empirical 92 evidence is limited temporally and taxonomically, and lacks a long-term correlation with climate warming.

94
In the Norwegian and Barents Seas, recent warming and increased Atlantic water inflow events have 95 been recorded, with a decline in sea-ice cover in the northern Barents Sea (18-20). As a consequence, 96 profound effects on the geographical distribution and productivity of commercial fishing stocks are 97 expected (21-23). In fact, species turnover was projected to increase in the area in the next decades, 98 resulting from some local species extirpations and the arrival of warmer-water species (16,24,25).

99
Here, we report three decades of field data to test these predictions.

122
Alpha diversity 123 A total of 193 demersal fish species were recorded between 1994 -2020 across the whole study area 124 ( Fig. 1, Table S2). The generalized additive models (GAM) used to account for the differences in 125 sampling effort across time, predicted a 68 % increase in the average number of species per trawl (alpha 126 diversity), from 8.1 species/trawl in 1994 (95 % CI 7.9, 8.3), to 13.6 species/trawl in 2020 (95% CI 13.3, 127 13.9) (Fig. 2, Table 1). The increase in species richness was correlated with changes in sea bottom 128 temperature (Pearson r = 0.59, p< 0.05, Supplementary Results). Increasing species richness was 129 also found for the adjacent areas in the North Sea and around Svalbard though no significant correlation 130 was detected with SBT at those regions (Table 1, Fig. S1).

134
Among the 17 explanatory variables included in the BRT model, "Depth" was the most relevant, followed 135 by Year and "Bottom temperature" (Fig. S2). The model predicted an increase in species richness at 136 higher latitudes, especially high in the Arctic region of the Barents Sea, where species richness more 137 than doubled in some regions, during the study period. (Fig. 3, Fig. S3).

139
Gamma diversity

140
The overall species richness in the main study area (gamma diversity), increased 45% during the study 141 period, from 65 species in 1994 to 94 species in 2020, following a similar temporal pattern to alpha 142 diversity (Fig. 4, Table 2, Fig. S4). Very similar increases were reported in adjacent areas (47 % in the 143 North Sea, 45 % in Svalbard, Table 2, Fig. S4, Fig. S4). GAM models for total species richness at each

149
Beta diversity 150 Pairwise mean total beta diversity increased from 1994 to 2002 in the main area, slightly declined until 151 2012, and slightly increased again until 2020, with an overall non-significant increase of 4 % (95 CI -3, 152 9) (Fig. 5). The turnover contribution to beta diversity, which is not affected by species richness, 153 increased until 2008, declined until 2014, and increased afterwards, with an overall increase of 16 % 154 (95 CI 6, 27) (Fig.5, Table 3). Similarly, total beta diversity in the adjacent region of Svalbard increased 155 until 2005 and declined afterwards. However, the turnover increased linearly across time in this region.

157
In the North Sea, both the total beta diversity and turnover declined with time (Table 3, Fig. S6). The

158
GAM models for temporal change in beta diversity and turnover, did not select trawl-swept area as a 159 relevant explanatory for the main study area or the adjacent regions. Thus, when the increasing species 160 richness is accounted for, species turnover beta diversity, as well as alpha and gamma diversity 161 increased over the years.

163
Species' trends 164 From the 193 species in our study, 99 species showed a significant temporal trend (increasing or 165 decreasing) in their probability of occurrence with time, in at least one of the studied areas (Fig. 6). Of 166 these, 71 species showed only positive trends, and 23 species showed only negative trends in at least 167 one studied area, while 5 species showed positive and negative trends depending on the study area significant consistent trend, increased. The number of species increasing was consistently higher than 170 the number of species decreasing over time, across all regions. Species declining mostly inhabited high 171 mean latitudes, and linear regression identified a negative effect of mean latitude on the temporal 172 change in species probability of occurrence in the main study area, and around Svalbard (Fig. 6). Among 173 the 67 species only found in the main region and/or Svalbard (Arctic region), 18 species showed a 174 significant change in their probability of occurrence with time, 6 increased and 12 declined.

178
Following global trends in warming, sea bottom temperature data suggests an increase of 0.3 °C per 179 decade in our study area from 1993 to 2019, and this increase was several times higher in some regions  thousands of species had shifted from equatorial to mid latitudes, and using independent data not 197 available to that study we confirm that this shift has continued into higher northern latitudes, leading to 198 two times more demersal fish species in the study area. Increases in alpha diversity were already 199 reported in some regions of the Arctic Ocean, albeit smaller areas and/or shorter time periods (14, 31).

200
With this study, we show that the increase in species richness is not localized, but widespread across 201 the region and concomitant with the increase in bottom temperature.

203
We identified sea bottom temperature as the most relevant temporal environmental variable in our

218
diversity is considered it has most commonly been measured as Jaccard's coefficient which is not 220 independent of species richness (48). Here, we find that while both alpha and gamma diversity are 221 significantly increasing, so is species-richness independent species turnover, but not total beta diversity 222 (Fig. 5). Thus, the spatial heterogeneity of biodiversity is increasing as well as biodiversity overall.

223
However, the difference between alpha and gamma diversity changes in our study is likely not large 224 enough to reflect a change in overall beta diversity.

226
Overall the study area, we found that 76% of the species showing a consistent significant change in 227 probability of occurrence increased with time. However, when restricted to species absent from the North

245
Study area

247
The database here analyzed comprises research trawl surveys restricted to the continental shelf and 248 slope from the north of the North Sea into the Arctic Ocean, mostly from 56 ° N to 81 ° N Latitude, and 249 from 2 °W to 50 °E Longitude (51). The whole area comprises a marked temperature gradient, with 250 average bottom water temperatures over 8 °C in the North Sea to -1 °C in the northern region of the 251 Barents Sea. Because of data temporal and spatial distribution, we considered one main study area, 252 with the longest time period (Norwegian & Barents Sea) and two adjacent areas (Svalbard and North 253 Sea), which contain a shorter temporal period (Fig. 1).

255
The database initially contained 60,355 research surveys, from 1980 to 2020. We excluded data which 256 were associated with broken gear, had incomplete metadata (data without reporting depth, or 257 coordinates, or type of gear), and were questionable (e.g., shrimp trawl opening of several km). We 258 restricted the analysis to shrimp trawling data using 20 mm mesh size, a maximum of 5 km trawling 259 distance and 60 m of trawl opening, from 30 m to 700 m depth. We only included fish species (classes 260 Actinopterygii, Elasmobranchii, Holocephali, Myxini and Petromyzonti). Invertebrates were not included 261 in the study. After data standardization and selection, 20,670 surveys from 1994 to 2020 remained.

263
Depending on the data spatial coverage over time, and based on the biogeographical realms that have 264 been described in the area (52), we divided the study area into a main area containing most of the data 265 and the longest temporal coverage, the Norwegian-Barents Sea, and two adjacent areas with more 266 limited data: the area around Svalbard and the North Sea (Fig. 1). Although the latitudinal and area of each trawl presented a weak negative trend with time for the Norwegian-Barents Sea and 269 Svalbard regions, which was accounted for within the statistical modelling (Fig. S7, Fig. S8).

273
Spatially explicit environmental co-variates were collated from three sources: Copernicus Marine 274 Service database, Bio-Oracle and MARSPEC database (53-55). Some of the environmental variables 275 were available as annual estimates (e.g., sea surface temperature) whereas others were available as 276 long-term averages (e.g., bottom nitrate concentration) ( Table S3).

278
Variables from the Copernicus Marine Service were available for each year from 1993 to 2019. These

285
Variables from Bio-Oracle and MARSPEC were downloaded using the "sdmpredictors" package from R 286 software (56).

288
A distance to coast layer was created using the "raster" package, also in R (57). All the environmental

294
The change in SBT with time was calculated as the annual mean across each study area (Fig. 2,

295
Supplementary information Methods). SBT thus reflects the likely average temperature individual fish 296 may have experienced over time rather than temperature when sampled.

300
Species richness is typically measured as alpha (local) and gamma (regional) diversity, while the extent 301 of differentiation of communities along habitat (spatial and/or temporal) gradients, is beta diversity (35, 302 60). The most commonly used measures of beta diversity, such as the Jaccard index, are influenced by 303 species richness, which should be accounted for where species richness significantly varies. We thus 304 report both total beta diversity and a species-richness independent index here called turnover, following 305 (34). The mathematical difference between beta diversity and turnover is called nestedness, and arises 306 where sites with less species are a subset of species from neighbouring sites with more species. Thus,

307
we report four measures of diversity: alpha, gamma, total beta, and turnover.

309
Alpha diversity

311
We assessed alpha diversity in the study area as the annual mean species richness per trawl 312 (species/trawl) using presence data only. Generalized Additive Models (GAMs) were used to explore 313 the changes in species richness with time, correcting for the different swept area by adding it as a logistic 314 explanatory variable. The model construction was done using the package "mgcv" within the R software 315 (61, 62).

317
We used boosted regression trees (BRT) to explore the drivers of these changes, and to make spatial 318 projections of the species richness in the study area. A threshold of 0.9 correlation coefficient was 319 selected due to the robustness of the procedure to collinearity (63) (Supplementary information 320 Methods). We restricted the modelling period to the years 1993 -2019 due to the restricted availability 321 of contemporaneous environmental data in the Copernicus Marine repository. The data were randomly 322 divided for model assessment: 75% was used for calibration, and 25% for model validation

325
To study the effect of each environmental variable on species richness, we built partial dependence 326 plots using the R package "devEMF" and basic R from the BRT model. These showed the mean species 327 richness predicted across each environmental variable gradient, while all the other variables were kept 328 at their means (61, 64).

330
To explore the spatial changes in alpha diversity over time, we statistically predicted the geographic 331 distribution of species richness in each year from 2017-2019 and 1994-1996 from the BRT model, and 332 calculated the difference between the mean of each period. We used the "predict" function of the "dismo" 333 package, and the package "raster" from R (57,61,63,65,66). Maps were created with QGIS(59).

335
Beta diversity

336
Beta diversity and its turnover component were calculated using the mean pairwise Jaccard dissimilarity 337 index (ßJ) from presence-only data, which can be divided into two components: species replacement

345
Here we report total beta diversity (ßJ) and the relative contribution of species replacement, to beta 346 diversity, which we refer to as turnover:

348
= . ß 349 To account for differing sample sizes between years we conducted a bootstrapping process. That is, we 350 randomly selected a subsample of the smallest sample size over years, calculated its total beta diversity 351 and turnover, and repeated this process 200 times for each year (34). We report the mean of these 352 bootstrapping calculations in Figure 5 and Figure S6 (34). All calculations were carried out using the 353 functions beta.div.comp, from the adespatial package (68).

355
To explore the trend of these changes with time and to account for the negative trend of swept area with 356 time (Fig. S3), we fitted a GAM model with Year and annual mean swept area as explanatory variables 357 to the annual mean pairwise total beta diversity and turnover.

359
Gamma diversity

360
To study changes in gamma diversity, we constructed species accumulation curves (SAC) for each year 361 across our study regions (69,70). To eliminate the bias that a selected order of the locations may have on the overall SAC, we used the function "specaccum", in the "vegan" package, set the parameter 363 method = "random" to randomize the order of site addition, and permutated the process 200 times (70, 364 71). We report the mean species richness per trawl of all random permutations. We then assessed the 365 change in gamma diversity with time, standardized to the minimum common number of sites, fitting a 366 GAM model with Year and annual swept area as the explanatory variables.

368
To test if the results were robust to alternative measures of richness, we used the package "SpadeR" to 369 calculate nine different species richness indices with time, each of them with different statistical 370 assumptions, including Chao indices, incidence based indices, and first and second order jackknife 371 estimators (72, 73).

375
To study which species drove the changes in biodiversity, we fitted generalized additive models to the 376 presence data of each species with smooth parameters for depth, swept area and latitude, and fixed 377 effect of Year, using a binomial distribution. Each species abundance-weighted mean latitude was 378 calculated from the complete area to arrange the species by their Arctic-affinity, and simple linear 379 regression was used to assess the effect of mean latitude on the temporal change in species probability  using Year and swept area as an offset (Table 1). Red line indicates changes in mean Sea Bottom 557 Temperature across the whole area (correlation with mean alpha diversity r = 0.59). 558  the Slope of the effect of Year on species probability of occurrence, using a GAM model including swept 568 area, Latitude and Depth as smooth predictors, and Year as a linear predictor. Species are ordered by 569 their mean latitude of occurrence during the study period. Linear regression on the effect of Year with 570 mean latitude was conducted and plotted, and asterisks indicate significance (p < 0.05). 571 predicted from a GAM with Year and Latitude as smooth predictors, and swept area as a logarithmic 574 predictor. 575 Table 2. Total species richness predicted from GAM with Year as the only explanatory variable at each 576 region (gamma diversity) in the first and last years of data, at minimum common annual trawls. 577 Table 3. Mean total beta diversity and turnover predicted from GAM at first and last year of data at 578 each region, at minimum common annual sites. 579