Successive extreme climatic events lead to immediate, large-scale, and diverse responses from fish in the Arctic.

The warming trend of the Arctic is punctuated by several record-breaking warm years with very low sea ice concentrations. The nature and reversibility of marine ecosystem responses to these multiple extreme climatic events (ECEs) are poorly understood. Here, we investigate the ecological signatures of three successive bottom temperature maxima concomitant with surface ECEs between 2004 and 2017 in the Barents Sea across spatial and organizational scales. We observed community-level redistributions of fish concurrent with ECEs at the scale of the whole Barents Sea. Three groups, characterized by different sets of traits describing their capacity to cope with short-term perturbations, reacted with different timing and intensity to each ECE. Arctic species co-occurred more frequently with large predators and incoming boreal taxa during ECEs, potentially affecting food web structures and functional diversity, accelerating the impacts of long-term climate change. On the species level, responses were highly diversified, with different ECEs impacting different species, and species responses (expansion, geographical shift) varying from one ECE to another, despite the environmental perturbations being similar. Past ECEs impacts, with potential legacy effects, lagged responses, thresholds, and interactions with the underlying warming pressure, could constantly set up new initial conditions that drive the unique ecological signature of each ECE. These results highlight the complexity of ecological reactions to multiple ECEs and give prominence to several sources of process uncertainty in the predictions of climate change impact and risk for ecosystem management. Long-term monitoring and studies to characterize the vertical extent of each ECE are necessary to statistically link demersal species and environmental spatial-temporal patterns. In the future, regular monitoring will be crucial to detect early signals of change and understand the determinism of ECEs, but we need to adapt our models and management to better integrate risk and stochasticity from the complex impacts of global change.

INTRODUCTION

Marine extreme climatic events (ECEs), such as marine heatwaves, are increasing in intensity and frequency with climate change (Frölicher et al., 2018; Laufkötter et al., 2020; Oliver et al., 2018; Perkins‐Kirkpatrick & Lewis, 2020). Recent reviews on the subject call for an effort toward better understanding of ECEs and their ecosystem consequences to improve their predictability, and support ecosystem‐based management to maintain marine systems resilience (Holbrook et al., 2020; Trebilco et al., 2022; van de Pol et al., 2017). These environmental “pulse” perturbations have often severe ecological impacts (Maxwell et al., 2019; Ummenhofer & Meehl, 2017), altering ecosystem structure and function (Huntington et al., 2020), and threatening the provision of ecosystem services (Cheung et al., 2021; Mills et al., 2013; Smale et al., 2019). The ecological effects of ECEs may accelerate climate‐driven trends or involve additional impacts, such as population collapses, decline of habitat‐forming species (Babcock et al., 2019; Wernberg et al., 2013), or increased sensitivity to multiple pressures (Collins et al., 2017). Two main challenges to understanding ECEs impacts on marine ecosystems are the diversity of ecological responses (van de Pol et al., 2017), challenging a mechanistic understanding of species responses, and the limited scope of empirical evidence. As marine ecosystems are increasingly exposed to successive ECEs over time, it is also critical to identify long‐term and potential cumulative impacts of successive ECEs (Bailey & Pol, 2016; Hughes et al., 2019; van de Pol et al., 2017). Evidence of cumulative impacts of successive ECEs on marine ecosystems is mainly confined to the studies of multiple bleaching events on coral reefs and experimental research on intertidal communities. Documented effects include long‐lasting impaired reproduction (Johnston et al., 2020), increased sensitivity to subsequent bleaching events (Dalton et al., 2020), and loss of protection mechanisms caused by acute stress (Ainsworth et al., 2016). The pace and intensity of successive extreme events are also important as manipulative studies have shown that single ECEs could have more severe impacts than multiple milder ECEs (Sanz‐Lázaro, 2016), but also that rapid succession of ECEs could trigger cumulative impacts (Hughes et al., 2019). However, little is known about the consequences of successive ECEs for other ecosystems and for mobile species.

Over the last two decades of rapid warming, the Arctic has faced longer, more frequent, and intense marine heatwaves. These ECEs were associated with reductions of the first‐year sea ice, changes in ocean stratification, higher sea surface temperatures (Hu et al., 2020), and extremely low sea ice extent (SIE) in 2007, 2012, and 2016 (Parkinson & Comiso, 2013; Petty et al., 2018). Arctic ecosystems are known to respond strongly to the long‐term warming trend, for example, with poleward displacements of fish species (Campana et al., 2020; Eisner et al., 2020; Fossheim et al., 2015), and it is thus likely that they respond to punctual heat stress such as MHWs. Species distribution shifts, changes in productivity regimes, and mass mortality events were indeed observed in association with MHWs in the Pacific Arctic in 2017 (Huntington et al., 2020). However, limited monitoring data challenge more detailed investigations on how arctic marine species and communities respond to successive ECEs as pulses on top of gradual climate trends.

The Barents Sea, a high‐latitude shelf‐sea and gateway from the Atlantic to the Arctic Ocean, is strongly affected by warming (Ingvaldsen et al., 2021; Screen & Simmonds, 2010) leading to an Atlantification of its northern Arctic part with increasing ocean heat content, loss of sea ice, declining sea ice inflows, and weakening ocean stratification (Lind et al., 2018). In this ecosystem, boreal and Arctic species meet at the northern and southern range edges of their distributions, respectively, along a zoogeographic transition. The rapid redistributions documented in the region lead to considerable alterations in fish communities (Fossheim et al., 2015), with associated changes in functional structure and food web properties, as small, benthivores fish in the Arctic are progressively replaced by larger, long‐lived piscivorous boreal fish (Frainer et al., 2017; Kortsch et al., 2015; Pecuchet et al., 2020). Since the middle of the 2000s, the region has experienced years of extreme ocean heat content on top of the warming trends, with up to 6 x standard deviations above the long‐term 1970–1999 average (Lind et al., 2018). Moreover, recent and future gradual expansion of Atlantic water inflow is a major driver of interannual variability in sea ice in the Barents Sea (Dörr et al., 2021; Sandø et al., 2014).

