Ocean temperature
The independent bottom water temperature estimates all showed variable but increasing autumn temperatures from 1996 to 2010, followed by a decline through to 2018 in the mid-depth offshore stations and the survey, but less of a decline in the inshore or deep offshore stations (Fig. 1, 2a). Both deep offshore stations warmed by ~ 0.2 °C over the time series. The ~ 1 °C temperature range evident in the southwestern mid-depth offshore FX8 autumn hydrographic series was larger than that of the less variable northern mid-depth offshore SI7 hydrographic station, but neither of the mid-depth offshore stations showed any net warming over the 22-year period. The shallow inshore stations both warmed considerably (1.5–2 °C) over the time series. Although no net increase was apparent in either of the mid-depth offshore hydrographic station autumn time series, 0.25–0.40 °C increases were observed in both the winter and spring measurement time series at both mid-depth hydrographic stations.
Map of 245 autumn groundfish survey core stations (black dots) sampled around Iceland every October between 1996 and 2018 (except 2011). Bottom-temperature profiles were collected independently at hydrographic stations FX3, FX8 and FX9 (red triangles) and SI1, SI7 and SI8 (blue triangles). Depth contours are 200, 500 and 1000 m.
The mean environmental temperature from the autumn survey (Tempe) lay midway between that of the two mid-depth offshore hydrographic stations, over a similar ~ 1 °C temperature range. Based on a fitted regression to the environmental temperature time series, there was a net warming of 0.33 °C across the survey area over the 22-year time period, although the increase was certainly not linear (Fig. 2c). There were large differences in both the temperature range and temperature-at-depth among some of the regions (Fig. 2b). Bottom water temperatures at 200 m depth were 3–4 °C warmer in the south (SW = 7.3 °C; SE = 8.2 °C) than in the north (NW = 5.4 °C; NE = 4.5 °C), and 4–5 °C warmer at a depth of 500 m (NW = 1.4 °C; NE = −0.5 °C; SW = 6.6 °C; SE = 3.9 °C). Although most stations deeper than 500 m in the NE were < 0 °C, all of the stations at that depth in the SW were > 3 °C.
Bottom-temperatures from hydrographic stations (A: left panel) and from the autumn groundfish survey (B: top right and C: bottom right panel). (A) Near-bottom temperatures at hydrographic stations FX3 and SI1 (~ 70 m), FX8 and SI7 (~ 400 m), and FX9 and SI8 (~ 1000 m) between 1996 and 2018 are shown for winter (January–March), spring (May–June) and autumn (October–December). A geometric smooth has been fitted to the annual means. Note different scales on the y-axes. (B) Mean depth-temperature profiles for each region based on bottom temperatures from the autumn survey. Fitted lines are loess regressions. (C) The mean annual environmental temperature time series estimated from a GLM of the autumn survey bottom temperature measurements (see text for details). The linear regression fitted to the annual means is intended only to show the overall rate of warming through the time series.
Survey catches
A total of 5390 tows were fished over the 22-year time series. Station depths ranged from 26 to 1203 m (mean of 377 m) while bottom temperatures ranged between − 2 and 11 °C (mean of 4.85 °C). Of the 7,246,474 fish representing 200 species that were caught, 82 species were caught during at least 19 of the survey years, and thus are the primary subject of this analysis (Suppl. Table 1). The mean standardized annual abundance of these 82 species ranged between 1 and 27,513 fish, with an overall mean across species of 4096; 70% of the species had a mean standardized annual abundance of less than 100.
There was no clear identifier of species that were rare, and thus caught sporadically, versus those that were newly immigrated to the survey area in response to the warming environment. Nor was it possible to exclude the possibility that certain rare species were incorrectly identified in some survey years. However, there were 6 species that were never caught prior to 1998, yet appeared in increasing numbers in subsequent years (Table 1). All of these species were caught in at least 8 survey years, and as many as 19. However, three of the species were deepwater (> 800 m), and thus unlikely to have been subject to temperature increases of > 0.2 °C. Of the remaining three species, Atlantic mackerel (Scomber scombrus), blackbelly rosefish (Helicolenus dactylopterus) and hollowsnout grenadier (Coelorinchus caelorhincus) were all warmwater (TB > 2) and stenothermal (Steno < 3.4), thus making warming temperatures a likely explanation for their recent arrival to the Icelandic fish community. There were no species that were originally present in the survey, but absent after 2015.
Fish thermal habitat
Fish were caught over the entire temperature range of the survey (− 2 to 11 °C). However, not all temperatures were equally represented by fish species. For the 82 species caught in at least 19 years, the mean thermal bias index (TB) was − 0.9 (range of − 6.1 to 3.2). Most species (n = 49) were associated with waters cooler than the overall environmental mean temperature of 4.85 °C, while the remaining 33 species were in waters warmer than the environmental mean. The mean Steno index of 4.0 (range of 0.8–9.2 °C) indicated that most species were relatively stenothermal, and thus intolerant of a broad range of temperatures. Almost all shallow water species (< 300 m) were warmwater (TB > 0), while the mid-depth species (300–800 m) were mainly coldwater; perhaps because of the presence of so many deepsea species in relatively warm southerly waters, all of the deepwater species were classified as coolwater, not coldwater (Fig. 3). Both shallow and deepwater species tended to be relatively stenothermal (Steno < 5); the most eurythermal species were the species from the mid-depths (300–800 m).
Temperature affinity indices of 82 fish species as a function of mean occupied depth (m). The Thermal Bias (TB) index of an individual species may be towards waters that are warmer (> 0) or colder (< 0) than that of the overall environment (i.e. a warm-water (> 0), cool-water (− 3 to 0) or cold-water (− 7 to − 3) species). The Stenothermal Index indicates the range of temperatures which are occupied: narrow (stenothermal) or broad (eurythermal).
The thermal habitat of most fish species could be predicted reasonably well based solely on thermal affinities. A GLM of temperature-at-capture with Species as the only factor explained much of the temperature variance (P < 0.001, r2 = 0.51) (Suppl. Table 2). Addition of Year as a covariate did not improve the explained variance, but all terms were significant (GLM, P < 0.001, r2 = 0.51) (Suppl. Table 2b). The slope of the Year term was statistically significant at 0.0026 (SE = 0.00014), but there is probably limited biological significance of residency in waters which are only 0.06 °C warmer by the end of the time series. Inclusion of an interaction term between Year and Species was significant, but did not improve model fit. A model which replaced the Year term with the annual environmental temperature (Tempe) and Depth (as covariates) produced a better overall model fit (GLM, P < 0.001, r2 = 0.62) (Suppl. Table 2). The parameter estimates for Tempe (0.718 ± 0.003) and Depth (− 0.007 ± 0.00001) suggested that many fish species tolerated an increased temperature as the environmental temperature increased or as they moved to shallower water. Again, the interaction term between Species and Tempe was significant, but did not improve overall model fit.
The influence of changing fish abundance on thermal habitat was tested by including standardized annual abundance (SAy) as a covariate in the above GLM, both as a main effect and with slopes nested within species (P < 0.001, r2 = 0.62) (Suppl. Table 2). All terms were significant, although SAy as a main effect was much less significant (P = 0.05) when the interaction term with species was present. It was easiest to see the influence of SAy on thermal habitat by regressing the standardized residuals of a GLM model with no abundance term on SAy. A total of 19 of the 82 species had significant relationships between the GLM residuals and standardized annual abundance (LM, P < 0.05), of which 13 of the 19 species had negative slopes (indicating a reduced sensitivity to increasing temperature with increased abundance). However, only 9 of the 19 species had mean annual abundances of more than 100 fish. Therefore, if there was an effect of within-species abundance on their thermal habitat, it was not strong.
Temporal shifts in depth distribution
There was no evidence of strong shifts in depth distribution across the time series of the survey. Of the 82 species, 22 species showed significant inter-annual trends in depth: 10 of these were negative (moving shallower) and 12 were positive (moving deeper). There was no significant relationship between the regression slope parameter of these 22 species and either the Steno or the TB indices (LM, P = 0.10, 21 df).
A GLM of fish depth at capture with Species as a factor and Year as a covariate resulted in a Year parameter of − 0.287 ± 0.008 (GLM, P < 0.001, r2 = 0.86), suggesting that there was a net shift to waters that were about 6 m shallower across the 22 survey years (Suppl. Table 3). The interaction term between Species and Year was significant, but did not appreciably improve model fit.
Shifts in spatial distribution
Spatial shifts in species distribution across years were often difficult to detect in distributional maps, due to routine inter-annual variance associated with the groundfish survey (Suppl. Figure 1). In addition, changes in population abundance often confounded any distributional patterns that might have been present. In species such as the grey gurnard (Eutrigla gurnardus), long-term distributional shifts were visible as a northwestwards shift in abundance over the 22-year survey period (Fig. 4a). A time series of the latitudinal centre of mass confirmed a significant northwards shift in distribution at a mean annual rate of 0.014° latitude·year−1 (LM, P < 0.001, 21 df, r2 = 0.64) (Fig. 4b). For most species however, the analyses reported below provided more robust indicators of distribution shifts.
Long-term shifts in species distribution were often difficult to see in abundance-weighted maps, but were visible in species such as the grey gurnard, Eutrigla gurnardus, between 1997 and 2018 in the autumn survey. (A) Survey catch locations in early and recent years, where symbol size is proportional to catch number, exemplified the long-term trend. All years are shown in Suppl. Figure 1; (B) time series of mean weighted latitude (°N), fitted with a linear regression, showing northwards shift across years.
Rose plots better showed the estimated direction and distance of net movement of each species. Of the 82 species analyzed, 41 of the species showed significant inter-annual trends in either latitude or longitude across 22 years. Almost all 41 species showed directed movement to the north and northwest (Fig. 5). This pattern was especially evident in shallow and mid-depth species, stenothermal species, and warmwater species. Deepwater species were most likely to move in any direction.
Rose plots showing estimated net movement in bearing (direction) and distance (log km) for 41 species with a significant trend in latitude or longitude over the 22-year duration of the autumn groundfish survey. (A) Colours indicate depth: red = shallow (< 300 m), green = mid-depth (300–800 m) and blue = deepwater (> 800 m) species. (B) Colours indicate Steno Index: red = most stenothermal (0 to 2.5); orange = 2.5 to 5.0; green = 5.0 to 7.5; blue = most eurythermal (7.5–9.5). (C) Colours indicate TB Index: blue = coldwater (− 7 to − 3); green = coolwater (− 3 to 0); red = warmwater (0 to 4).
A multivariate GLM of location at capture (latitude, longitude and depth) with Species as a factor, and Year and SAy (species’ standardized annual abundance) as covariates, allowed the effects of time trends and species abundance trends on species’ distributions to be evaluated (GLM, P < 0.001) (Suppl. Table 4). Inclusion of interaction terms between each of the covariates and species allowed for different trends in each species. The model was highly significant (GLM, P < 0.001), explaining 26%, 34% and 86% of the variance for latitude, longitude and depth, respectively (Suppl. Table 4b). Many of the interaction terms were significant, indicating the presence of species-specific trends. However, overall patterns were evident by examining the same model without interaction terms (Suppl. Table 4a). The positive Year slope parameters for latitude (0.011 ± 0.001), longitude (0.035 ± 0.001) and depth (0.139 ± 0.005) indicated that there was an overall distributional shift towards the northwest and deeper waters through the survey time series. Increasing abundance produced a similar northwestern shift (latitude: 0.013 ± 0.001; longitude: 0.327 ± 0.002), but into shallower waters (− 13.6 ± 0.1).
Temperature and abundance effects on distribution
There is no logical reason why the survey year should influence species’ distribution, except through correlation with a direct effect such as temperature or species’ annual abundance. This hypothesis was tested through use of a model similar to that described above, but with the Year term replaced by environmental temperature (Tempe). Thus the model was a multivariate GLM of location at capture (latitude, longitude and depth) with Species as a factor, and Tempe and SAy as covariates, plus the interaction terms between each of the covariates and species. The resulting model was highly significant (GLM, P < 0.001), explaining 25%, 34% and 86% of the variance for latitude, longitude and depth, respectively (Suppl. Table 5). Most of the interaction terms were significant, indicating the presence of species-specific trends in distribution with both temperature and abundance.
Environmental temperature as a covariate in the above model was a significant predictor of the latitude, longitude and depth of almost all species (Suppl. Table 5). As a main effect, the Tempe parameter indicated an overall increase in latitude across species of 0.267 ± 0.003° for every 1° increase in temperature. When combined with the interaction term, which allowed for species-specific movements, 70% of the 82 species were predicted to move northwards as temperatures increased. A Tempe parameter of 6.78° ± 0.29° for longitude indicated an overall shift to the west with warming. Again, most Species by Tempe interaction terms were significant, indicating species-specific shifts in longitude with temperature. A total of 60% of all species would be predicted to move westward with warming water temperature. The Depth parameter, plus many of the interaction terms, were also significant, predicting an overall movement to shallower waters of 29.5 ± 6.4 m for every 1 °C increase in temperature; 60% of the species were predicted to move into shallower water.
Standardized species-specific annual abundance (SAy) was a significant term in the GLM, both as a covariate and in interaction with Species (Suppl. Table 5). However, abundance did not provide a consistent change in latitude, longitude and depth across species. With a doubling of SAy, 52% of species would be predicted to move northwards, 49% would move westwards, and 36% would move deeper.
Spatially-explicit predictions of distribution shifts with warming and abundance
More detailed predictions of distribution shifts in response to temperature increases were possible with a multivariate GLM of the location of each species with Region as a factor, Tempe and SAy as covariates, and interaction terms between each of the covariates and Region. Thus the model allowed for species-specific shifts in latitude, longitude and depth within each region, without the geographical land mass barriers that might be imposed by a non-regional model. The model was fitted separately for each species (Suppl. Table 6). All models were highly significant (GLM, P < 0.001), with 72% of the 82 species showing significant parameter estimates for Tempe as either a main effect or in interaction with Region. Rose plots for these 59 species showed clear distribution shifts in response to a 1 °C increase in temperature (Fig. 6). However, the distribution shifts differed among regions, whereby most species (except those living in deepwater) moved either offshore or along-shelf; only one species moved closer to shore by more than 30 km. Of the 20 species predicted to shift distributions by more than 100 km, 45% were warmwater. Although more coldwater species tended to move long distances in the north, coldwater species were also more prevalent in the north. A similar pattern was observed in the south, where warmwater species were more likely to move long distances to the south, but were also more prevalent in the south. The mean predicted distribution shift over all regions was 38 km (range of 1–326 km), with 7% of the species predicted to shift their centre of mass by more than 100 km.
Regional rose plots showing estimated net movement in bearing (direction) and distance (km) for 59 species with a significant trend in latitude or longitude over the 22-year duration of the autumn groundfish survey. Each vector represents a species colour-coded by their Thermal Bias Index: blue = coldwater; green = coolwater; red = warmwater. Circular contours represent 50-km distances, constrained to a maximum of 150 km.
Although both temperature and abundance were found to influence regional species distribution in the GLM, logic suggests that temperature increases would tend to shift populations northwards, while abundance increases would be unbiased with respect to direction. Our analysis indicated that the mean predicted increase in latitude with a 1 °C increase in temperature for all species and regions with a significant Tempe parameter estimate would be a 0.074° ± 0.031° shift to the north (~ 8.2 km) in conjunction with a shallowing of 126 ± 75 m. Of the 100 species/region combinations which were significant, 62% would be predicted to move northwards. The predicted latitude shifts did not differ significantly among the four regions (P > 0.08), but the two western regions accounted for most of the northward-moving species. The northward prediction differed significantly from a zero net shift (P < 0.05), while the depth prediction did not (P > 0.10).
Species abundance (SAy) was a significant effect in the GLM for 51 of the 82 species (Suppl. Table 6). The mean predicted increase in latitude with a doubling of annual species abundance was 0.036° ± 0.025°, a mean which did not differ significantly from zero (P > 0.10). Overall, 55% of the regional species would be predicted to move northwards, a percentage which does not differ significantly from a random orientation (Chi-squared test, P > 0.10). The difference in regional latitude shifts was not significant (P > 0.10), but virtually all of the predicted northwards movement took place in the northeast region. Most species (61%) were predicted to move to waters which were 15 ± 8 m shallower.
A matched pair comparison for each species/region combination allowed an evaluation of the relative contribution of temperature and abundance to distribution shifts wherever both parameters were statistically significant. Of the 40 species/region combinations, 55% would be predicted to move northwards more after a 1 °C increase in temperature than after a doubling of species abundance. The mean predicted latitudinal shift was 3.4 times greater for the temperature effect than the abundance effect.
Functional predictors of warming-induced distribution shifts
The statistical models described earlier were essential to identify those species sensitive to warming-induced distribution shifts. However, it is more interesting to identify the overall physiological and ecological characteristics that made those species sensitive. A GLM of the predicted distance shifted in response to a 1° increase in water temperature included those biological variables with the potential to influence distribution: Depth (0–300, 300–800 and 800 + m), the Steno and TB indices, regional species’ coverage (Areasr) and abundance (Ar), all in interaction with Region (Suppl. Table 7). All factors except Ar were significant, although the significant effect of Steno was through its interaction with Region. Overall, the strongest influences on the predicted distance shift were Region, Steno, Depth and Areasr. Similar results were obtained if the analysis was restricted only to deepwater or shallow species, or if the predicted latitudinal or depth shift was used as the dependent variable rather than distance. The marginal effects of each of these biological characteristics is shown in Fig. 7. Distance shifted was considerably reduced in deepwater species and in species with broad spatial coverage. Stenothermal species shifted greater distances than eurythermal species. However, the pattern with TB was particularly interesting: coldwater and warmwater species were more likely to shift distributions than those species residing in intermediate temperatures.
Predicted distance moved by each species in four regions around Iceland as a function of species’ ecological characteristics. The 67 species shown were those that showed significant regional temperature responses in the spatially-explicit model described in the text. (A) spatial coverage (number of stations occupied); (B) Steno index; (C) TB index; (D) depth.
Source: Ecology - nature.com