Here, we assessed Barents Sea fish community responses to three consecutive ECEs using observations from the annual Barents Sea ecosystem survey from 2004 to 2017 (Eriksen et al., 2018). Because ECEs have been shown to cause isotherm displacements (Jacox et al., 2020) to which marine species respond quickly (Burrows et al., 2011), we expected each ECE to trigger immediate, behaviorally mediated redistributions accelerating the long‐term effects of climate change. Furthermore, as the fish responses to the gradual warming of the Barents Sea are linked to their traits, with, for example, large, piscivore species with high mobility (Frainer et al., 2017, 2021), we also expected species responses to ECEs to depend on their traits. Our objective was, therefore, to characterize fish responses to ECEs at three different organizational levels (the community, the functional group, and the individual species), with a focus on species behavioral responses (longitudinal and latitudinal displacements, geographical expansions or contractions, and subsequent changes in density), and associated changes in environmental niche. We used these results to discuss how the interplay between climate trends and successive ECEs act on the Barents Sea fish community.

MATERIALS AND METHODS

STUDY AREA

The Barents Sea is a shelf‐sea at the doorstep to the Arctic Ocean. Warm and saline Atlantic Water enters in the southwest and transits through the southern Barents Sea, while the northern part is dominated by colder and fresher Arctic Water (Loeng, 1991). The northern part is covered by sea ice during winter, and the extent varies with sea ice inflows from the northeast (Lind et al., 2018) and the Atlantic Water inflow (Årthun et al., 2012). The sea ice adds freshwater to the upper ocean when it melts in summer, and is the primary freshwater source that maintains the ocean stratification in the northern Barents Sea (Lind et al., 2018). The region is currently transitioning toward Atlantic Ocean climate, with increased salinity and less sea ice (Ingvaldsen et al., 2021; Lind et al., 2018), a process known as Atlantification. It causes major changes in the ecosystem in the Barents Sea, including northward displacement of fish communities (Fossheim et al., 2015), increased primary production (Dalpadado et al., 2020), incoming temperate species such as coccolithophid microalgae Emiliana huxleyi (Oziel et al., 2020), and deteriorating conditions for sympagic fauna in the arctic food web, including marine mammals (Bogstad et al., 2015; Dalpadado et al., 2020).

Data collection

Survey description

Data on fish densities and distributions were extracted from the joint Norwegian–Russian (Norwegian Institute of Marine Research (IMR) and the Knipovich Polar Research Institute of Marine Fisheries and Oceanography (PINRO, since 2019—Polar Branch of Russian Federal Research Institute of Fisheries and Oceanography) ecosystem surveys designed to monitor the entire Barents Sea (Eriksen et al., 2018) each year during summer–autumn between 2004 and 2017. The survey starts in the south and reaching the northernmost areas during the short window of time when it is ice free. Several research vessels are mobilized, so that the central and southwestern areas are covered in 2 weeks and the total area is covered in only 1.5 months. An average of 278 bottom trawls per year (range 230–316 trawls) was sampled with a typical distance of 35 nautical miles, on fixed stations across years (Figure S1). We included only stations that were sampled at least during 7 of the 14 years of the survey to ensure sufficient regional representation over time while avoiding a bias in observed interannual variability due to difference in sampling effort (Figure S1). Furthermore, to ensure equal sampling across the study area, we defined a regular grid allowing only one trawling station per grid cell each year, thus removing repeated samples from the same location.

Fish sampling

Demersal fish were captured with a Campelen 1800 bottom trawl (22 mm mesh size at cod end). The towing speed was approximately three knots. Trawling operations that lasted less than 15 min or more than 1 hour were removed from the dataset, as well as all stations in water shallower than 50 m or deeper than 500 m, to remove coastal and slope stations that harbor different ecosystems than the Barents Sea shelf. On board, the fish caught were identified to the highest possible taxonomic resolution, counted, and weighed. Counts and biomasses were standardized by unit area based on trawling distance and a mean trawl opening of 25 m. Only the taxa present in ≥5% of all sampled stations were kept (33 species/species groups).

Environmental sampling

Ocean temperature and salinity observations were obtained from Conductivity–Temperature–‐Depth (CTD) profiles obtained at each trawling station during the Barents Sea ecosystem survey between 2004 and 2017. We used the temperature and salinity measurements at 10 and 50 m (to define oceanographic domains, Figure S2, see below) and at the deepest level observed. Depth (m) was extracted from NOAA bathymetry raster for the region (Jakobsson et al., 2012). The number of days with sea ice cover in each grid cell was estimated from integrating the number of days with more than 15% daily sea ice concentration from the National Snow and Ice Data Centre (Cavalieri et al., 1996). Chlorophyll a concentrations (mg/m3) were estimated at the sample location through kriging of the NASA ocean color dataset (NASA OBPG, 2018). These data were used at the station level to describe environmental conditions in the Barents Sea and relative changes in environmental conditions from year to year or in comparison with the climatic trend (Table 1a; Figures S2, S3, S4 and S9).

TABLE 1

(a) Environmental variables, their short name and use in the different analyses and exploration of the study. (b) Variables of species capacity to respond to change in their environment, their short name and use in the different analyses and exploration of the study. (c) Variables of species response to change in their environment, their short name and use in the different analyses and exploration of the study

Variable	Short name	Plot or analysis
(a)
Temperature (°C, 10 m, and bottom)	T.10 m and T.bottom	Description of environmental conditions in the Barents Sea (Figures S3 and S8); QGAM analysis for species potential niche (see below)
Salinity (p.s.u., 10 m, and bottom)	S.10 m and S.bottom	Description of environmental conditions in the Barents Sea (Figures S3 and S8); QGAM analysis for species potential niche (see below)
Temperature (°C) and salinity (p.s.u.) at 50 m	T.50 m and S.50 m	Definition of Arctic and Atlantic oceanographic domains
Depth (m)	Depth	QGAM analysis for species potential niche (see below)
Number of days with ice cover	Ice	Description of environmental conditions in the Barents Sea (Figures S3 and S8); QGAM analysis for species potential niche (see below)
Chlorophyll a (mg/m3)	chla	Description of environmental conditions in the Barents Sea (Figures S3 and S8); QGAM analysis for species potential niche (see below)
Heat content in the upper water column (MJ/m², 0–100 m)	—	Detection of ECEs
Freshwater content in the upper water column (m/m²,0–100 m)	—	Detection of ECEs
Sea ice extent (millions km²)	—	Detection of ECEs
(b)
Species preferendum for surface and bottom salinity, chlorophyll a, days of ice, surface temperature	Mode_S.10m, Mode_S.bottom, Mode_chla, Mode_ice, Mode_T.10m	Correlation analysis before PCA. Not included in the PCA because of high significant correlations
Species tolerance for surface and bottom salinity, chlorophyll a, days of ice	Range_S.10m, Range_S.bottom, Range_chla, Range_ice	Correlation analysis before PCA. Not included in the PCA because of high significant correlations
Species preferendum for bottom temperature and depth	Mode_T.bottom, Mode_depth	Correlation analysis and PCA
Species tolerance for bottom and surface temperature and depth	Range_T.10m, Range_T.bottom, Range_depth	Correlation analysis and PCA
Maximum length (cm)	Length.max	Correlation analysis and PCA
Longevity (years)	Age.max	Correlation analysis and PCA
Fecundity (number of batches per female per year)	fecundity	Correlation analysis and supplementary variable in PCA
Offspring.size (mm)	Offspring.size	Correlation analysis and supplementary variable in PCA
Trophic level	tl	Correlation analysis and PCA
Feeding mode	Generalists, specialists, planktivorous, piscivorous, benthivorous	Correlation analysis and supplementary variable in PCA
Body shape	Fusiform, elongated, eel‐like, flat	Correlation analysis and supplementary variable in PCA
Fin shape	Truncated, rounded, forked, pointed	Correlation analysis and supplementary variable in PCA
(c)
Mean densities (ind/m²)	Mean density	Description of species response (Figures 2,4; Figure S4)
Mode of distribution along the longitude (°E)	Mode of longitude	Description of species response (Figures 2,4; Figure S4)
Mode of distribution along the latitude (°N)	Mode of latitude	Description of species response (Figures 2,4; Figure S4)
Geographical extent (number of cells where the species has been observed)	Geographical extent	Description of species response (Figures 2,4; Figure S4)

Long time series for heat and freshwater content in the upper water column (0–100 m) were updated from previous studies (Lind et al., 2016, 2018), from long‐term CTD sampling (end of summer, 1970–2019) in a 85,000 km² area (hereafter referred to as ‘polygon’) representative for the northern Barents Sea, comprised between 18–44°E and 77–79°N (Figure 1). The mean sea ice extent for the whole Barents Sea in June and September was estimated from monthly sea ice concentration fields from NSIDC (70–80.5°N, 22–60°E, Figure 2a) using 15% concentration as a cutoff. We included both the June sea ice extent and the sea ice minimum in September because early retreat of sea ice in early summer increases the potential for ocean heat uptake and increases light conditions supporting phytoplankton blooms and primary production. We chose to compare sea ice extent during the study period (2004–2017) to a reference period, 1999–2009, which was with good coverage in the northern Barents Sea, and represents the start of the warmer period, which is suitable for showing anomalies in the warm period (Figure 1; Figure S3). As the data from the whole Barents Sea ecosystem survey covered only 14 years, we used these longer time series data to detect ECE (see below and Table 1a).

FIGURE 1

ECE in 2016 in the Barents Sea. Ocean heat content anomaly in the upper 100 m in August–September compared to the reference period 1999–2009 (shading) and June sea ice extent in 1979–1999, 1999–2009, and 2016 (blue, orange, and red curves, respectively). The black outline shows the area where the time series of the heat and freshwater content were calculated for the upper 100 m of the water column, shown in Figure 2 and corresponding with Lind et al. (2018)

FIGURE 2

Identification of Barents Sea ECEs and concomitant ecosystem response. Left panels: Environmental conditions in the Barents Sea. (a) Sea ice extent (SIE, millions km²), (b) ocean freshwater content (m/m²), (c) ocean heat content (MJ/m²) in the 100 first meters (Lind et al., 2018), and (d) average bottom temperatures in the whole study area. In A) line type indicates the month: June (solid) or September (dashed). Red‐dotted lines show the 5th and 95th percentile relative to the variability around the nonlinear average of the whole period (see Methods). Red dots indicate years with extreme conditions (<5th or >95th percentile). Red‐shaded ribbons show the confidence interval (95%) around the mean. Dotted vertical lines show peaks in bottom temperature studied in this paper. Right panels: Species‐level (gray lines) and community‐level (black line) standardized (e) mean density (original units: ind/km²), (f) geographical extent (original units: number of cells in which the species has been found) and mode of their spatial distribution along (g) the longitude (original units: °E) and (h) the latitude (original units: °N)

Functional traits

To investigate how fish species responses were associated with their functional traits, we characterized each of the species using traits collected in Beukhof et al. (2019). We selected traits that were available for all species that could inform their responses to changes in habitat and traits with strong implications for community structure. The final list was composed of life‐history traits such as maximum length and age (length.max and age.max), fecundity, and offspring size, indicators of plasticity such as trophic level (tl) and feeding mode (e.g., generalists, specialists, planktivorous, piscivorous, benthivorous) and indicators of mobility: body shape, fin shape (Table 1b; Table S1). We used trait information obtained from the Barents Sea only, except for Icelus spp., for which information was taken from the East Bering Sea (as no data were recorded for the Barents Sea, Beukhof et al. (2019)).

Data preprocessing

Environmental time series per water domains

The polar front that separates the Arctic and Atlantic domains in the Barents Sea is topographically anchored and highly stationary in the western Barents Sea, but more variable between years in the eastern Barents Sea (Loeng, 1991). To determine which water‐mass domain each sampling station belonged to each year we used the salinity at 50 m depth from the CTD profile, which is the typical depth of the core of the Arctic water mass (Lind et al., 2016). Even if the temperature of the Arctic water mass has increased substantially during the 2000s, it is still distinguishable from Atlantic Water by salinity. Therefore, stations with salinity at 50 m < 34.7 (in practical salinity units) and above 72°N latitude were assigned to the Arctic domain, while stations with salinity >34.9 (p.s.u) was assigned to the Atlantic domain (Lind & Ingvaldsen, 2012; Loeng, 1991; Pfirman et al., 1994). If salinity was between those two thresholds, the station was assigned as intermediate water mass domain (Figures S2 and S4).

Species response time series

Mean densities (abundance/m²) across the Barents Sea were calculated for each taxon each year. All average densities were then scaled (calculated as (value‐mean)/standard deviation, Table 1c; Figure S5) per species over their entire time series so that every taxon has the same weight in the community level analysis. The series were further averaged across taxa per year to highlight the peaks in abundance at the community scale. Spatial redistributions were represented by species displacements (the mode of distribution along the longitude and latitude, Table 1c), and expansions or contractions of their geographical extent (the proportion of sampled grid cells in which the species was present each year, Table 1c).

Environmental niche descriptors

To understand if species could cope with short‐term changes in their environment, we assessed the environmental preferences for each taxon (potential niche, assessed from all sampled years), by applying quantile regression on generalized additive models (QGAM, Cade et al., 1999) to each species—environmental variable combination, following Husson et al. (2020). When applying this method to fit high quantiles of the response of a species to an abiotic gradient (here we fitted the 99th), we obtain an approximation of the potential niche regarding the considered environmental variable. We extracted two descriptors: (i) the mode, that is, the preferendum, the region of the environmental gradients one can expect to find the highest abundances and (ii) the range of environmental conditions the species can inhabit. Here, we defined the range as the difference between the values of the environmental variable at the 2.5th and 97.5th quantiles of the gradient where the species was present. Each taxon was thus associated with a mode and a range of bottom and surface salinity and temperature, chlorophyll a concentrations, days of ice, and depth for each species.

Analyses

Extreme climatic event definition

Extreme climatic events are commonly identified when a given environmental condition goes beyond a certain threshold (usually the 90th or 95th percentile for heatwaves, Hobday et al., 2018; Sen Gupta et al., 2020). In our case, as we are both investigating extreme highs (heat content) and lows (freshwater content, sea ice extent), we characterized ECEs as years for which environmental conditions went beyond the 95th percentile or below the 5th percentile of a reference period, respectively (red dots in Figure 2a–c). The reference period is usually either chosen to be a fixed baseline or variation of detrended anomalies (e.g., Holbrook et al., 2019; Jacox et al., 2020). As identified in a recent publications about the Barents Sea (Ingvaldsen et al., 2021), the climate change trend is clearer since 2000, which made it difficult to apply these methods with a linear trend or baseline. Indeed, all latest years appeared as ECEs when applying a fixed baseline as reference period, and the climatic trend is not linear when incorporating data from before 2000 (which is necessary to have enough point to calculate percentiles, Figure2a–c). Therefore, we detected ECEs from the variation in anomalies relative to nonlinear trends (red line, Figure2a–c) fitted on each of the long time series available. The nonlinear trend was defined by locally weighed smoothing (calculated by the function “loess” in r package “stats” (R Core Team, 2019)) with an alpha parameter of 0.75 (i.e., size of the neighborhood used to fit a polynomial surface).

Multivariate analyses

The available time series (14 years) were too short for applying ordinal time series analysis (Hardison et al., 2019). To investigate how species with different traits and from different habitats responded to each ECE, we applied a principal component analysis (PCA) followed by a hierarchical clustering on PC scores (HCPC, from the R package FactoMineR (Le et al., 2008; R Core Team, 2019)) to group species based on similarity in traits and their environmental niche. The objective of joining these two datasets was to cover possible descriptors of species capacity to cope with short‐term changes (trophic plasticity, mobility, life‐history traits, tolerances). We included relevant and uncorrelated variables (correlation absolute value <0.65, Figure S6) from both the trait data and the potential niche descriptors. Among the niche descriptors, significant correlations were found between bottom salinity, depth, surface salinity, temperature, and days of ice. As we expected ECEs to impact limits of suitable habitats, we prioritized the most common parameters limiting spatial distributions in the Barents Sea, as identified by Husson et al. (2020): bottom temperature, depth, and surface temperature modes and ranges. However, as mode of surface temperature was correlated to mode of bottom temperature, it was removed from the analysis (Table S2). In the traits dataset, quantitative variables were not strongly correlated.

Quantitative traits and potential niche descriptors were scaled and centered before applying PCA. Qualitative traits (body and fin shape, feeding mode) and traits that would a priori not be relevant to short‐term response mechanisms (offspring size, fecundity) were added to the analysis as supplementary variables, that is, they were not used to construct the principal components but are indicated in the final plot.

RESULTS

Immediate and large‐scale responses of fish to multiple ECEs

Since the 1970s, there has been substantial increase in interannual variability of ocean heat content and June SIE (Figure 2a–b; Figure S7). On top of this increased natural variability, ECEs involving various environmental drivers occurred more frequently since the beginning of the 2000s (Figure 2a–c). During that period, we identified a sequence of three ECEs in 2006, 2013, and 2016, characterized by high ocean heat content at the end of summer (Table 2; Figure 2b). Also, we identified an ECE of low SIE in June 2012 (Table 2; Figure 2a; Figure S3). These anomalies were associated with a decline in freshwater content in the upper water column (Figure 2c) and high temperatures at the bottom (Figure 2d). The magnitude of changes in heat content between an ECE and its preceding year amount to 6‐ to 12‐year gradual change of the current climatic trend in the region (Figure S8).

TABLE 2

Values of ECE in the 2004–2017 period identified in Figure 2a–c

Year	Extreme environmental condition	Value (unit)
2006	High heat content	319.6 (MJ/m² in 0–100 m)
2012	Low sea ice extent in June	49 652 (km²)
2013	High heat content	538.0 (MJ/m² in 0–100 m)
2014	Low heat content	59.6 (MJ/m² in 0–100 m)
2016	High heat content	637.4 (MJ/m² in 0–100 m)
2017	Low heat content	117.5 (MJ/m² in 0–100 m)

Bottom temperatures peaked with the Barents Sea ECEs of 2006, 2012, and 2016 (Figure 2d, dotted vertical lines in all panels), and we observed concomitant strong peaks in fish species mean density, longitudinal and latitudinal position, and/or geographical extent (Figure 2e–h). The 3 years were extreme, and equivalent in terms of magnitude of changes in bottom temperatures and other environmental conditions relative to the preceding years, and to the variability around the climatic trend (Figure 2d; Figure S9). In contrast, the species responses are not equivalent among the peaking years. The first two ECEs (2006 and 2012) were associated with geographic expansions and changes in densities (Figure 2e–f), while the third ECE (2016) was associated with a sudden northwestward shift (Figure 2g–h). This shift came on top of a positive trend in latitudinal displacement across the study period. These average responses to the ECEs were not driven by a few species displaying stronger responses in peak years, but by a collective response across many of the studied taxa (gray lines in Figure 2e–h).

Three functional groups in the Barents Sea

The first two axes of the PCA on species traits and niche descriptors explained 55% of the variability in the data, while the third axis captured 16% of the remaining variability (Figure 3a). The first axis contrasted species according to their traits. PC1 was positively correlated with large, long‐lived, upper trophic level taxa, and with generalist and fecund taxa. The second axis contrasted species found across a large range of conditions (eurytherms, eurybathic) and with elongated body shapes to those found only in more specific habitats (stenothermic, stenobathic relative to conditions sampled in the Barents Sea) and with flat body shapes. The third axis contrasted deep‐sea fish, confined in cooler, more stable bottom temperatures and species with shallower habitats, that are found in a wider range of thermal conditions. Feeding modes and body shapes also varied along this axis, with bathypelagic fish (fusiform, planktivorous) being associated with PC3 lower scores.

FIGURE 3

Hierarchical clusters and distribution of species based on their traits and habitat. (a) Results from principal component analysis grouping species on environmental niche and quantitative traits. Supplementary qualitative traits are indicated with gray point and text while quantitative traits are indicated with a gray‐dashed arrow and gray text. Colored polygons indicate species grouped by cluster analyses on the principal components. (b), (c), and (d) Spatial distribution of species grouped by the cluster analysis

The cluster analysis based on the three principal components identified three clusters coarsely corresponding to their zoogeographic affiliation (Figure 3a–d; Table 3). The first cluster grouped small, short‐lived, low trophic level fish found in restricted depth ranges at cooler bottom temperatures (Tables S1, S2). Most of these species are characterized as Arctic (Fossheim et al., 2015; Table 3), like polar cod (Boreogadus saida) and polar and bigeye sculpins (Cottunculus microps and Triglops nybelini) or as species living in deeper, colder, more stable waters like the deep water spotted barracudina (Arctozenus risso). The second cluster was significantly associated with low values of PC2 (lower ranges in environmental variables) and grouped species that are found mostly in warmer waters, with low ranges (relative to conditions in the Barents Sea) of surface temperature, and no sea ice cover (Table S2). This group includes boreal taxa that have their main habitat in the Norwegian Sea but expand their distributions into the Barents Sea, such as Norway pout (Trisopterus esmarkii) or Norway redfish (Sebastes viviparus, Table 3). The last cluster grouped predominantly large, eurybathic, long‐lived, high trophic level species, with a high fecundity, such as North‐East Atlantic cod (Gadus morhua), saithe (Pollachius virens), and all wolffish, although herring (Clupea harnegus), a small pelagic fish was also included in this group (Table 3; Tables S1, S2).

TABLE 3

List of fish species used in the analyses. English names and biogeographic group following Fossheim et al. (2015), organized by the groups identified by the cluster analyses. A: Arctic, B: Boreal, AB: Arctic‐boreal. * Missing biogeography for benthopelagic/pelagic taxa that were removed from Fossheim et al. (2015)

Cluster	Species	Abbreviation	English name	Biogeography
Arctic‐like	Arctozenus risso	A. risso	Spotted baracudina	*
	Artediellus atlanticus	A. atlanticus	Atlantic hookear sculpin	AB
	Aspidophoroides olrikii	A. olrikii	Northern alligatorfish	A
	Boreogadus saida	B. saida	Polar cod	*
	Cottunculus microps	C. microps	Polar sculpin	A
	Icelus	Icelus	Twohorn/Spatulate sculpin	A
	Leptagonus decagonus	L. decagonus	Atlantic poacher	AB
	Leptoclinus maculatus	L. maculatus	Daubed shanny	AB
	Liparidae	Liparidae	Snail fishes	A
	Lumpenus lampretaeformis	L. lampretaeformis	Snakeblenny	AB
	Mallotus villosus	M. villosus	Capelin	*
	Triglops murrayi	T. murrayi	Moustache sculpin	B
	Triglops nybelini	T. nybelini	Bigeye sculpin	A
	Triglops pingelii	T. pingelii	Ribbed sculpin	AB
Boreal‐like	Argentina silus	A. silus	Greater argentine	B
	Gadiculus argenteus	G. argenteus	Silvery pout	B
	Micromesistius poutassou	M. poutassou	Blue whiting	AB
	Sebastes viviparus	S. viviparus	Norway redfish	B
	Trisopterus esmarkii	T. esmarkii	Norway pout	B
Widespread predators	Amblyraja hyperborea	A. hyperborea	Arctic skate	A
	Amblyraja radiata	A. radiata	Thorny skate	AB
	Anarhichas denticulatus	A. denticulatus	Northern wolffish	AB
	Anarhichas lupus	A. lupus	Atlantic wolffish	AB
	Anarhichas minor	A. minor	Spotted wolffish	AB
	Clupea harengus	C. harengus	Herring	*
	Gadus morhua	G. morhua	North‐East Atlantic cod	AB
	Hippoglossoides platessoides	H. platessoides	Long rough dab	AB
	Melanogrammus aeglefinus	M. aeglefinus	Haddock	B
	Pollachius virens	P. virens	Saithe	B
	Reinhardtius hippoglossoides	R. hippoglossoides	Greenland halibut	A
	Sebastes mentella	S. mentella	Deepwater redfish	AB
	Sebastes norvegicus	S. norvegicus	Golden redfish	AB

Responses to ECEs varied between and within functional groups

The time series of the three clusters’ responses to the ECEs revealed the same general signals of geographical expansion and latitudinal shifts as seen for the whole community (Figure 4). However, the groups reacted to the ECEs with different intensities and timing. The first ECE of 2006 affected more strongly the arctic and boreal groups, as seen by the peaks in density and geographical extent (Figure 4) while widespread predators reacted mainly to the 2012 peak. Geographic expansion peaked in 2012 for all three groups (Figure 4b). Latitudinal shifts differed between the groups, with peaks in northward displacement in 2006–2008 for the boreal species, in 2016 for the widespread predators, while arctic species showed a more gradual displacement with limited peaks. Among species of the same functional group, individual time series of responses to ECEs displayed highly variable trajectories (Figure 4; Figure S5). While patterns summarized across the whole community indicated a return to the initial state after each ECE for all responses except latitudinal displacements, (Figure 2e–g), the cluster analysis revealed decreasing trends in geographical extent for the Arctic cluster (F = 6.0, 12 df, p‐value: .031), contrasted by the increasing trends of both boreal and widespread predators (F = 5.0, 12 df, p‐value: .046 and F = 6.2, 12 DF, p‐value: .029, respectively) (Figure 4).

FIGURE 4

Average response time series of each cluster. (a) Standardized mean density (original units: ind/km²), (b) geographical extent (number of cells in which the species has been found) and mode of their spatial distribution along (c) the longitude and (d) the latitude. Dotted lines indicate minimum and maximum values of each time series

DISCUSSION

Demersal fish respond to peaks in bottom temperature associated with multiple ECEs in the Barents Sea

We detected multiple ECEs in recent years in the heat content of the upper water column in 2006, 2013, and 2016 as well as in June SIE in 2012. These results are confirmed by a recent study focusing specifically on marine heatwaves and using daily data of sea surface temperatures (SST) from 1982 to 2020 identified the same heatwaves in the Barents Sea in 2013 and 2016, negatively correlated with sea ice concentrations (Mohamed et al., 2022). Marine heatwaves in 2006 and 2012 were frequent in the north, but with a much lower cumulative intensity than 2013 and 2016. These years, in particular 2012, were also associated with extreme ice extent and concentrations, relative to previously recorded conditions in the Arctic (Kirchmeier‐Young et al., 2016; Orsolini et al., 2012).

Marine heat waves are usually triggered in mid‐July to early August, and have increased in duration as the end date occurred progressively later, from mid‐August in 1982–2000 versus late September since 2000 across the Arctic (Huang et al., 2021). The Barents Sea ecosystem survey and the long‐term monitoring data used to detect ECEs in our study are thus concomitant with the heatwave season. The observed responses from the fish communities were matching the patterns in bottom temperature maxima (2006–2007, 2012, 2016), which slightly differed from the heatwaves identified from surface. Future research needs to investigate the vertical extent of those heatwaves to explore if and how they impact the bottom temperature and stratification depth (Huang et al., 2021) to better understand how they affect demersal fish habitats. For example, the 2013 heatwave, although particularly intense, might have been limited to the upper water column, and therefore have little effect on bottom fish habitats. The different responses from the whole community to the three ECEs could have many causes. They could be linked to the difference in intensity of the heatwaves, that is, 2016 being more intense than 2006 and 2012, it could have elicited a different response, including, for example, threshold effects. They could also be linked to the different ecological context in which each event occurs, caused for example by other conditions not included in the analysis, such as interspecific interactions. Future studies could update the data used here to include more data points and ECEs (e.g., the 2020 heatwave, Mohamed et al., 2022)), to statistically link the observed temporal and spatial patterns across fish communities with the environmental conditions. This was unfortunately not possible in our study, owing to the low number of data points in the time series.

The spatial imprint of ECEs and functional group tolerances modulate specific responses

Consistent ecological signatures of the ECEs were observed at the community level, yet there were divergences in functional group responses to the different ECEs. This could emerge from different processes. First, stronger responses of arctic and boreal species to the first ECE in terms of density and geographic extent could be explained by their respective tolerances to environmental conditions in the Barents Sea, and their resulting distribution in the region. Indeed, marine heatwaves displace temperature distributions (Jacox et al., 2020) and thus species’ suitable habitats, and behavioral thermoregulation is a well‐known first buffer against (or means to take advantage of) long‐term warming as well as extreme events for mobile taxa (Burrows et al., 2019; Huey & Tewksbury, 2009; Kearney et al., 2009). As a doorstep area between the Atlantic and the Arctic domains and given its strong thermal gradient, the Barents Sea includes the upper (leading) edge and the lower (trailing) edge of boreal and arctic species’ ranges, respectively. Hence, while the widespread predator group may be more tolerant to the perturbed environmental conditions, the ECEs are likely to have changed habitat availability for the arctic and boreal species, triggering species relocations. Our results may seem counterintuitive as one would expect arctic species to retract during warm years, while they rather seem to expand and increase in abundance. Our two main hypotheses are (i) that the potentially early and extended productive period might have benefitted also arctic species, allowing for a better survival, especially during the two first milder heatwaves of 2006 and 2012 (Mohamed et al., 2022), (ii) that the observed increases in density and geographical extent could be linked to arctic species expanding and aggregating in cooler, more suitable habitats, while still being present in lesser abundances at the trailing edge of their distribution, thus artificially inflating their density in the hauls.

Second, differences in timing and intensity of response from the various functional groups could be linked to the various nature and strength of the ECEs in the different oceanographic domains. Mohamed et al. (2022) showed that marine heatwaves tended to be more frequent in the North of Barents Sea, and more intense and longer in the South. The 2006, 2012, and 2016 heatwaves were more frequent but less intense in the Northern Barents Sea than in the South (Mohamed et al., 2022). In our study, the relative change in environmental conditions in those years was also less extreme in the Atlantic than in the Arctic (Figure S9). Hence, ECEs intensity and timing varied spatially, which might have modulated functional group responses at the local scale.

ECEs “consensus years” hide a high diversity of species‐specific responses

At the species scale, each ECE seemingly affected different sets of species, and some extreme spatial shifts occurred also during non‐ECE years. Previous studies have reported similar species‐specific population trajectories punctuated by “consensus years,” synchronized with extreme climatic conditions. During those years, a multitude of taxa with differing behaviors, life histories and adaptive capacities would display extreme population responses (e.g., Lepidoptera and birds in Palmer et al., 2017). In our study, the fact that the ecological signals occurred the same year as the ECEs suggest that these immediate responses were primarily large‐scale redistributions rather than demographic responses. These redistributions are not stronger during years of ECEs but coordinated across most of the species. The strong diversity of individual species trajectories likely reflects a mixture of species‐specific responses to multiple abiotic and biotic drivers operating on different temporal and spatial scales (Boyd et al., 2016; Smale et al., 2017), with different impacts on their life stages depending on sensitivity bottlenecks (Dahlke et al., 2020; Madeira et al., 2020), as well as indirect responses through species interactions. For instance, for polar cod, a key species in the arctic food web, the life cycle is tightly dependent on temperature, ice, or salinity conditions of the Arctic ecosystem (Gjøsæter et al., 2020; Huserbråten et al., 2019). The warmer Barents Sea has led to both declining extent in the region (Figure S5) as well as major changes in their spawning areas (Aune et al., 2021; Gjøsæter et al., 2020). Similarly, recent studies have shown that ECEs could set unfavorable or favorable conditions for reproduction, recruitment, and growth depending on species thermal sensitivities or that of their preys (Rogers et al., 2021; Smith et al., 2019). In contrast, eurytherm species are less sensitive to variations in temperature, and their responses could be driven less by direct than indirect, food web‐mediated effects of ECEs, stock status, and long‐term warming trends. Cod biomass, for example, is tightly linked to long‐term changes in temperature (Årthun et al., 2018; Kjesbu et al., 2014) but does not seem directly affected by the specific ECEs (Figure S5). Eriksen et al. (2017) have also noticed generalist predators like cod were less affected by short‐term climatic fluctuations, and several studies have linked cod biomass increase in the Barents Sea between 2008 and 2013 to reduced variability in prey availability (e.g., krill, Dalpadado et al., 2020; Orlova et al., 2013, 2015) and improved recruitment (Bogstad et al., 2015; Johannesen et al., 2012; Johansen et al., 2013). Thus, responses to ECEs involve numerous species‐specific processes that challenge our ability to provide common mechanistic explanations to observed ecological signatures.

Possible hidden impacts of ECEs

In our study, the ECEs were primarily associated with large‐scale displacement, while variations in density were emergent patterns of species relocation, such as arctic fish aggregating in refuge areas, or boreal species entering the Barents Sea in higher abundances. Our analysis ignores the species occurring in <5% of the total number of hauls, which means that we might not account for rare boreal species entering the Barents Sea during warm years, or arctic and boreal species not properly sampled with the bottom trawl. Therefore, the study includes mainly species that are relatively often found in the hauls and would be less sensitive to conditions in the Barents Sea. It is thus likely that we underestimate the effect of ECEs on fish communities. In Oceania and the Pacific Artic, previous studies on ECE have also identified redistributions of fish following ECEs (Huntington et al., 2020; Oliver et al., 2017; Wernberg et al., 2013). Moreover, the authors observed lower abundances of copepods and epibenthic biomass, high mortalities in mollusks, low‐energy content in forage fish, and seabirds reproductive failures and die‐offs (Arimitsu et al., 2021; Huntington et al., 2020; Piatt et al., 2020). Drastic changes in species density most often affect habitat‐forming species (e.g., Babcock et al., 2019; Filbee‐Dexter et al., 2020; Straub et al., 2019). The milder ecological signatures to ECEs in our study, together with the absence of reports on mass mortalities or other extreme responses in other functional groups in recent literature of the region, could be due to relatively moderate increases in temperature in the Barents Sea relative to the ECEs in, for example, Pacific Arctic. In the Atlantic part, the temperature increase was between 0.5 and 1°C, while in the Arctic part, the ECEs reflected temperature anomalies between 0.5°C and 1.4°C at the surface and around 0.5°C at the bottom. In comparison, Huntington et al. (2020) reported on temperature anomalies of 3°C near bottom in 2017 in the Pacific Arctic. Future research should further investigate the responses of other functional groups to the ECEs. The extensive monitoring program for the Barents Sea includes annual surveys of both plankton, benthos, marine mammals, and seabirds, and hence offer good opportunities for such studies.

Single ECEs can also have legacy effects, for example, causing delayed effects through species life cycles or interspecific interactions (Johnston et al., 2020; Nowicki et al., 2019). A recent study in the Gulf of Alaska demonstrated that changes across trophic levels lasted up to at least 5 years, following the 2014–2016 “Blob” (Suryan et al., 2021). Also, impacts of the 2017 ECE in the Pacific Arctic were sustained over years, potentially causing a regime shift with a lasting transformation of the ecosystem (Huntington et al., 2020). In the Barents Sea, Johannesen et al., 2012) identified a change in the ecosystem's dynamic after 2006–2008, that is, at the time of the first ECE of our study, and recent work suggest that this transformation is still ongoing (Frainer et al., 2021; Ingvaldsen et al., 2021). However, based on the current analyses, we cannot assess the potential role of the ECE in causing these changes.

Press and successive pulse drivers may interact in the Barents Sea

The combination of the long‐term press with the successive pulse perturbations likely impacts the communities through a combination of acute, gradual, and cumulative effects (Harris et al., 2018). The gradual redistribution and establishment of boreal species in the historically Arctic part of the Barents Sea have been suggested to be detrimental to the Arctic residents through increased habitat use, competition, and predation (Kortsch et al., 2015; Pecuchet et al., 2020). In our study, each ECE led to similar effects, although intensified responses over shorter time periods. In addition, the large‐scale redistribution of fish communities increased the co‐occurrences of smaller arctic species with large predators (Figure S10). This finding corroborates the hypothesis that ECEs might accelerate otherwise slow‐paced, gradual shifts caused by climate change (Smith et al., 2019).

The fact that each ECE had a different ecological signature despite the similarity of their environmental signals (at least in the Arctic domain) suggests that the impacts of press and pulse environmental drivers interact in the Barents Sea. Indeed, ecosystem memory of past events and pressures can cause contingencies between successive ECEs (Hughes et al., 2019). Experimental and modeling studies have revealed that multiple successive pulse perturbations occurring in a short amount of time, or intense events, could alter system state and resilience through collapses or legacy effects, and concluded that historical context and past event characteristics (frequency, duration, intensity) could be important to include to better understand and predict how ECEs could lead to regime shifts (Bello et al., 2019; Fabina et al., 2015; Sanz‐Lázaro, 2016). Mechanisms involved could be, for example, increased sensitivity of impacted species to repeated or long‐term increases in temperature (e.g., Dalton et al., 2020), modification of survival and growth rates (Hughes et al., 2019; Ogle et al., 2015; Saha et al., 2020), relaxed or strengthened competition for habitat with species relocations, or species adaptation to new conditions. Indeed, it has been suggested that mild ECEs could facilitate acclimatization (Clarke et al., 2019; Coleman & Wernberg, 2020; DeCarlo et al., 2019; Grant et al., 2017; Sanz‐Lázaro, 2016). However, fish are mobile and can relocate to avoid the cost of acclimatization, which could delay adaptation (Fox et al., 2019). In our study, the extreme arctic conditions (temperature, ice, photoperiod) can also impose some strong constraints on species adaptation (Ljungström et al., 2021). However, to date, little is known about species adaptive capacity in subarctic and arctic seas. Nevertheless, the high response diversity between species within the different functional groups could point to a community‐level adaptive capacity of the system. Limited synchronicity among ecologically similar species may support the maintenance of ecological functions, also through periods with extreme climatic conditions, thus contributing to ecosystem resilience (e.g., Fauchald et al., 2011; Mori et al., 2013).

Response diversity calls for more flexible and dynamic management

Our study highlighted that bottom temperature maxima could be linked to surface ECEs. Although the environmental changes between a peak year and the previous year exhibit similar environmental signatures, a high diversity of responses emanated from a wide range of possible processes operating at different time and spatial scales in the Arctic. This challenges our ability to reach a mechanistic understanding of ECEs consequences at the individual species scale, where most of the management strategies still operate to date. As ECEs increase in frequency, duration, and intensity (Frölicher et al., 2018; Oliver et al., 2019), managing their stochastic, complex, multiscale consequences is a key condition to achieve sustainable exploitation of marine resources (Trebilco et al., 2022). The Arctic is a region with particularly high stakes as global models anticipate increased production and northward shift of important commercial stocks (Bryndum‐Buchholz et al., 2019; Lotze et al., 2019). However, high uncertainties emanating from disagreements between model predictions and the exclusion of important cross‐scale processes question the possibility of a successful Blue Growth in the region (Bryndum‐Buchholz et al., 2020; Bryndum‐Buchholz et al., 2019; Niiranen et al., 2018). In addition, although long‐term trends seem beneficial to most commercial stocks, more frequent, intense, and durable ECEs could exacerbate climate change impacts (Caputi et al., 2017; Cheung & Frölicher, 2020), for example, through affecting vulnerable life stages (Rogers et al., 2021) or energy flows (Arimitsu et al., 2021).

To better understand ECEs drivers and consequences in the Arctic, regular monitoring of Arctic ecosystems, such as the Barents Sea ecosystem survey, is required. In particular, recent studies have suggested that the detection of temperature hotspots, indicators of community changes, and monitoring of early life stages could serve as warning signals of ongoing changes linked to ECEs or regime shifts (Caputi et al., 2016; Pedersen et al., 2020; Rogers et al., 2021). While monitoring data are often too short to address lagged and cascading effects, models can link ECEs’ properties (e.g., intensity, longevity) to diverse species and functional group responses over time. Yet such modeling is challenged by a range of uncertainties, including unknown thresholds and adaptive capacities (Jacox et al., 2020; Maxwell et al., 2019; Nowicki et al., 2019), as well as the high singularities of species responses and ECEs’ ecological impacts. There is thus a strong need to adapt multispecies and ecosystem models to address ECEs adequately by gaining insight on the triggers and mechanisms (Raymond et al., 2020), but also by better integrating stochasticity in the frequency, duration, and intensity of disturbance processes (Fraterrigo et al., 2020). Finally, we need to adapt our management models to fully integrate risk assessments and provide more proactive, flexible, and dynamic strategies in response to rapid unforeseeable changes (Barbeaux et al., 2020; Gladstone‐Gallagher et al., 2019; Holsman et al., 2019; Yuan et al., 2017).

CONFLICT OF INTEREST

The authors declare no competing financial interests.

AUTHOR CONTRIBUTIONS

BH wrote the paper together with SL, MF, MSM, and RP. SL did oceanographic data preparation. BH and HKS did the analyses and HKS wrote the methods for the statistical analyses. LP, RI, and AD brought their expertise and reviewed the manuscript before submission.

Fig S1–S10

Click here for additional data file.

Table S1

Click here for additional data file.